C  STGF  VERSION 4.5-5.2  22.7.87-22.5.88(MEU'92AVR8WE)/22.4.93-JUL28
C  USES IOPT1=10 AND CHANNELS 11 AND 12 TO PUNCH JJOM OR JJXI INPUT
C  IOPT1=12 IF TERM COUPLING COEFFICIENTS SHALL BE READ AND TRANSFERED
C  PUNCHES REACTANCE MATRICES FOR IMESH=1 AND XI FOR IMESH =2. HES JAN92
C  CALCULATIONS FOR FREE ELECTRONS.
C       + UPDATES (OUTJJM ETC) + CORRECTIONS (REACT,OUTJJ) WE'92AUG-SEP
C  MODIFICATIONS SINCE STG 2                           ...99FEB08/MAY20.
C        INCLUDES IRDEC
C        HANDLES NEW BUTTLE FIT
C        OPENS FILES B00, B01, ....
C        NUMEROV INTEGRATIONS REQUIRED IN STGBB, STGBF
C        NO UNIT 7 WRITES FOR IRAD = 2
C
C
C     IMPLICIT REAL*8(A-H,O-Z)
C
C     PARAMETER (MXTST=( 109* 109- 109)/2)
      PARAMETER (A0=1.0,A1=1.E-9, ZERO=0.0,TINY=-1.E-6)
C
      LOGICAL QDT, WARN,WARNE
C
      COMMON/CQDT/R2ST( 500),QDT,NQ
      COMMON/CEN/ETOT,MXE,NWT,NZ
      COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,LACC
      COMMON/CPOINT/RZERO,RONE,RTWO,H,KP0,KP1,KP2
      COMMON/CDEGEN/NASTD,NASTR,NLEV( 109)
      COMMON/CINPUT/ CF( 500, 500,  20),COEFF(3,  25),ENAT(0: 109),
     5 WMAT( 500,7000),VALUE(7000),RA,BSTO,AZAZ,
     1 LAMAX,LRANG2,LRGL2,MNP2,NAST,NCHAN,NELC,NPTY2,NSPN2,NZED,MORE2,
     2 ISAT( 109),LAT( 109),NCONAT( 109)
      COMMON/CHAN/ECH( 500),EPS( 500),FKNU( 500),CC( 500)
     1  ,RINF( 500),ITARG( 500),LLCH( 500),NCHF,NCHOP,NCHOP1
      COMMON/CMESH/EMAX,EMIN,DEOPEN,DQN,QNMAX,EMESH(9000),IMESH
      COMMON/COMEGA/IE,IRAD0,LOM
      COMMON/CNTRL/IPRINT,IRAD,IPERT,KP2X
      COMMON/CBODE/WBODE(2000)
      COMMON/CDEC/ARAD(5886),ARDEC( 109),IRDEC,IPAR( 109),KCH( 500,3)
      COMMON /CTOP/ TOPS( 500, 500),LRGLAM,LITLAM(5886),NTOP(5886,2),
     * NTCHAN(0: 109),INDM
C
      DIMENSION OMEGA(0:5886), MSLP( 200), LMESH(9000), OMDIP(0:5886,2)
C
      CHARACTER DATEX*9,TIMEX*8, LMAT*4,NAME*3,FMT*11,NUM(0:9)*1
      DATA NUM,TIMEX/'0','1','2','3','4','5','6','7','8','9','NO TIME?'/
C
      COMMON/CWARN/WARN
      WARN=.TRUE.
      WARNE=.TRUE.
C
C  CALL BLOCK DATA AS SUBROUTINE TO AVOID LINKAGE PROBLEMS
CX    CALL BLOCK  -- OUT QUB'90FEB1, WE.
C
      CALL DATE(DATEX)
      CALL TIME(TIMEX)
C
C  INITIALISE INDM
      INDM=0
C
C  I/O UNITS
C  *********
C
C          1 FOR SCRATCH FILE (USED FOR IRAD.LT.2)
C          5 FOR PROGRAM CONTROL
C          6 FOR PRINTED OUTPUT
C          7 FOR VALUES OF OMEGA (USED FOR IRAD.LT.2). FILE OMEGA.
C          8 FOR DIPOLE DATA INDEPENDENT OF SLPI. FILE F00.
C          9 FOR DIPOLE DATA DEPENDENT ON SLPI. FILES F01 ETC.
C         10 FOR INPUT R-MATRIX DATA. FILE H.
C         11 FOR K-MATRICES, JJOM FORMAT, OR XI MATRICES IF IMESH=2
C         12 FOR CASE DESCRIPTION AS REQUIRED BY JJOM
      open(5,file='stgf.inp',status='old')
      open(6,file='stgf.out',status='unknown')
      OPEN(10,FILE='H.DAT',STATUS='OLD',FORM='UNFORMATTED')
C OUT REWIND(10)
C
C
C  WRITE NEWS
C  **********
C
      WRITE(6,6000)
      WRITE(6,6001)
      WRITE(6,6002)
C
C
C  READ DATA FROM UNIT 5
C  *********************
C
C     IPRINT =-3 FOR MINIMUM PRINT
C            = 3 FOR MAXIMUM PRINT
C     IRAD   = 0 FOR NO RADIATIVE DATA ON UNIT 9
C             -1 SAME BUT OMEGA FILE NON-SPARSE (OMS IN NATURAL ORDER)
C            = 1 FOR RADIATIVE DATA ON UNIT 9
C            = 2 FOR RADIATIVE DATA ON UNIT 9 AND NO COLLISION STRENGTHS
C MOD(IPERT,3)= 0 FOR OMISSION OF LONG-RANGE MULTIPOLE POTENTIALS
C             = 1 FOR THEIR INCLUSION
C             = 2 FOR " AND TOP-UP OF ALL TRANSITIONS, E1: VMB-MJS STYLE
C     /KP2X/ = NUMBER OF POINTS UP TO RTWO FOR IRAD=1 ONLY: MUST AGREE
C                WITH STGB OUTPUT! NORMALLY =0, DEFAULTING TO (AMP)PTS
C     AC     = ACCURACY REQUIRED
C     RONE   = DEBUG PARAMETER (NORMALLY TAKE RONE = 1.)
C              NEGATIVE R1 FOR TRUNCATING F (KAB'94JAN) NOT IMPLEMENTED
C     IMESH  = 1 FOR FIXED INCREMENT IN ENERGIES
C                FOLLOWED BY
C                MXE, E0, EINCR FOR NUMBER OF ENERGIES,
C                               FIRST ENERGY AND INCREMENT
C            = 2 FOR FIXED INCREMENT IN EFFECTIVE QUANTUM
C                NUMBERS FOLLOWED BY
C                DQN, QNMAX, EMIN, EMAX, DEOPEN
C                DEFINED IN SUBROUTINE  MESH
C            = 3 TO READ ENERGY MESH, FOLLOWED BY
C                MXE, QNMAX
C                (EMESH(M),M=1,MXE)
C            = -S TO CHOOSE MESH APROPRIATE FOR CASES WITH TOTAL
C               SPIN 2S+1
C              FOLLOWED BY DATA AS FOR IMESH=2
C     IOPT1 = 1 FOR ALL SLPI CASES
C           = 2 FOR SELECTED SLPI CASES FOLLOWED BY
C               IS, IL, IP FOR CASES SELECTED TERMINATING WITH
C              -1,-1,-1 OR EOF
C           = -1 OR -2 AS ABOVE FOR + CASES BUT FIRST FOLLOWED BY
C             NASTD,(NLEV(N),N=1,NASTD); NOT EFFECTIVE IF IRDEC=1
C KAB/HES   = 10,11 for JAJOM with algebraic recoupl. without/with topup
C             10,11 REDUNDANT, NOW SPECIFYING T OUTPUT. WE/JAT'93APR2-3.
C FEB 91    = 12,13 FOR JAJOM WITH TERMCOUPLING (COEFF. READ IN STGF)
C DEC 97    = 14,15 KAB'S JBIN BINARY FORM; 11, 13 AND 15 SETTING ITOP=2
C MAY92       TOPUP NOT YET WORKING WITH JAJOM. TCC MAY BE SUPPLIED AT
C             END OF STGF INPUT FOR TRANSFER. SEE ALSO SUBR. OUTJJ.
C     IRDEC = 0 FOR NO RADIATIVE DECAYS
C           = 1 FOR GALITIS-AVERAGE CALCULATED ALLOWING FOR
C               RADIATIVE DECAYS (BUT NOT YET IN JJXI)
C     LRGLAM= LARGEST VALUE OF LRGL2 ON H-FILE
C             (READ ONLY FOR IOPT1=1 AND IPERT.GE.2)
C
C  READ IPRINT, IRAD, IPERT, KP2
      READ(5,*)IPRINT,IRAD0,IPERTR,KK
      IRAD=MAX(IRAD0,0)
      IPERT=MOD(IPERTR,3)
C  READ AC AND RONE
      READ(5,*)AC
      READ(5,*)RONE
C  SELECT TYPE OF ENERGY MESH
      READ(5,*)IMESH
      IF (IMESH.EQ.2  .OR.  IMESH.LT.0) THEN
C  CASE OF IMESH = 2 OR = -S
      READ(5,*)DQN
      READ(5,*)QNMAX
      READ(5,*)EMIN
      READ(5,*)EMAX
      READ(5,*)DEOPEN
      READ(5,*)IRDEC
      ELSE
CZHL
C      MULTI-LINE INPUT
C
C      IF(IMESH.EQ.1) THEN
C       READ(5,*)MXE,E0,EINCR
C      ELSE
C       READ(5,*)MXE,QNMAX
C      ENDIF
C      IF(MXE.GT.9000)THEN
C       WRITE(6,699)MXE
C       STOP
C      ENDIF
C  CASE OF IMESH = 1
C
       IF(IMESH.EQ.1) THEN
        MXE=0
        IE=0
        DEOPEN=1.0E+10
   81   READ(5,*)MXE1,E0,EINCR,MOREE
        E=E0-EINCR
        IF(DEOPEN.GT.EINCR) DEOPEN=EINCR
        DO 100 I=1,MXE1
         IE=IE+1
         E=E+EINCR
  100    EMESH(IE)=E
        write(6,'(" mxe1=",i5)') mxe1
        write(6,'(1p5e16.9)') (emesh(i),i=1+mxe,mxe1+mxe)
        MXE=MXE+MXE1
        IF(MXE.GT.9000)THEN
        WRITE(6,699)MXE
        STOP
        ENDIF
        IF(MOREE.GT.0) GO TO 81
        write(6,'(" mxe =",i5)') mxe
       ELSE
C
C  CASE OF IMESH = 3 (AND GREATER -- RUB'94JUN16)
C       DO 101 M=1,MXE; 101 READ(5,*)EMESH(M) -- CHANGED RUB'94OCT25:
C
        READ(5,*)MXE,QNMAX
        IF(MXE.GT.9000)THEN
        WRITE(6,699)MXE
        STOP
        ENDIF
        READ(5,*) (EMESH(M),M=1,MXE)
       ENDIF
CZHL END
      ENDIF
      ISP=0
      IF(IMESH.LT.0) ISP=ABS(IMESH)
C   SELECT IOPT1
      READ(5,*)IOPT1
      NASTD=0
      IF(IOPT1.LT.0) THEN
      IOPT1=-IOPT1
C FOLLOW CODING OF VAL.B IN MAIN, READ2 AND SCALE1
      READ(5,*) NASTD
      IF (NASTD.GT. 109) THEN
      WRITE(6,'(16H TOO MANY STATES,I3)') NASTD
      STOP
      ENDIF
      READ(5,*)(NLEV(N),N=1,NASTD)
      WRITE(6,'(25H TARGET TERMS SQUASHED TO,I3,1H:/(25I3))')
     +NASTD,(NLEV(N),N=1,NASTD)
      ENDIF
      IF (IOPT1.GT.9) THEN
       CALL OUTJJ(-1,-1)
       LMAT='XMAT'
       IF(IMESH.EQ.1.OR.IMESH.EQ.3) LMAT='KMAT'
       FMT = 'FORMATTED'
       IF(IOPT1.GE.14) FMT = 'UNFORMATTED'
       OPEN(11,FILE=LMAT,STATUS='UNKNOWN',FORM=FMT)
       IF(IRAD.NE.0) WRITE(6,*)' IRAD SET ZERO BECAUSE IOPT1.GT.9'
       IRAD=0
      ENDIF
      LRGLAM=0
C  CASE OF IOPT1 = 2
      IF(IOPT1.EQ.2) THEN
        WRITE(6,602)
        KSLP=0
   10    READ(5,*,END=3999)IS,IL,IP
         IF(IPERT.EQ.2)LRGLAM=MAX(IL,LRGLAM)
         IF(IL.GE.0)THEN
          KSLP=KSLP+1
          IF(KSLP.GT. 200)THEN
           WRITE(6,605)KSLP
           STOP
          ENDIF
          ISLP=(ABS(IS)*100+IL)*100+IP
          IF(IS.LT.0)ISLP=-ISLP
          MSLP(KSLP)=ISLP
          WRITE(6,'(52X,I7.3)')ISLP
          IF (ISP.LE.0) GO TO 10
          IF (ISP.EQ.ABS(IS)) GO TO 10
             WRITE(6,'('' SPIN FOR THIS CASE .NE. IMESH'')')
             STOP
         ENDIF
        IF(KSLP.EQ.0)THEN
          WRITE(6,600)
          STOP
        ENDIF
      ELSE
      IF(IPERTR.GE.2)READ(5,*,END=4000)LRGLAM
      ENDIF
 3999 WRITE(6,601) IPRINT,IRAD0,IPERTR,KK, AC,RONE, IMESH,IOPT1
      IF(KK.EQ.0) KK=2000
      KP2X=SIGN(MIN(ABS(KK),2000),KK)
      IF(IMESH.EQ.2 .OR. IMESH.LT.0) WRITE(6,'(10X,8HIRDEC  =,I3)')IRDEC
C
C  RADIATIVE TRANSITION PROBABILITIES (CASE OF IRDEC=1)
C  **********************************
C
      DO 2 I=1, 109
       ARDEC(I)=0.
    2  IPAR(I)=-1
C     USED IN READ2/RAD AND OUTJJ.
      IF(IRDEC.NE.0)THEN
        CALL READ1
        K=0
        DO 4 I=2,NAST
         DO 4 J=1,I-1
          K=K+1
    4     ARAD(K)=-LAT(J)-LAT(I)
        NASTR=NAST
C       AS REQUIRED BY READ2 - RUB'96MAR23; MOREOVER SECURE
        NASTD=0
    6   CALL READ2
C       TEMPORARY IRDEC TEST FACILITY 1996MAR-MAY:
        IF(IRDEC.GT.1) THEN
         IRDEC=IRDEC-1
         GO TO 6
        ENDIF
        CALL RAD
        IF(IRDEC.EQ.0) THEN
          IF(MORE2.NE.0) GOTO 6
            WRITE(6,606)
            STOP
        ENDIF
        REWIND 10
      ENDIF
C
C  START MAIN CALCULATIONS
C  ***********************
C
C-    WRITE(6,'(1H1)') -- 650 FORMAT(1H1,...) IN SCALE1 SHOULD SUFFICE!
C
C     CALL ACSUB:
      ACNUM=(24.*AC)**.166666667
      ACJWBK=(6.*AC)**.2
      ACZP=16.*AC
      LACC=0
      IF(AC.LT.1.E-3)LACC=2
      IF(AC.LT.1.E-4)LACC=4
C
C
C  READ R-MATRIX DATA FOR TARGET
C
      CALL READ1
      CALL SCALE1(IOPT1)
      IF(IMESH.EQ.2 .OR. IMESH.LT.0) CALL MESH
C
C  INITIAL WRITES TO UNIT 8 -- KP2X FOR (AMP)PTS, REF VKL - RUB'94OCT25
C
      IF(IRAD.GT.0)THEN
      H=ACNUM
      KP2=((ABS(KP2X)-1)/4)*4+1
      RTWO=(KP2-1)*H+RZERO
      CALL BODE
      OPEN(8,FILE='F00',STATUS='UNKNOWN',FORM='UNFORMATTED')
      WRITE(8)NZED,NELC
      WRITE(8)NAST,(ENAT(I),I=1,NAST)
      WRITE(8)KP2
      WRITE(8)RZERO,H
      WRITE(8)(WBODE(I),I=1,KP2)
      WRITE(8)IPERT
      WRITE(8)MXE
      WRITE(8)(EMESH(I),I=1,MXE)
      WRITE(8)BSTO
      ENDIF
C
C  START CALCULATION OF OMEGAS
C
      IF (IOPT1.GT.9) GO TO 25
C  INITIALISE OMEGAS TO ZERO
      IF(IRAD.GE.2) GO TO 25
      NOMT=((NAST-1)*NAST)/2
      L=0
      DO 23 K=1,MXE
       IF(IRAD0.LT.0) THEN
        N=0
        DO 21 I=1,NAST
         IF(ENAT(I)+TINY.GT.EMESH(K)) GO TO 22
   21    N=I
   22   NOMT=((N-1)*N)/2
       ENDIF
       L=MAX(L,NOMT)
   23  LMESH(K)=NOMT
      KK=4
      E=A0+A1
      IF(E.GT.A0) KK=8
      OPEN(1,STATUS='SCRATCH',ACCESS='DIRECT',RECL=(L+1)*KK)
C EVT      DECLARE LARGE BLOCK SIZE!
C OLD OPEN(1,FILE='SCF1',STATUS='UNKNOWN',ACCESS='DIRECT',RECL=IRECL)
      DO 24 K=1,MXE
       DO 24 J=0,1
   24   WRITE(1,REC=K*3-J) EMESH(K),(ZERO,I=1,L)
C
C  START LOOP ON SLPI CASES
C  ************************
C
C  READ R-MATRIX DATA FOR NEXT SLPI CASE
C  *************************************
   25 KASE=0
      NRANG2=999
 1000 CALL READ2
      NRANG2=MIN(MNP2/NCHAN,NRANG2)
C     JUST TO HELP IDENTIFICATION -- RUB'96FEB28.
      NKAS=KASE+1
C  FOR IOPT1=2, CHECK WHETHER REQUIRED VALUE OF SLPI HAS BEEN FOUND
      IF(IOPT1.EQ.2)THEN
       ISLP=SIGN((ABS(NSPN2)*100+LRGL2)*100+NPTY2,NSPN2)
       DO 30 K=NKAS,KSLP
        KK=K
        IF(ISLP.EQ.MSLP(K))GOTO 40
   30   CONTINUE
C      CASE NOT FOUND:
       GOTO 2000
C      CASE FOUND:  RE-ORDER STORED UNIT 5 DATA
   40  IF(KK.NE.NKAS) MSLP(KK)=MSLP(NKAS)
      ENDIF
C
C  SCALE R-MATRIX DATA FOR SLPI CASE
C  *********************************
      KK=IPRINT
      IF(IPERTR.EQ.2 .AND. LRGLAM.EQ.LRGL2 .AND. IPRINT.EQ.-2) IPRINT=-1
      CALL SCALE2
      IPRINT=KK
      IF(IPERTR.GE.2) THEN
       IF(LRGL2.GT.LRGLAM) THEN
        WRITE(6,"(' CASE SKIPPED AS BEYOND LRGLAM.')")
        GO TO 2000
       ENDIF
       IF(IPERT.EQ.2) CALL TOP(LRGLAM)
      ELSE
       LRGLAM=MAX(LRGL2,LRGLAM)
C      MERELY FOR IDENTIFICATION ON THE OMEGA OUTPUT FILE.
       IF(IPERT.EQ.0)WRITE(6,641)
      ENDIF
      IF(IPRINT.EQ.0.AND.IRAD.NE.2) WRITE(6,640)
      IF(IRAD.NE.0 .AND. IOPT1.LE.9) THEN
        WRITE(8)NSPN2,LRGL2,NPTY2
        NAME='F'//NUM(NKAS/10)//NUM(NKAS-10*(NKAS/10))
        OPEN(9,FILE=NAME,STATUS='UNKNOWN',FORM='UNFORMATTED')
        WRITE(9)NSPN2,LRGL2,NPTY2
        WRITE(9)MNP2,NCHF
C
        WRITE(9) (ECH(I),I=1,NCHF),
C --- EXTRA OUTPUT NEEDED FOR BETA PARAMETER VERSION OF STGBF
C    KAB MARCH 1991; MODIFIED FOR BP-QN K (CC OUT) -- RUB'96APR1-3WE:
     *         (ITARG(I),I=1,NCHF),   ((KCH(I,J),I=1,NCHF),J=1,3)
        WRITE(9) (VALUE(I),I=1,MNP2),((WMAT(I,J),J=1,MNP2),I=1,NCHF)
        IF(IPERTR.GT.0) WRITE(9)(((CF(I,J,L),I=1,NCHF),J=1,NCHF),L=1,2)
      ENDIF
C
C  START ENERGY LOOP
C  *****************
C
      DO 50 IE=1,MXE
      IPERT=MOD(IPERTR,3)
      ETOT=EMESH(IE)
      IF(WARNE.AND.ETOT.GT..5*VALUE(1))THEN
            WRITE(6,690)ETOT,.5*VALUE(1)
            WARNE=.FALSE.
      ENDIF
      IF(IPRINT.GT.0.AND.IRAD/2.EQ.0) WRITE(6,680)ETOT
C
      CALL POINTS
      IF(NCHOP.EQ.0) THEN
        IF (IOPT1.LT.10) GO TO 50
        IF (.NOT.QDT) GO TO 50
      ENDIF
C
C  CONTINUE WITH CALCULATION OF OMEGAS
      LOM=LMESH(IE)
      CALL REACT(IOPT1)
   50 CONTINUE
C
      IF(IRAD.NE.0)CLOSE(9,STATUS='KEEP')
      KASE=KASE+1
C
C  END OF ENERGY LOOP, GO TO NEXT SLPI CASE
C  ******************
C
C..  MODIFICATION 4.2.87 -- AND 12.5.96 MOD(,), MOVED FROM 1000 (OUTJJ!)
      IPERT=MOD(IPERTR,3)
C..  END MODIFICATION
 2000 IF(IOPT1.NE.2)THEN
        IF(MORE2.NE.0) GOTO 1000
      ELSE
C     CASE OF IOPT1=2
        IF(KASE.EQ.KSLP)GOTO 3000
        IF(MORE2.NE.0)GOTO 1000
C       CASES NOT FOUND ON UNIT 10
        WRITE(6,660) (MSLP(K),K=KASE+1,KSLP)
        WRITE(6,'(10X,30(1H*)/)')
      ENDIF
C
C  END OF SLPI LOOP
C  ****************
C
C
C  PRINT TOTAL OMEGAS
C  ******************
C
 3000 IF(IRAD.EQ.2) GOTO 3200
      IF (IOPT1.GT.9) THEN
      OPEN(12,FILE='JJDAT',STATUS='UNKNOWN')
      OPEN(13,FILE='JOMSPECS',STATUS='UNKNOWN')
      CALL OUTJJ(IOPT1,NRANG2)
      CLOSE(12)
      CLOSE(13)
      IF(IOPT1.LT.14) WRITE(11,'(75X,5HEEEEE)')
      CLOSE(11,STATUS='KEEP')
         ELSE
      OPEN(7,FILE='OMEGA',STATUS='UNKNOWN')
      WRITE(7,'(I5,I7,34X,A9,1X,A8)') NZED,NELC, DATEX,TIMEX
      N=MXE
      IF(IRAD0.LT.0) N=-N
      IF(IPRINT.GT.-2)WRITE(6,650) N
      WRITE(7,'(I5,I7,28H =NAST,MXE  RA,IPERT,LRGLAM=,F6.3,I2,I3,I4,
     * 23H SYMMETRIES  NRANG2.LE.,I2.2/1P,(5E16.6))')
     *   NAST,N, RA,IPERTR,LRGLAM,KASE,NRANG2, (ENAT(I),I=1,NAST)
      NFAIL = 0
      NRARN = 0
      DO 3100 K=1,MXE
       NOMT=LMESH(K)
       READ(1,REC=K*3) (OMEGA(I),I=0,NOMT)
       IF(IPERT.EQ.2) THEN
        DO 3010 J=1,2
 3010    READ(1,REC=K*3-J) (OMDIP(I,J),I=0,NOMT)
        DO 3020 I=1,NOMT
         IF(OMDIP(I,1).LE.0.) THEN
C         DIPOLE CASE OR NO TOPPING (OMDIP=0 IN RESONANCE REGIONS):
          X = OMDIP(I,1)+OMDIP(I,2)
          IF(-X.LE.OMEGA(I)*1.005) GO TO 3020
c         if(K.eq.2) print'(I5,1P,3E12.4,I5)', I,OMEGA(I),X,OMDIP(I,1),K
          NFAIL=NFAIL+1
          OMEGA(I)=X
         ELSE
          X = OMDIP(I,1)-OMDIP(I,2)
          DEL = OMEGA(I)-OMDIP(I,1)
          IF(X+DEL.LT.OMEGA(I)*1.E-5) GO TO 3020
           IF(DEL.LT.X*.9) THEN
            DEL = (X/(X-DEL))*DEL
            X = OMEGA(I)
            OMEGA(I) = DEL + OMDIP(I,1)
            IF(OMEGA(I).LE.X*1.12) GO TO 3020
           ENDIF
           OMEGA(I) = -OMEGA(I)
           NWARN = NWARN+1
         ENDIF
 3020    CONTINUE
c       -start- test *****
c       IF(IPERTR.GT.2) WRITE(7,700) (OMDIP(I,1),I=0,NOMT)
c       -end- test *****
       ENDIF
       IF(IPRINT.GT.-2) WRITE(6,700) (OMEGA(I),I=0,NOMT)
 3100  WRITE(7,700) (OMEGA(I),I=0,NOMT)
      IF(NFAIL.GT.0) print 4004, NFAIL
      IF(NWARN.GT.0) print 4009, NWARN
      CLOSE(7)
      WRITE(6,695)
         ENDIF
C
 3200 IF(IRAD.NE.0) THEN
        CLOSE(8,STATUS='KEEP')
        WRITE(6,'(//10X,22HRADIATIVE FILE WRITTEN/)')
      ENDIF
      STOP
 4000 WRITE(6,670)
      STOP
C
C  FORMATS
C  *******
C
  600 FORMAT(//10X,'READS IOPT1 = 2 FOLLOWED BY TERMINATOR'/)
  601 FORMAT(//5X,'DATA READ FROM UNIT 5'/
     1  10X,'IPRINT =',I3/
     2  10X,'IRAD   =',I3/
     3  10X,'IPERT  =',I3/
     *  10X,'KP2X =',I5/
     4  10X,'AC   =',1P,E9.1/10X,'RONE   =',0P,F7.2/
     5  10X,'IMESH  =',I3/
     6  10X,'IOPT1  =',I3/)
  602 FORMAT(/5X,'VALUES OF  10000*IS+100*IL+IP  READ FOR IOPT1=2')
  605 FORMAT(//10X,30(1H*)//10X,'NUMBER OF SLPI CASES =',I4,
     + ' EXCEEDS MAXIMUM OF (AMPERSAND)SLP =  200'//10X,30(1H*)/)
  606 FORMAT(//5X,'***** DATA ON UNIT 10 INSUFFICIENT'
     +  /5X,'***** FOR CALCULATION OF ALL RADIATIVE PROBABILTIES'/)
  640 FORMAT(//'   ETOT   QDT  IPERT',2X,
     * 'INITIAL AND FINAL TARGET STATES, AND COLLISION STRENGTHS')
  641 FORMAT(/10X,10(1H+),' RUN WITH IPERT = 0 ',10(1H+)/)
  650 FORMAT(///2(1X,79(1H*)/)/20X,'ENERGIES AND TOTAL OMEGAS',
     *  9X,'MAXE =',I5/20X,25(1H*)/)
  660 FORMAT(//10X,30(1H*)/11X,'NO DATA ON R-MATRIX FILE FOR'/
     * 11X,'10000*IS+100*IL+IP ='/(28X,I10))
  670 FORMAT(//5X,5(1H*),' FOR IPERT=2 AND IOPT1=1, VALUE OF',
     +  ' LRGLAM REQUIRED ON UNIT 5 ',5(1H*)/)
  680 FORMAT(/5X,'ETOT =',1P,E14.6/5X,20(1H=)/)
  690 FORMAT(//5X,56(1H*)/5X,1P,'* ETOT =',E10.2,
     * ' LARGER THAN .5*VALUE(1) =',E10.2,' *'/5X,'* RESULTS MAY BE',
     + ' INACCURATE. NO MORE SIMILAR WARNINGS. *'/5X,56(1H*)/)
  695 FORMAT(//10X,45(1H*)/10X,
     + '* COLLISION STRENGTHS WRITTEN TO FILE OMEGA *'/10X,45(1H*)/)
  699 FORMAT(10X,30(1H*)/10X,'MXE = ',I5,
     + ' IS LARGER THAN (AMPERSAND)MSH = 9000'/10X,30(1H*)/)
 6000 FORMAT(/1X,78(1H+)//10X,'STGF-JJ VERSION 5.2, 22.5.88/20.05.99'/
     + 10X,28(1H*)//5X,'MODIFICATIONS SINCE VERSION 1.0'//
     + 10X,'(1) BETTER ACCURACY FOR LARGE L (IN VERSION 1.1)'/
     + 10X,'(2) SUM TO INFINITY (TOP-UP) FOR OPTICALLY ALLOWED'/
     + 14X,'TRANSITIONS (IN VERSION 1.2), FOR ALL SINCE 1998'/
     + 10X,'(3) USE OF QDT TO GIVE GAILITIS AVERAGE (IN VERSION 1.3)'/
     + 10X,'(4) USE OF IRDEC TO GIVE RADIATIVE DECAYS ',
     + '(IN VERSION 2.0)'/10X,'(5) HANDLES NEW BUTTLE FIT',
     + ' (IN VERSION 3.0)'/10X,'(6) OPENS FILES B00, B01, ....',
     + ' (IN VERSION 4.0)'/10X,'(7) DATA ON',
     + ' B FILES REQUIRES NUMEROV INTEGRATIONS IN STGBB, STGBF'/10X,
     * '(8) CORRECTION OF ROUNDING-ERRORS AT THRESHOLDS (VERSION 4.5)'/
     + 10X,'(9) PRINTS DIMENSIONS (VERSION 4.5)'/
     * 9X,'(10) COMBINES NEARLY DEGENERATE TARGET LEVELS (VERSION 5.1)'/
     * 9X,'(11) KP2X (DEFAULT (AMP)PTS) ALLOWS TO LIMIT RTWO AS IN STGB'
     */9X,'(12) PROVIDES K, T, OR CHI MATRICES (THIS VERSION STGFJJ)'/
     */5X,'FOR TOP-UP, READ IPERT = 2'/5X,'IF TOP-UP USED WITH',
     + ' IOPT1 = 1 (ALL SLPI CASES ON H-FILE)'/5X,'LAST RECORD '
     + ,'READ FROM 5 IS LRGLAM (LARGEST VALUE OF LRGL2)'/)
 6001 FORMAT(5X,'FOR IMESH = 2, THE PARAMETERS READ FROM 5 ARE'/
     + 15X,'DQN = INTERVAL IN EFFECTIVE QUANTUM NUMBER'/
     + 15X,'QNMAX = MAXIMUM EFFECTIVE QUANTUM NUMBER'/15X,
     + 'EMIN, EMAX, DEOPEN'/15X,'IRDEC = 0 FOR NO RADIATIVE DECAYS'/21X,
     + '= 1 FOR RADIATIVE DECAYS'//5X,'THE GAILITIS AVERAGE IS ',
     + 'CALCULATED FOR ENERGIES GIVING'/5X,'EFFECTIVE QUANTUM',
     + ' NUMBERS LARGER THAN QNMAX. IN THIS REGION'/
     +  5X,'RADIATIVE DECAYS ARE INCLUDED IF IRDEC = 1'//
     +  5X,'THE STANDARD PRINT OPTION IS IPRINT = 0.'/
     +  5X,'MINIMUM PRINT FOR IPRINT = -3, MAXIMUM FOR IPRINT = +3'//
     +  5X,'IRAD CONTROLS CALCULATION OF COLLISION STRENGTHS AND ',
     + 'RADIATIVE DATA,'/5X,'= 0,-1 FOR COLLISION STRENGTHS ONLY,',
     * ' 1 FOR BOTH, 2 FOR RADIATIVE DATA ONLY'/8X,
     *  '< 0 FOR NON-SPARSE OMEGA FILE, OMEGAS IN NATURAL ORDER'//
     * /5X,'FOR IOPT1 SET NEGATIVE, PARAMETERS READ FROM 5 ARE'/10X,
     * 'NASTD = NO. OF TARGET STATES AFTER COMBINING SELECTED STATES'
     +/10X,'NLEV(N),N=1,NASTD = NO. OF STATES TO BE COMBINED INTO ',
     * 'STATE N'// 5X,'IOPT1= (10,11) OR (12,13) FOR T OR K FORMAT; ',
     * 'CHI MATRICES IF IMESH=2.'//1X,78(1H+)/)
 6002 FORMAT(/10X,'COMPILED FOR DIMENSIONS -'//
     + 15X,'CHANNELS                   (AMP)CHF =  500'/
     + 15X,'TARGET STATES              (AMP)TAR =  109'/
     + 15X,'MULTIPOLES                 (AMP)LMX =   20'/
     + 15X,'SMALL L VALUES             (AMP)LP1 =   25'/
     + 15X,'OUTER-REGION RADIAL POINTS (AMP)PTS = 2000'/
     + 15X,'ENERGY-MESH POINTS         (AMP)MSH = 9000'/
     + 15X,'R-MATRIX POLES             (AMP)MNP = 7000'/
     + 15X,'SYMMETRIES S,L,PI OR J,PI  (AMP)SLP =  200'/
     + 15X,'COEFFICIENTS FOR THETA     (AMP)TET =   50'/
     + 15X,'DEGENERATE CHANNELS        (AMP)DEG =   20'/
     + 15X,'TERMS IN BUTTLE FIT        (AMP)NRG =   48'/)
  700 FORMAT(F11.7,1P,(T12,6E11.3))
 4004 FORMAT(//' FAILURE: ',I4,' DIPOLE ALLOWED COLLISION STRENGTHS HAVE
     * BEEN TAGGED WITH'/ ' =======  A MINUS SIGN, AS TOPPING WOULD MAKE
     * OMEGA SMALLER!'/10X,
C    * AS TOPPING AT LAMDA-1 CHANGES OMEGA BY OVER 2%:'/10X,
     *'INCLUDE MORE PARTIAL WAVES (BUT AVOID THIS AT LOW ENERGIES)!')
 4009 FORMAT(//' WARNING:',I4,' NON DIPOLE COLLISION STRENGTHS HAVE BEEN
     * TAGGED WITH A MINUS SIGN:' / ' ------- TOPPING HAS INCREASED THEM
     * BY MORE THAN 12% (OR IS IT HIGH L .5% NOISE)')
      END
C
C***************************************************************
C
      SUBROUTINE READ1
C
C  READS DATA INDEPENDENT OF SLPI FROM R-MATRIX FILE, FILA
C
C  THE FOLLOWING DATA ARE READ _
C     NZ = NUCLEAR CHARGE
C     NELC = NUMBER OF ELECTRONS IN TARGET
C     NAST = NUMBER OF TARGET STATES
C     LRANG2 = TOTAL NUMBER OF SMALL L VALUES
C     LAMAX = MAXIMUM LAMBDA FOR MULTIPOLE POTENTIALS
C     RA = R-MATRIX RADIUS
C     BSTO = LOGARITHMIC DERIVATIVE
C     FOR I = 1,NAST -
C          ENAT(I) = TARGET ENERGIES
C          LAT(I) = TARGET ORBITAL ANGULAR MOMENTA
C          ISAT(I) = VALUES OF (2*S+1) FOR TARGET STATES
C     FOR I = 1,3 AND L = 1,LRANG2 -
C          COEFF(I,L) = BUTTLE CORRECTION
C
C     IMPLICIT REAL*8(A-H,O-Z)
C
      COMMON/CINPUT/ CF( 500, 500,  20),COEFF(3,  25),ENAT(0: 109),
     5 WMAT( 500,7000),VALUE(7000),RA,BSTO,AZAZ,
     1 LAMAX,LRANG2,LRGL2,MNP2,NAST,NCHAN,NELC,NPTY2,NSPN2,NZ,MORE2,
     2 ISAT( 109),LAT( 109),NCONAT( 109)
C
      READ(10)NELC,NZ,LRANG2,LAMAX,NAST,RA,BSTO
      IF(NAST.GT. 109)THEN
      WRITE(6,610)NAST
      STOP
      ENDIF
      READ(10)(ENAT(I),I=1,NAST)
      READ(10)(LAT(I),I=1,NAST)
      READ(10)(ISAT(I),I=1,NAST)
      IF(LRANG2.GT.  25)THEN
      WRITE(6,620)LRANG2
      STOP
      ENDIF
      READ(10)((COEFF(I,L),I=1,3),L=1,LRANG2)
C
  610 FORMAT(//20X,'TOO MANY TARGET STATES'//
     1 10X,'VALUE READ FOR NAST IS',I4//
     2 10X,'MAXIMUM ALLOWED BY DIMENSIONS IS TAR =  109'/)
  620 FORMAT(//20X,'TOO MANY BUTTLE COEFFICIENTS'//
     1 10X,'VALUE READ FOR LRANG2 IS',I3//
     2 10X,'MAXIMUM VALUE ALLOWED BY DIMENSIONS IS LP1 =   25'/)
C
      RETURN
      END
C
C***************************************************************
      SUBROUTINE READ2
C
C  READS R-MATRIX DATA FOR ONE SLPI CASE, FROM FILA
C
C  THE FOLLOWING DATA ARE READ -
C     LRGL2 = TOTAL ORBITAL ANGULAR MOMENTUM
C     NSPN2 = TOTAL (2*S+1)
C     NPTY2 = TOTAL PARITY
C     NCHAN = NUMBER OF CHANNELS
C     MNP2 = NUMBER OF R-MATRIX POLES
C     MORE2 = ZERO TO TERMINATE SLPI CASES
C     FOR I = 1,NAST -
C          NCONAT(I) = NUMBER OF CHANNELS FOR TARGET STATE I
C     FOR I = 1,NCHAN -
C          L2P(I) = SMALL L FOR CHANNEL I
C     FOR I = 1,NCHAN AND N = 1,NCHAN AND M = 1,LAMAX -
C          CF(I,N,M) = COEFFICIENTS IN MULTIPOLE POTENTIALS
C     FOR I = 1,MNP2 -
C          VALUE(I) = R-MATRIX POLE ENERGIES
C     FOR K = 1,NCHAN AND I = 1,MNP2 -
C          WMAT(K,I) = R-MATRIX AMPLITUDES
C
C     IMPLICIT REAL*8(A-H,O-Z)
C
      COMMON/CINPUT/ CF( 500, 500,  20),COEFF(3,  25),ENAT(0: 109),
     5 WMAT( 500,7000),VALUE(7000),RA,BSTO,AZAZ,
     1 LAMAX,LRANG2,LRGL2,MNP2,NAST,NCHAN,NELC,NPTY2,NSPN2,NZ,MORE2,
     2 ISAT( 109),LAT( 109),NCONAT( 109)
      COMMON/CDEC/ARAD(5886),ARDEC( 109),IRDEC,IPAR( 109),KCH( 500,3)
      COMMON/CDEGEN/NASTD,NASTR,NLEV( 109)
C OUT DIMENSION NCNATR( 109)
      DIMENSION L2P( 500)
C
      READ(10)LRGL2,NSPN2,NPTY2,NCHAN,MNP2,MORE2
C     READ(10)(NCONAT(I),I=1,NAST) -- NASTD.GT.0 FIX  WE  RUB'96FEB29:
      READ(10)(NCONAT(I),I=1,NASTR)
      IF(NCHAN.GT. 500)THEN
      WRITE(6,600)NSPN2,LRGL2,NPTY2,NCHAN
      STOP
      ENDIF
      READ(10)(L2P(I),I=1,NCHAN)
      IF(LAMAX.GT.  20)THEN
      WRITE(6,620)
      STOP
      ENDIF
      READ(10)(((CF(I,N,M),I=1,NCHAN),N=1,NCHAN),M=1,LAMAX)
      IF(MNP2.GT.7000)THEN
      WRITE(6,610)NSPN2,LRGL2,NPTY2,MNP2
      STOP
      ENDIF
      READ(10)(VALUE(I),I=1,MNP2)
      READ(10)((WMAT(K,I),K=1,NCHAN),I=1,MNP2)
C
C     RECONSTRUCT BP CHANNEL QN K (FOR RAD AND TOP) -- RUB'96MAR25
C
      N2=0
      DO 33 I=1,NASTR
C      IF (NSPN2.EQ.0) IPAR(I) = ISAT(I)
       IF(NCONAT(I).EQ.0)GOTO 33
       N1=N2+1
       N2=N2+NCONAT(I)
       IF(IPAR(I).LT.0)IPAR(I)=ABS(NPTY2-L2P(N1)+2*(L2P(N1)/2))
C      AS IDENTIFIERS IN RAD AND OUTJJ;
       DO 32 N=N1,N2
        KCH(N,3)=LAT(I)*2
        KCH(N,1)=-1
        IF(NSPN2.NE.0) GO TO 32
        KCH(N,3)=LAT(I)
        KCH(N,1)=MAX(ABS(L2P(N)*2-LAT(I)),LRGL2-1)
        DO 31 K=N1,N-1
         IF(L2P(N).NE.L2P(K)) GO TO 31
         KCH(N,1)=LRGL2+1
         GO TO 32
   31    CONTINUE
   32   KCH(N,2)=L2P(N)*2
   33  CONTINUE
C
C  IOPT1<0?..MAKE SOME TERMS ENERGETICALLY DEGENERATE BUT KEEP ID
   40 IF (NASTD.NE.NAST) GO TO 44
C I.E.   (NASTD.EQ.0 .OR. IOPT1.GT.9)
        N2=0
        DO 42 I=1,NASTD
         K=0
         N1=N2+1
         N2=NLEV(I)+N2
         DO 41 N=N1,N2
   41     K=K+NCONAT(N)
   42    NCONAT(I)=K
C
  600 FORMAT(//20X,'TOO MANY CHANNELS FOR (IS, IL, IP) = (',
     1  3I3,')'//10X,'VALUE READ FOR NCHAN IS',I4//
     2  10X,'MAXIMUM ALLOWED BY DIMENSIONS IS CHF =  500'/)
  610 FORMAT(//20X,'TOO MANY R-MATRIX STATES FOR (IS, IL, IP) = (',
     * 3I3,')'//10X,'VALUE READ FOR MNP2 IS',I5//
     3 10X,'MAXIMUM ALLOWED BY DIMENSIONS IS MNP = 7000'/)
  620 FORMAT(//20X,'TOO MANY MULTIPOLES'//
     1 10X,'VALUE READ FOR LAMAX IS',I3//
     2 10X,'MAXIMUM ALLOWED BY DIMENSIONS IF LMX =   20'/)
C
   44 RETURN
      END
C
C***************************************************************
C
      SUBROUTINE SCALE1(IOPT1)
C
C   CONVERTS R-MATRIX TARGET DATA TO Z-SCALED FORM
C
C     IMPLICIT REAL*8(A-H,O-Z)
      LOGICAL NEWBUT
C
C  COMMON BLOCKS FROM ASYMPTOTIC ROUTINE
      COMMON/CPOINT/RZERO,RONE,RTWO,H,KP0,KP1,KP2
C     COMMON/CENAT1/ENAT1  --> ENAT(0) FOR USE IN SCALE2
      COMMON/CLOGB/NEWBUT
      COMMON/CBUT/FKN(0:  48),UKN(0:  48)
C
C  COMMON BLOCK FROM SUBROUTINE READ
C  NOTE USE OF NZED IN PLACE OF NZ (FROM /CEN/)
      COMMON/CINPUT/ CF( 500, 500,  20),COEFF(3,  25),ENAT(0: 109),
     5 WMAT( 500,7000),VALUE(7000),RA,BSTO,AZAZ,
     1 LAMAX,LRANG2,LRGL2,MNP2,NAST,NCHAN,NELC,NPTY2,NSPN2,NZED,MORE2,
     2 ISAT( 109),LAT( 109),NCONAT( 109)
C
C  VMB LEVEL FACILITY
      COMMON/CDEGEN/NASTD,NASTR,NLEV( 109)
      DIMENSION AVENGY( 109),INDEX( 109), ENATR( 109)
C
C  NUCLEAR CHARGE
      AZ=MAX(NZED-NELC,1)
      AA=1./AZ
      AZ2=AA*AA
      AZAZ=AZ*AZ
C
C  Z-SCALED TARGET ENERGIES, RELATIVE TO ZERO FOR TARGET GROUND STATE
      ENAT(0)=ENAT(1)
      DO 2 I=1,NAST
       ENAT(I)=(ENAT(I)-ENAT(0))*2.*AZ2
    2  ENATR(I)=ENAT(I)
      NASTR=NAST
C
C  WRITE TARGET PROPERTIES
C     IF (IOPT1.LE.9 .OR. NASTD.EQ.0)
      WRITE(6,650) NZED,NELC, (I,LAT(I),ISAT(I),ENAT(I),I=1,NAST)
C
C     IF(NAST.EQ.1)GOTO 9
      IF (NASTD.EQ.0) GO TO 9
C TEST FOR DEGENERACY AND FORM AVERAGE ENERGY FOR BUNCHED TERMS
      N=0
      DO 3 J=1,NASTD
    3  N=N+NLEV(J)
      IF (N.NE.NAST) THEN
        WRITE(6,660) N,NAST,(NLEV(I),I=1,NASTD)
        STOP
      ENDIF
      N2=0
      DO 5 J=1,NASTD
       N1=N2+1
       N=NLEV(J)
       N2=N+N2
       AA=0.0
       DO 4 I=N1,N2
        INDEX(I)=J
    4   AA=ENAT(I)+AA
    5  AVENGY(J)=AA/N
C     END INSERTIONS FROM VALF
C
      IF(IOPT1.GT.9) THEN
C       MORE CODING FROM VALF (AND HES IOPT1.GT.9)
C       IF THERE ARE BUNCHED TERMSPUT (?)  ENERGIES IN CINPUT
        WRITE(6,665)
        DO 6 I=1,NAST
         J=INDEX(I)
         ENAT(I)=AVENGY(J)-AVENGY(1)
    6    WRITE(6,'(I17,F16.8,F19.6)') I,ENAT(I),ENATR(I)
      ELSE
        N2=0
        DO 7 J=1,NASTD
         N1=N2+1
         N2=N2+NLEV(J)
         IF(N2.GT.N1) WRITE(6,664) N1,N2
    7    INDEX(J)=N2
C       REPLACE ENERGY BY AVERAGE FOR BUNCHED TERMS BUT KEEP ALL ID'S
        WRITE(6,666)
        DO 8 J=1,NASTD
         ENAT(J)=AVENGY(J)-AVENGY(1)
    8    WRITE(6,'(3X,2I14,F21.6)') INDEX(J),J,ENAT(J)
      NAST=NASTD
      ENDIF
C
C
C  RA AND BSTO
    9 RZERO=RA*AZ
      BSTO=BSTO/RZERO
C
C  BUTTLE CORRECTION
      NEWBUT=COEFF(3,1).LE.-10000.
      IF(.NOT.NEWBUT)THEN
        AA=RZERO*AZ2
        DO 80 M=1,3
        AA=AA*AZAZ
        DO 80 L=1,LRANG2
   80   COEFF(M,L)=COEFF(M,L)*AA
      ELSE
        DO 90 L=1,LRANG2
        NBUT=-INT(COEFF(3,L))/10000
        IF(NBUT.GT.  48)THEN
          WRITE(6,699)NBUT
          STOP
        ENDIF
        COEFF(3,L)=NBUT
   90   COEFF(1,L)=RZERO*COEFF(1,L)
C       INITIALISE FKN AND UKN
        G=-1.570796327
        DO 100 I=0,  48
        G=G+3.141592654
        FKN(I)=G
  100   UKN(I)=G*G
      ENDIF
C
  650   FORMAT(//10X,'NUCLEAR CHARGE =',I3,
     * ', NUMBER OF TARGET ELECTRONS =',I3/10X,52(1H*)//
     */20X,'TARGET STATES -'/20X,15(1H*)//
     * 10X,'INDEX    TOTAL L   (2*S+1)     ENERGY SCALED AS RY*Z**2'/
     * 20X,'OR 2*J    AND P'/(3X,3I10,7X,F12.6))
  660 FORMAT(//5X, '***** INCORRECT DATA *****'/5X,'EXPECT',I4,
     *' STATES FROM NLEV DATA BUT NAST IS',I4/5X,'NLEV=',(20I3))
  664 FORMAT(/' STATES',I4,' TO',I4,' ARE COMBINED')
  665 FORMAT(/15X,'EQUIVALENT TARGET STATES -'/15X,26(1H*)//
     C 12X,'INDEX     AVER. ENERGY',7X,'ORIGINAL ENERGY'/)
  666 FORMAT(//15X,'EQUIVALENT TARGET STATES -'/15X,26(1H*)//
     * 5X,'UP TO  OLD INDEX     NEW INDEX',7X,'SCALED ENERGY'/)
  670 FORMAT(//5X,'RZERO =',F11.4,3X,'BSTO =',1P,E12.4)
  699 FORMAT(//' ******* NBUT =',I5,
     + '  LARGER THAN (AMPERSAND)NRG =   48'/)
C
      WRITE(6,670) RZERO,BSTO
C
      RETURN
      END
C
C***************************************************************
C
      SUBROUTINE SCALE2
C
C   CONVERTS R-MATRIX DATA TO Z-SCALED FORM
C
C     IMPLICIT REAL*8(A-H,O-Z)
C
C  COMMON BLOCKS FROM ASYMPTOTIC ROUTINE
      COMMON/CEN/ETOT,MXE,NWT,NZ
      COMMON/CHAN/ECH( 500),EPS( 500),FKNU( 500),CC( 500)
     1  ,RINF( 500),ITARG( 500),LLCH( 500),NCHF,NCHOP,NCHOP1
      COMMON/CPOINT/RZERO,RONE,RTWO,H,KP0,KP1,KP2
      COMMON/CPOT/BW( 500, 500),LAMP( 500, 500)
C
C  COMMON BLOCK FROM SUBROUTINE READ
C  NOTE USE OF NZED IN PLACE OF NZ
      COMMON/CDEC/ARAD(5886),ARDEC( 109),IRDEC,IPAR( 109),KCH( 500,3)
      COMMON/CINPUT/ CF( 500, 500,  20),COEFF(3,  25),ENAT(0: 109),
     5 WMAT( 500,7000),VALUE(7000),RA,BSTO,AZAZ,
     1 LAMAX,LRANG2,LRGL2,MNP2,NAST,NCHAN,NELC,NPTY2,NSPN2,NZED,MORE2,
     2 ISAT( 109),LAT( 109),NCONAT( 109)
      COMMON/CNTRL/IPRINT,IRAD,IPERT,KP2X
C
C  STATISTICAL WEIGHT * 2 (OR * 4 FOR NON-EXCHANGE CASE)
      NWT=(2*LRGL2+1)
      IF(NSPN2) 10,30,20
   10   NWT=-2*NWT
   20   NWT=NWT*NSPN2*2
      GO TO 33
   30 NWT=NWT+1
C
C  NUCLEAR CHARGE
   33 AZ=MAX(NZED-NELC,1)
      AZ1=1./AZ
      AZ2=AZ1*AZ1
      TAZ2=2.*AZ2
      AZHR=1./SQRT(AZ)
      AZAZ=AZ*AZ
C
C
C  CHANNELS
      NCHF=NCHAN
      DO 40 I=1,NCHF
      LL=KCH(I,2)/2
      CC(I)=LL*(LL+1)
   40 LLCH(I)=LL
      I=0
      DO 50 J=1,NAST
      K=NCONAT(J)
      IF(K.EQ.0)GOTO 50
      DO 45 L=1,K
       I=I+1
       ITARG(I)=J
   45  ECH(I)=ENAT(J)
   50 CONTINUE
C
C  WRITE SLPI AND CHANNEL DATA
C
      WRITE(6,600)NSPN2,LRGL2,NPTY2
  600 FORMAT(//1X,77(1H+)//15X,'S L P =',3I3/15X,16(1H*)/)
  670 FORMAT(/12X,'CHANNEL  TARGET    SMALL'
     *       /12X,2('INDEX    '),'   L'/28X,'OR  2K   AND  L')
  680 FORMAT(5X,2(I10,I9))
      IF(IPRINT.GE.-1) THEN
       WRITE(6,670)
       DO 55 I=1,NCHF
        IF(NSPN2.NE.0) THEN
         WRITE(6,680) I,ITARG(I),LLCH(I)
        ELSE
         WRITE(6,680) I,ITARG(I),KCH(I,1),LLCH(I)
        ENDIF
   55  CONTINUE
      ELSE
       WRITE(6,'(17X,I5,9H CHANNELS)') NCHF
      ENDIF
C
C  R-MATRIX
C  VALUE AND WMAT
      IF(IPRINT.GT.2) WRITE(6,"(/' N, VALUE(N) AND WMAT(I,N)'/)")
      DO 70 N=1,MNP2
       DO 60 I=1,NCHF
   60   WMAT(I,N)=WMAT(I,N)*AZHR
       VALUE(N)=(VALUE(N)-ENAT(0))*TAZ2
       IF(IPRINT.GT.2)WRITE(6,640)N,VALUE(N),(WMAT(I,N),I=1,NCHF)
   70  CONTINUE
  640 FORMAT(I5,1P,E12.3/(6X,6E11.3))
C
C  COEFFICIENTS IN POTENTIAL
      IF(IPERT.EQ.0)GOTO 160
  655 FORMAT(/' PERTURBATION P FOR MULTIPOLE POTENTIAL'//
     1 '  - DIPOLE PART'/)
      LMX=LAMAX
      IF(IPRINT.GT.2)WRITE(6,655)
C  LAMP(I,J) AND BW(I,J)
      DO 120 J=1,NCHF
       DO 120 I=1,NCHF
  120   LAMP(I,J)=1
      IF(LMX.EQ.0) GO TO 160
      A1=1./RZERO
      A2=.5*A1
      IF(NCHF.EQ.1)GOTO 140
      DO 130 J=2,NCHF
       DO 130 I=1,J-1
        BU1=-CF(I,J,1)
        IF(BU1.EQ.0.)GOTO 130
        LAMP(I,J)=2
        BW(I,J)=BU1
        IF(IPRINT.GT.2) WRITE(6,'(2I5,1P,E12.3)') I,J, A2*BU1
  130   CONTINUE
  140 IF(LMX.EQ.1) GO TO 160
      IF(IPRINT.GT.2)WRITE(6,"(/'  - QUADRUPOLE PART'/)")
      A2=A2*A1
      DO 150 J=1,NCHF
       DO 150 I=1,J
        BU2=-CF(I,J,2)
        IF(BU2.EQ.0.)GOTO 150
        LAMP(I,J)=3
        BU2=BU2*AZ
        BW(I,J)=BU2
        CF(I,J,2)=AZ*CF(I,J,2)
        IF(IPRINT.GT.2) WRITE(6,'(2I5,1P,E12.3)') I,J, BU2*A2
  150   CONTINUE
C
  160 RETURN
      END
C
C***************************************************************
C
      SUBROUTINE MESH
C
C  CALCULATES ENERGY MESH FOR CASE OF IMESH=2
C  OR IMESH=-S WHERE S=2*TOTAL SPIN +1
C
C
C     IMPLICIT REAL*8(A-H,O-Z)
C
      LOGICAL CUT
      COMMON/CDEGEN/NASTD,NASTR,NLEV( 109)
      COMMON/CEN/ETOT,MXE,NWT,NZ
      COMMON/CINPUT/ CF( 500, 500,  20),COEFF(3,  25),ENATK(0: 109),
     5 WMAT( 500,7000),VALUE(7000),RA,BSTO,AZAZ,
     1 LAMAX,LRANG2,LRGL2,MNP2,KAST,NCHAN,NELC,NPTY2,NSPN2,NZED,MORE2,
     2 ISAT( 109),LAT( 109),NCONAT( 109)
      COMMON/CMESH/EMAX,EMIN,DEOPEN,DQN,QNMAX,EMESH(9000),IMESH
      COMMON/CNTRL/IPRINT,IRAD,IPERT,KP2X
      DIMENSION ENAT( 109)
C
C
C  PARAMETERS READ FOR IMESH=2
C  ***************************
C
C  DQN     INTERVAL FOR EFFECTIVE QUANTUM NUMBER
C
C  QNMAX   LARGEST ALLOWED VALUE OF EFFECTIVE QUANTUM NUMBER
C
C  EMIN    LOWEST VALUE OF ETOT
C
C  EMAX    HIGHEST VALUE OF ETOT
C
C  DEOPEN  INTERVAL IN ETOT FOR ALL CHANNELS OPEN
C
C  IRDEC   RADIATIVE DECAYS INCLUDED FOR IRDEC=1
C
      WRITE(6,600)DQN,QNMAX,EMIN,EMAX,DEOPEN
C
C  CASE OF IMESH =  -TOTAL SPIN  (ADDED 9.06.88)
C
C  NOTE CHANGES TO /CINPUT/
C       NAST REPLACED BY KAST
C       ENAT   ..     .. ENATK
C  FORM NEW VALUE NAST FOR USE ONLY BY SUBROUTINE MESH WHERE
C       NAST=NO. OF STATES WHICH CAN FORM CHANNELS FOR TOTAL
C            SPIN GIVEN BY ABS(IMESH)
C       ENAT=CORRESPONDING ARRAY OF TARGET ENERGIES
C  NOTE SPECIAL TREATMENT IF DEGENERATE LEVELS HAVE BEEN COMBINED
C  IN SUBROUTINE SCALE1
C  IF IMESH IS .GT. 0 NAST, ENAT COPIED FROM KAST, ENATK
C
      IF (IMESH.LT.0) THEN
         ISP=ABS(IMESH)
         ISP1=ISP-1
         ISP2=ISP+1
         IMESH=2
         NAST=0
         IF (NASTD .EQ. 0) THEN
           DO 1   I=1,KAST
            IF (ISAT(I) .LT. ISP1 .OR . ISAT(I) .GT. ISP2) GO TO 1
               NAST=NAST+1
               ENAT(NAST)=ENATK(I)
    1       CONTINUE
         ELSE
           N2=0
           DO 2   I=1,KAST
            N1=N2+1
            N2=NLEV(I)+N2
            DO 3 N=N1,N2
             IF (ISAT(N) .LT. ISP1 .OR. ISAT(N) .GT. ISP2) GO TO 3
               NAST=NAST+1
               ENAT(NAST)=ENATK(I)
               GO TO 2
    3        CONTINUE
    2       CONTINUE
         END IF
      ELSE
         NAST=KAST
         DO 4   I=1,KAST
    4     ENAT(I)=ENATK(I)
      END IF
C     WRITE(6,'('' NAST='',I4/(1X,5F12.6))') NAST,(ENAT(I),I=1,NAST)
C
C  CASE OF EMIN.GE.ENAT(NAST)  (ADDED 9.12.87) GO TO 300 OF 6.5.94
      CUT=.FALSE.
      IE=0
      E=EMIN
      IF(EMAX.LT.E) GO TO 400
      IF(E.GE.ENAT(NAST)) GO TO 300
C
C  FIND NC1, LOWEST LEVEL ABOVE EMIN
      DO 10 N=1,NAST
       NC1=N
       IF(E.LT.ENAT(N)) GO TO 20
   10  CONTINUE
   20 WRITE(6,610)NC1
C
C  INITIALISATIONS
      QNINT=1./DQN
      SE=1./(QNMAX*QNMAX)
      E=MAX(EMIN,ENAT(1)-SE)
      IF(EMAX.LT.E)GOTO 400
C
C  FIND NC2, HIGHEST LEVEL BELOW EMAX
      DO 30 N=NAST,1,-1
       NC2=N
       IF(ENAT(N).LE.EMAX) GO TO 40
   30  CONTINUE
C     GO TO 400
   40 WRITE(6,620)NC2
C
C     RPS-BRAKE - WE'93MAY2:
C     IF(IOPT1.LE.9) GO TO 50
      IF(ENAT(NC1)-SE.GE.EMIN) GO TO 50
      WRITE(6,605)
      EMIN=ENAT(NC1)+.000005
C
C   ENERGIES IN RANGE EMIN TO ENAT(NC2)
C  ************************************
C
   50 IF(IPRINT.GT.0)WRITE(6,625)
      EN=-1
      DO 160 N=NC1,NC2
      ENM=EN
      EN=ENAT(N)
      IF(EN.EQ.ENM)GOTO 160
C
C   RANGE UP TO (EN-SE)
      IF(EN-SE.GT.E)THEN
        F=1./SQRT(EN-E)
        D=(ABS(QNMAX)-F)*.999
        J=D*QNINT+1
        IF(IPRINT.GE.0) WRITE(6,*)
        CUT=IE+J.GE.9000
        D=D/J
        IF(CUT) J=9000-IE-1
        if(iprint.gt.0)print *,' enter DO110:'
        DO 110 K=0,J
         IE=IE+1
         EMESH(IE)=EN-1./(F*F)
         IF(IPRINT.GT.0) WRITE(6,630)IE,EMESH(IE),N,N,F
  110    F=F+D
        E=EN-SE
C TMP  IF(QNMAX.LT.0.) E=EN+.000005
      ENDIF
C
C    RANGE UP TO EN
C...FIND NEXT HIGHER THRESHOLD
      DO 120 M=N+1,NAST
       MM=M
       EM=ENAT(M)
       IF(EM.NE.EN) GO TO 140
  120  CONTINUE
C... NO NEXT THRESHOLD, USE DEOPEN
      F=(ENAT(N)-E)*.999
      J=F/DEOPEN+1
      IF(IPRINT.GT.0)WRITE(6,690)
      CUT=IE+J.GE.9000
      F=F/J
      IF(CUT) J=9000-IE-1
      if(iprint.gt.0)print *,' enter DO130:'
      DO 130 K=0,J
       IE=IE+1
       EMESH(IE)=E
       IF(IPRINT.GT.0) WRITE(6,630)IE,E,N
  130  E=E+F
      GO TO 290
C... USING NEXT HIGHER THRESHOLD
  140 F=1./SQRT(EM-E)
      D=(1./SQRT(EM-EN)-F)*.999
      J=D*QNINT+1
      IF(IPRINT.GT.0)WRITE(6,*)
      CUT=IE+J.GE.9000
      D=D/J
      IF(CUT) J=9000-IE-1
      if(iprint.gt.0)print *,' enter DO150: F,D =',F,D
      DO 150 K=0,J
       IE=IE+1
       EMESH(IE)=EM-1./(F*F)
       IF(IPRINT.GT.0) WRITE(6,630)IE,EMESH(IE),N,MM,F
  150  F=F+D
      IF(IPRINT.GT.0)WRITE(6,*)
      E=EN
  160 CONTINUE
C
C
C  ENERGIES IN RANGE ENAT(NC2) TO EMAX
C  ***********************************
C
C  FIND NEXT HIGHER THRESHOLD
      DO 210 M=NC2+1,NAST
      IF(ENAT(M).EQ.EN) GO TO 210
        EM=ENAT(M)
C
C  RANGE UP TO MIN(EMAX,EM-SE)
  220 IF(EM-SE.GT.E)THEN
        IF(EM-SE.GT.EMAX)THEN
C...  CASE OF EMAX.LT.(EM-SE)
          F=1./SQRT(EM-E)
          D=1./SQRT(EM-EMAX)-F
          J=D*QNINT+1
          IF(IPRINT.GT.0) WRITE(6,*)
          CUT=IE+J.GE.9000
          D=D/J
          IF(CUT) J=9000-IE-1
          if(iprint.gt.0)print *,' enter DO230:'
          DO 230 K=0,J
           IE=IE+1
           EMESH(IE)=EM-1./(F*F)
           IF(IPRINT.GT.0) WRITE(6,630)IE,EMESH(IE),N,M,F
  230      F=F+D
          GO TO 400
        ELSE
C... CASE OF EMAX.GT.(EM-SE)
          F=1./SQRT(EM-E)
          D=(ABS(QNMAX)-F)*.999
          J=D*QNINT+1
          IF(IPRINT.GT.0) WRITE(6,*)
          CUT=IE+J.GE.9000
          D=D/J
          IF(CUT) J=9000-IE-1
          if(iprint.gt.0)print *,' enter DO240:'
          DO 240 K=0,J
           IE=IE+1
           EMESH(IE)=EM-1./(F*F)
           IF(IPRINT.GT.0) WRITE(6,630)IE,EMESH(IE),N,M,F
  240      F=F+D
          IF(IPRINT.GT.0)WRITE(6,*)
          E=EM-SE
        ENDIF
      ENDIF
C
C  FIND NEXT HIGHER THRESHOLD
C
      N=M
      EN=EM
      DO 250 J=N+1,NAST
       IF(ENAT(J).EQ.EN) GO TO 250
        EM=ENAT(J)
        MM=J
        GOTO 260
  250  CONTINUE
      GOTO 280
  260 F=1./SQRT(EM-E)
      D=1./SQRT(EM-EMAX)-F
      J=D*QNINT+1
      IF(IPRINT.GT.0) WRITE(6,*)
      CUT=IE+J.GE.9000
      D=D/J
      IF(CUT) J=9000-IE-1
      DO 270 K=0,J
       IE=IE+1
       EMESH(IE)=EM-1./(F*F)
       IF(IPRINT.GT.0) WRITE(6,630)IE,EMESH(IE),N,MM,F
  270  F=F+D
      GO TO 400
C
  210 CONTINUE
C
C
C  UP TO EMAX USING DEOPEN
  280 IF(IPRINT.GT.0) WRITE(6,690)
      GO TO 320
C
C
C  ALL CHANNELS OPEN
C  *****************
C
  290 IF(IPRINT.LE.0) GO TO 320
  300 WRITE(6,'(31X,8HALL OPEN)')
  320 F=EMAX-E
      J=F/DEOPEN+1
      CUT=IE+J.GE.9000
      F=F/J
      IF(CUT) J=9000-IE-1
      DO 330 K=0,J
       IE=IE+1
       EMESH(IE)=E
       IF(IPRINT.GT.0) WRITE(6,630)IE,E
  330  E=E+F
C
C  TASK COMPLETED
C  **************
C
  400 IF(CUT) WRITE(6,670) EMESH(IE)
      MXE=IE
      WRITE(6,'(/'' NUMBER OF ENERGIES  MXE ='',I6/)') IE
      RETURN
C
C
  600 FORMAT(//1X,70(1H*)//20X,'ENERGY MESH'/20X,11(1H*)
     *//' DQN    =',F11.6/' QNMAX  =',F11.6/
     *  ' EMIN   =',F11.6/' EMAX   =',F11.6/' DEOPEN =',F11.6/)
  605 FORMAT(/' EMIN INCREASED TO E(NC1) BECAUSE INCOMPATIBLE WITH DQN')
  610 FORMAT('  LOWEST LEVEL ABOVE EMIN  NC1 =',I4)
  620 FORMAT(' HIGHEST LEVEL BELOW EMAX  NC2 =',I4)
  625 FORMAT(/' VALUES OF  -  IE, E = EMESH(IE), N = LOWEST'
     + ,' LEVEL ABOVE E,'/15X,'M = LEVEL USED FOR EFFECTIVE',
     +' QUANTUM NUMBER'/15X,'AND FNU = EFFECTIVE QUANTUM NUMBER'/
     + /3X,'IE    EMESH',17X,'N    M    FNU')
  630 FORMAT(I5,F12.6,10X,2I5,F12.6)
  670 FORMAT(//10X,54(1H*)/10X,
C    +'NUMBER OF ENERGIES EXCEEDS MAXIMUM OF 9000 ALLOWED BY DIMENSIONS'
     *'NUMBER OF ENERGIES REDUCED TO (AMP)MSH AT E =',F9.5/10X,54(1H*)/)
  690 FORMAT(/36X,'USING DEOPEN')
C
      END
C
C***************************************************************
C
      SUBROUTINE POINTS
C
C  LAST MODIFIED 26.2.86 & '93MAY09
C
C  CALCULATES CHANNELS ENERGIES, NUMBER OF OPEN CHANNELS
C  AND TABULAR POINTS
C     IMPLICIT REAL*8(A-H,O-Z)
C
      PARAMETER(TINY=-1.E-6)
      LOGICAL QDT,WARN
C
      COMMON/CQDT/R2ST( 500),QDT,NQ
      COMMON/CMESH/EMAX,EMIN,DEOPEN,DQN,QNMAX,EMESH(9000),IMESH
      COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,LACC
      COMMON/CPOINT/RZERO,RONE,RTWO,H,KP0,KP1,KP2
      COMMON/CHAN/ECH( 500),EPS( 500),FKNU( 500),CC( 500),
     1 RINF( 500),ITARG( 500),LLCH( 500),NCHF,NCHOP,NCHOP1
      COMMON/CEN/ETOT,MXE,NWT,NZ
      COMMON/CNTRL/IPRINT,IRAD,IPERT,KP2X
      COMMON/CINPUT/ CF( 500, 500,  20),COEFF(3,  25),ENAT(0: 109),
     5 WMAT( 500,7000),VALUE(7000),RA,BSTO,AZAZ,
     1 LAMAX,LRANG2,LRGL2,MNP2,NAST,NCHAN,NELC,NPTY2,NSPN2,NZED,MORE2,
     2 ISAT( 109),LAT( 109),NCONAT( 109)
      COMMON/CWARN/WARN
C
C
C
      KP0=1
      QDT=.FALSE.
C
C     FOR IRAD.GT.0 HOLD OLD KP2 RTWO AND H
      IF(IRAD.GT.0) THEN
       KP2OLD=KP2
       RTWOLD=RTWO
       HOLD=H
      ENDIF
C
C  CHANNEL ENERGIES EPS AND NUMBER OF OPEN CHANNELS NCHOP
      NCHOP=0
      DO 100 I=1,NCHF
       E=ETOT-ECH(I)
       IF(E.LE.TINY) GO TO 90
       IF(E.LT.0.) E=0.
C      HES'93MAY5 FIX FOR QDT,NQ NO LONGER NEEDED IN OUTJJM - RUB'96APR2
       FKNU(I)=SQRT(E)
       NCHOP=NCHOP+1
       GO TO 100
   90  FKNU(I)=1./SQRT(-E)
  100  EPS(I)=E
      NCHOP1=NCHOP+1
C
C  SET QDT FOR CASE OF IMESH.EQ.2
C  ASSUMES  LL IN ASCENDING ORDER
C  NO QDT CHANNEL FOR FNU.LT.LL+0.5
      IF(IMESH.GE.2.AND.NCHOP.LT.NCHF)THEN
        DO 120 I=NCHOP+1,NCHOP+NCONAT(ITARG(NCHOP+1))
         IF(FKNU(I)+0.00002.LE.MAX(QNMAX,LLCH(I)+.5)) GO TO 121
         QDT=.TRUE.
  120    NQ=I
  121   IF(QDT)THEN
          NCHOP1=NQ+1
          IPERT=0
          IF(IPRINT.GT.0) WRITE(6,630)ETOT,NCHOP,NQ
        ENDIF
      ENDIF
C
C
C  CALCULATION OF RTWO
C  INCLUDES CALCULATION OF INNER POINTS OF INFLECTION RINF
      RTWO=RZERO
C
C  (1) OPEN CHANNELS
C  FOR OPEN CHANNELS RTWO IS DEFINED BY CONVERGENCE CRITERION
C  FOR THE JWBK METHOD.
      DO 139 I=1,NCHOP
      E=EPS(I)
      C=CC(I)
      L=LLCH(I)
      EC=E*C
      IF(EC.GT.AC)THEN
        RINF(I)=(SQRT(1.+EC)-1.)/E
      ELSE
        RINF(I)=.5*C
      ENDIF
      IF(L-1)131,132,133
  131 CONST=56.
      IF(AC.GE.1.E-3) CONST=12.
      CE=CONST*E
      IF(CE.LT..1)THEN
        R2=.5*CONST
      ELSE
        R2=(SQRT(1.+CE)-1.)/E
      ENDIF
      GOTO 138
  132 CONST=16.
      IF(AC.GE.1.E-3) CONST=3.9
      GOTO 135
  133 IF(AC.GE.1.E-3) THEN
        CONST=5.7/L+1.2
      ELSE
        CONST=9.8/L+1.4
      ENDIF
  135 R2 = RINF(I)*CONST
  138 IF(RTWO.LT.R2)RTWO=R2
  139 R2ST(I)=R2
C
C  (2) CLOSED CHANNELS (STRONGLY CLOSED FOR QDT.EQ..TRUE.)
C  FOR CLOSED CHANNELS RTWO IS EQUAL TO THE OUTER POINT OF
C  INFLECTION, EXCEPT FOR THE CASE OF FNU.LT.(LL+1)
      DO 150 I=NCHOP1,NCHF
      FNU=FKNU(I)
      IF(FNU.LT.REAL(LLCH(I)+1)) THEN
        RINF(I)=0.
        R2=RZERO
      ELSE
        A1=SQRT(FNU*FNU-CC(I))
        R2=FNU*(FNU+A1)
        RINF(I)=FNU*(FNU-A1)
      ENDIF
      IF(RTWO.LT.R2)RTWO=R2
  150 R2ST(I)=R2
      IF(RONE.GT.RTWO)RTWO=RONE
C
C  CASE OF IRAD.GT.0
      IF(IRAD.GT.0)THEN
C  CHECK VALUE OF RTWO
      IF(RTWO.GT.RTWOLD)THEN
        IF(WARN)WRITE(6,600)ETOT,RTWOLD
        WARN=.FALSE.
      ENDIF
C  RE-INSTATE OLD KP2,RTWO AND H
      KP2=KP2OLD
      RTWO=RTWOLD
      H=HOLD
      RETURN
      ENDIF
C
C  FIND INTERVAL AND TABULAR POINTS BETWEEN RZERO AND RTWO
      IF(RTWO.LE.RZERO) THEN
         KP2=1
         H=0.
         RETURN
      ENDIF
      WM=0.
      DO 170 I=1,NCHF
       C=CC(I)
       E=EPS(I)
       X=1./RZERO
       W=ABS(E+X*(2.-C*X))
       X=1./RTWO
       WM=MAX(ABS((C*X-2.)*X-E),W,WM)
       IF(C.LT.RZERO.OR.C.GT.RTWO)GOTO 170
       W=ABS(1./C+E)
       IF(W.GT.WM)WM=W
  170  CONTINUE
C
      H=ACNUM/SQRT(WM)
C..... (AMB'97JUN: FATEFUL) MODIFICATION 2.2.87
C         IF(ABS(RTWO-RZERO).LT.2.*H) THEN  RTWO=RZERO ETC
C         SOLUTION OF IPERT=2 CATASTROPHE -- WE'97AUG5:
          RTWO = MAX(RZERO+H*2.,RTWO)
C.....    END MODIFICATION
      KP2 = ((INT((RTWO-RZERO)/H)-1)/4)*4+5
C
C  CHECK DIMENSIONS FOR NUMBER OF OUTER-REGION POINTS
      IF(KP2.GT.2000)THEN
        L=KP2
        R2=RTWO
        KP2=((2000-1)/4)*4+1
        RTWO=(KP2-1)*H+RZERO
C       IF(IPRINT.GT.0.AND.
        IF(IPERT.NE.0) THEN
         WRITE(6,610) R2,L, RTWO,KP2, ETOT
         DO 200 I=1,NCHF
          IF(R2ST(I).GT.RTWO) WRITE(6,620) I,R2ST(I)
  200     CONTINUE
        IPERT=0
        ENDIF
      ENDIF
C
      H=(RTWO-RZERO)/(KP2-1)
C
C  WEIGHTS FOR BODE RULE INTEGRATION
      CALL BODE
C
C
  600 FORMAT('  **WARNING** FOR ETOT =',
     1 F10.5,' RTWO REDUCED TO MAXIMUM VALUE OF',F8.2
     + /11X,'NO MORE SIMILAR WARNINGS GIVEN'/)
  610 FORMAT(/10X,'USE OF PERTURBATION REQUIRES RTWO =',
     *  F9.2,', KP2 =',I5/10X,
     * 'WHICH IS LARGER THAN MAXIMUM OF 2000 ALLOWED BY DIMENSION:'/
     *  10X,'SETTING IPERT = 0, RTWO =',F8.2,', KP2 =',I4,
     * '  AT ETOT =',1P,E9.3E1)
  620 FORMAT(10X,'CHANNEL',I4,' REQUIRES    =',F8.2,
     * ' FOR USE OF PERTURBATION')
  630 FORMAT(/' ***** QDT USED  AT ETOT =',F10.6/
     * 7X,'NCHOP =',I4/10X,'NQ =',I4/)
C
      RETURN
      END
C
C***************************************************************
C
      SUBROUTINE COUL(INOUT)
C
C  MODIFIED 21.11.85
C  MODIFIED 22.5.86 FOR CASE OF IRAD.NE.0
C
C  CALCULATION OF COULOMB FUNCTIONS
C
C     IMPLICIT REAL*8(A-H,O-Z)
C
      LOGICAL QDT
C
      COMMON/CQDT/R2ST( 500),QDT,NQ
      COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,LACC
      COMMON/COULSC/FS( 500,2000),FSP( 500),FC( 500,2000),FCP( 500)
      COMMON/CPOINT/RZERO,RONE,RTWO,H,KP0,KP1,KP2
      COMMON/CHAN/ECH( 500),EPS( 500),FKNU( 500),CC( 500)
     1  ,RINF( 500),ITARG( 500),LLCH( 500),NCHF,NCHOP,NCHOP1
      COMMON/CNTRL/IPRINT,IRAD,IPERT,KP2X
      COMMON/CEN/ETOT,IE,NWT,NZ
C
      DIMENSION FST(2000), INOUT( 500)
C
C
C  CALCULATE OPEN-CHANNEL SOLUTIONS AT RTWO
C  AND PERFORM INWARDS INTERGATIONS TO RZERO
C  (FOR CASE OF RINF.GT.RZERO, THE REGULAR SOLUTION IS CALCULATED
C   AT RZERO USING SERIES AND INTEGRATED OUTWARDS).
      IF(NCHOP.GT.0)THEN
        HM=-H
        J=0
        DO 190 I=1,NCHOP
C...    CASE OF R2ST(I).LE.RTWO, PERTURBATION MAY BE USED
        IF(R2ST(I).LE.RTWO)THEN
        CALL INJWBK(EPS(I),LLCH(I),J)
        CALL JWBK(RTWO,J,FS(I,KP2),FSP(I),FC(I,KP2),FCP(I))
        IF(KP2.GT.1)THEN
          IF(RINF(I).LE.RZERO)THEN
            CALL NUMSC(EPS(I),CC(I),RTWO,HM,KP2,I)
                                    INOUT(I)=0
          ELSE
            IF(IPERT.GT.0)THEN
C MJS         IPERT=-IPERT --- MJS'95JAN19: DELETE STATEMENT.
C CM FE16     RESTORED AGAINST PARABOLIC PROPAGATION -- RUB'96FEB27:
              IPERT=-IPERT
              IF(IPRINT.GT.0)WRITE(6,720)
            ENDIF
            CALL COULS(LLCH(I),EPS(I),RZERO,S,SP)
                                    INOUT(I)=1
            FSP(I)=SP
            CALL NUMS(EPS(I),CC(I),RZERO,H,1,KP2,S,SP,FST)
            W2=FC(I,KP2)*SP-FCP(I)*S
            DO 170 K=1,KP2
  170       FS(I,K)=FST(K)
            C=FC(I,KP2)
            CP=FCP(I)
            CALL NUMS(EPS(I),CC(I),RTWO,H,KP2,1,C,CP,FST)
            W0=C*FSP(I)-CP*FS(I,1)
            IF(IPRINT.GT.1)WRITE(6,600)I,RINF(I),W0,W2
            FCP(I)=CP
            DO 180 K=1,KP2
  180       FC(I,K)=FST(K)
          ENDIF
        ENDIF
C...    CASE OF R2ST(I).GT.RTWO, PERTURBATION CANNOT BE USED
        ELSE
          IPERT=0
C         print *,' AMB test: IPERT=0'
          IF(IRAD.EQ.0) THEN
            CALL SC(EPS(I),LLCH(I),RZERO,AC,
     +          FS(I,1),FSP(I),FC(I,1),FCP(I))
          ELSE
            IF(RINF(I).LT.RZERO) THEN
              CALL SC(EPS(I),LLCH(I),RZERO,AC,FS(I,1),FSP(I),FC(I,1),
     +          FCP(I))
              CALL NUMSC(EPS(I),CC(I),RZERO,H,KP2,I)
                                    INOUT(I)=0
            ELSE
              WRITE(6,620)ETOT,I,R2ST(I),RTWO,I,RINF(I),RZERO,KP2
              STOP
            ENDIF
          ENDIF
        ENDIF
  190   J=J+15
      ENDIF
C
C  CALCULATE QDT SOLUTIONS AT RZERO
      IF(QDT)THEN
        DO 200 I=NCHOP+1,NQ
        CALL SC(EPS(I),LLCH(I),RZERO,AC,
     +          FS(I,1),FSP(I),FC(I,1),FCP(I))
        IF(IRAD.GT.0) CALL NUMSC(EPS(I),CC(I),RZERO,H,KP2,I)
                                    INOUT(I)=0
  200   CONTINUE
      ENDIF
C
C  CALCULATE CLOSED CHANNEL SOLUTIONS AT RTWO
C  AND PERFORM INWARDS INTEGRATIONS TO RZERO
      DO 210 I=NCHOP1,NCHF
C...  CASE OF R2ST(I).LE.RTWO, PERTURBATION MAY BE USED
      IF(R2ST(I).LE.RTWO)THEN
      CALL THETA(RTWO,I,FS(I,KP2),FSP(I),FC(I,KP2),FCP(I),ICONV)
      IF(ICONV.EQ.0)THEN
        CALL NUMT(EPS(I),CC(I),RTWO,H,KP2,KP0,I)
                                    INOUT(I)=0
      ELSE
        IPERT=0
        IF(IRAD.GT.0) THEN
          WRITE(6,630)ETOT,I,LLCH(I),FKNU(I)
          STOP
        ENDIF
      ENDIF
C...  CASE OF R2ST(I).GT.RTWO, PERTUBATION CANNOT BE USED
      ELSE
      IPERT=0
      CALL SC(EPS(I),LLCH(I),RZERO,AC,FSA,FSPA,FCA,FCPA)
      SINF=SIN(3.141592654*FKNU(I))
      COSF=COS(3.141592654*FKNU(I))
      FS(I,1)=FCA*SINF-FSA*COSF
      FSP(I)=FCPA*SINF-FSPA*COSF
      IF(IRAD.GT.0) THEN
        S=FS(I,1)
        SP=FSP(I)
        CALL NUMS(EPS(I),CC(I),RZERO,H,1,KP2,S,SP,FST)
                                    INOUT(I)=1
        DO 220 K=1,KP2
  220   FS(I,K)=FST(K)
      ENDIF
      ENDIF
  210 CONTINUE
C
C
      IF(IPRINT.GT.1)  WRITE(6,700) ETOT,
     *    RTWO,KP2,H, (J,FS(J,1),FSP(J),FC(J,1),FCP(J),J=1,NCHF)
C
      RETURN
C
  600 FORMAT(24X,'I =',I3,', RINF =',F7.2,', W0 =',F9.6,', W2 =',F9.6/)
  620 FORMAT(//10X,30(1H*)//10X,'FOR ETOT =',E13.6/
     + 10X,'(R2ST(',I2,')=',F8.2,').GT.(RTWO =',F9.2,')  AND'/
     + 10X,'(RINF(',I2,') =',F9.2,').GT.(RZERO=',F8.2,')'/
     + 10X,'KP2 =',I5,', (AMPERSAND)PTS = 2000'/
     + 10X,'CANNOT CALCULATE RADIATIVE DATA FOR THIS CASE'/
     + 10X,'TRY LARGER VALUE OF (AMPERSAND)PTS'/)
  630 FORMAT(//10X,30(1H*)//10X,'FOR ETOT =',E13.6/
     + 10X,'CHANNEL',I4,' HAS'/
     + 10X,'CHANNEL ANGULAR MOMENTUM QUANTUM NUMBER =',I3/
     + 10X,'CHANNEL EFFECTIVE QUANTUM NUMBER =',F9.2/
     + 10X,'CANNOT CALCULATE RADIATIVE DATA FOR THIS CASE WITH ',
     + '(AMPERSAND)TET=  50'/10X,
     * 'TRY SMALLER VALUE OF QNMAX OR LARGER VALUE OF (AMPERSAND)TET'/)
  700 FORMAT(//11X,'E =',F11.6/11X,14('-')/
     *       /'   RTWO =',F11.6,',  KP2 =',I5,',  H =',F11.6/
     *       /10X,'COULOMB FUNCTIONS S,SP,C AND CP'/1P,(I5,4E14.5))
  720 FORMAT(5X,5(1H*),' REGULAR COULOMB FUNCTION FROM SERIES ',5(1H*)/)
C
      END
C***************************************************************
C
      SUBROUTINE INJWBK(E,L,J)
C
C  COMPUTES ARRAY D WHICH IS HELD IN COMMON/CJWBK/ AND
C  USED FOR CALCULATION OF JWBK FUNCTIONS.
C
C     IMPLICIT REAL*8(A-H,O-Z)
C  D HAS DIMENSIONS OF 15*NCHF
      COMMON/CJWBK/D(7500)
      COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,LACC
C
      D(J+1)=E
      D(J+4)=0.
C
C     TIDIED RUB'94JUN16:
      IF(E.EQ.0) THEN
       IF(L.EQ.0) GO TO 30
C       CASE OF L.GT.0 AND E.EQ.0
        C=L*(L+1)
      ELSE
       FK=SQRT(E)
       D(J+2)=FK
       D(J+3)=1./FK
       IF(L.EQ.0) GO TO 30
C       CASE OF L.GT.0 AND E.GT.0
        C=L*(L+1)
        A=1.+E*C
        D(J+7)=A
        A=3.*A
        D(J+8)=A-1.
        D(J+9)=A+1.
        D(J+10)=FK*C
        D(J+11)=-4.*E
        D(J+12)=-9.+2.*A
      ENDIF
      SC=SQRT(C)
      D(J+4)=C
      D(J+5)=SC
      D(J+6)=(C+.125)/SC
      D(J+13)=6.*C
      D(J+14)=-C*C
C
C  TERM IN ARG GAMMA ETC
   30 D(J+15)=ARGC(E,L,AC)
C
      RETURN
      END
C
C***************************************************************
C
      SUBROUTINE JWBK(R,J,S,SP,C,CP)
C
C  COMPUTES FUNCTIONS S AND C AND THEIR DERIVATIVES SP AND
C  CP USING IJWBK METHOD.
C  USES DATA IN ARRAY D WHICH IS HELD IN COMMON/CJWBK/
C  AND SHOULD HAVE BEEN COMPUTED IN SUBROUTINE INJWBK.
C
C     IMPLICIT REAL*8(A-H,O-Z)
C  D HAS DIMENSION 15*NCHF
      COMMON/CJWBK/D(7500)
      COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,LACC
C
C
C
      E=D(J+1)
      C=D(J+4)
      X=1./R
C
      IF(C.EQ.0)GOTO 30
      IF(E.EQ.0)GOTO 70
C
C  CASE OF C.GT.0 AND E.GT.0
      W=E+X*(2.-C*X)
      WH=SQRT(W)
      Z=R*WH
      FK=D(J+2)
      RK=R*FK
      RMC=R-C
      ALP=Z+RK
      CK=D(J+10)
C  COMPUTE PHASE
      P=Z+D(J+15)
C  LOG TERM
      B=FK*ALP
      IF(B.GT.ACJWBK)GOTO 10
      B=-B
      P=((((.2*B+.25)*B+.33333333)*B+.5)*B+1.)*ALP+P
      GOTO 20
   10 P=P+D(J+3)*LOG(1.+B)
C  ARCTAN TERM
   20 P = ATAN2((Z-FK*RMC)*D(J+5),CK*Z+RMC) * D(J+6) + P
C  CAP. PHI TERM
      P=P+((5.*RMC/(Z*Z))-(Z*D(J+9)+RK*D(J+8)+CK)/(ALP*D(J+7)))/(24.*Z)
C  COMPUTE AMPLITUDE
      A1=.0625*(X/W)**3
      CC=A1*(((D(J+14)*X+D(J+13))*X+D(J+12))*X+D(J+11))
      BB = A1*(((6.*D(J+14)*X+5.*D(J+13))*X+4.*D(J+12)) * X+3.*D(J+11))
      GOTO 100
C
   30 IF(E.EQ.0)GOTO 60
C
C  CASE OF C.EQ.0 AND E.GT.0
      W=2.*X+E
      WH=SQRT(W)
      Z=R*WH
      FK=D(J+2)
      RK=R*FK
      ALP=Z+RK
C  COMPUTE PHASE
      P=Z+D(J+15)
      B=FK*ALP
      IF(B.GT.ACJWBK)GOTO 40
      B=-B
      P=P+ALP*((((.2*B+.25)*B+.33333333)*B+.5)*B+1.)
      GOTO 50
   40 P=P+D(J+3)*LOG(1.+B)
   50 P=P+1/(4.*ALP)+(5.*R/(Z*Z)-2.*(Z+ALP)/ALP)/(24.*Z)
C  COMPUTE AMPLITUDE
      A1=.0625*(X/W)**3
      CC=A1*(-4.*E-3.*X)
      BB=-12.*A1*(E+X)
      GOTO 100
C
C  CASE OF C.EQ.0 AND E.EQ.0
   60 W=2.*X
      WH=SQRT(W)
      Z=R*WH
      P=2.*Z*(1.+.046875*X)+D(J+15)
      WMQ=1./SQRT(WH)
      ET=(1.+.0234375*X)*WMQ
      ZET=(1.-.046875*X)*WH
      ETP=.25*(1.-.0703125*X)*X*WMQ
      GOTO 110
C
C  CASE OF E.EQ.0 AND C.GT.0
   70 W=X*(2.-C*X)
      WH=SQRT(W)
      Z=R*WH
      RMC=R-C
C  COMPUTE PHASE
      P = ATAN2(D(J+5)*Z,RMC) * D(J+6) + 2.*Z+D(J+15)
     *   -(3.*R+C)/(24.*(RMC+R)*Z)
C  COMPUTE AMPLITUDE
      A1=.0625*(X/W)**3
      CC=((D(J+14)*X+D(J+13))*X-3.)*X*A1
      BB=((6.*D(J+14)*X+5*D(J+13))*X-12.)*X*A1
C
C  COMPLETE CALCULATION OF S,SP,C AND CP
 100  WMQ=1./SQRT(WH)
      ET=(1.-CC)*WMQ
      ETP=(.5*(X*X/W)*(1.-C*X)*(1.-13.*CC)+X*BB)*WMQ
      ZET=(1.+2.*CC)*WH
  110 SI=SIN(P)
      CO=COS(P)
      S=ET*SI
      C=ET*CO
      SP=ETP*SI+C*ZET
      CP=ETP*CO-S*ZET
C
      RETURN
      END
C***************************************************************
C
      FUNCTION ARGC(E,L,AC)
C
C  CALCULATES ARG(GAMMA(L+1-I/K)) -1/K -(1/K)*LN(K) - L*PI/2
C
C     IMPLICIT REAL*8(A-H,O-Z)
C
      IF(E.GT.0)GOTO 10
      ARGC=-(L+.25)*3.141592654
      RETURN
C
   10 FK=SQRT(E)
      ET=1./FK
      IP=L+1
      P=IP
      PP=IP*IP
C
      IF(AC.LT.1.E-4)GOTO 100
      A1=10.*SQRT(ET)-ET*ET
      IF(A1.GT.PP)GOTO 20
      X=PP*E
      XP1=X+1.
      XH=P*FK
      A=-1.5707963327*(P+L-.5)
      GOTO 200
   20 L1=IP
      IP=1.+SQRT(A1)
      P=IP
      PP=IP*IP
      X=PP*E
      XP1=X+1.
      XH=P*FK
      A=-1.570796327*(P+L-.5)
      L2=IP-1
      DO 30 I=L1,L2
   30 A=A+ATAN(ET/I)
      GOTO 200
C
  100 A1=35.*ET**.25-ET*ET
      IF(A1.GT.PP)GOTO 120
      X=PP*E
      XP1=X+1.
      XH=P*FK
      A=-1.570796327*(P+L-.5)
      GOTO 140
  120 L1=IP
      IP=1.+SQRT(A1)
      P=IP
      PP=IP*IP
      X=PP*E
      XP1=X+1.
      XH=P*FK
      A=-1.570796327*(L+P-.5)
      L2=IP-1
      DO 130 I=L1,L2
  130  A=A+ATAN(ET/I)
  140 A = (((5.*X-10.)*X+1.)*E*2.-(3.*X-1.)*XP1*XP1*7.)*E*
     *    FK*.000396825548*XP1**(-5) + A
C
  200 IF(FK*X*X*.1667*PP.GT.AC) GO TO 210
      A=(X-2.)*FK*.25*PP+A
      GOTO 220
  210 A=A-.5*ET*LOG(XP1)
  220 A1=(P-.5)*XH
      IF(A1*X*X.GT.AC) GO TO 230
      A=(1.-X*.33333333)*A1+A
      GOTO 240
  230 A=(P-.5)*ATAN(XH)+A
  240 ARGC=FK/((X+1.)*12.)+A
C
      RETURN
      END
C
C***************************************************************
C
      FUNCTION FCHI(E,L,AC)
C
C  CALCULATES REAL PART OF PSI(L+1+I*GAM) - LN(GAM)
C  WHERE E = 1/(GAM**2).
C  THIS IS REQUIRED FOR CALCULATION OF SCRIPT G.
C
C     IMPLICIT REAL*8(A-H,O-Y)
      IMPLICIT COMPLEX(Z)
C
      FCHI=0.
      IF(E.EQ.0)RETURN
C
      AC1=(20.*AC)**.333
C
      IF(E.GT.AC1)GOTO 100
C
      C=0.
      IF(L.EQ.0)GOTO 20
      A1=1.
      A2=-E
      A3=E+E
      DO 10 I=1,L
      A2=A2+A3
      A1=A1+A2
   10 C=C+I/A1
   20 FCHI = ((((1.05*E+1.)*E+2.1)*E+21.)*.003968253968+C)*E
      RETURN
C
  100 AC1=1./SQRT(AC1)
      FL=L+1
      IF(FL.GT.AC1)GOTO 300
C
      N=AC1
      FL=N+1
      DO 210 I=L+1,N
  210  FCHI = I/(I*I*E+1.)+FCHI
      FCHI=-FCHI*E
C
  300 X1=FL*E
      X = X1*FL+1.
      ZE=CMPLX(FL,1./SQRT(E))
      ZE=-1./(ZE*ZE)
      FCHI = REAL((((1.05*ZE+1.)*ZE+2.1)*ZE+21.)*ZE)*.003968253968 +
     * (LOG(X)-X1/X)*.5 + FCHI
C
      RETURN
      END
C**********************************************************
C
      SUBROUTINE THETA(R,I,T,TP,TD,TDP,ICONV)
C
C     IMPLICIT REAL*8(A-H,O-Z)
C
      COMMON/CHAN/ECH( 500),EPS( 500),FKNU( 500),CC( 500)
     1  ,RINF( 500),ITARG( 500),LLCH( 500),NCHF,NCHOP,NCHOP1
      COMMON/CEN/ETOT,MXE,NWT,NZ
      COMMON/CACC/AX,ACNUM,ACJWBK,ACZP,LACC
      COMMON/CTHET/BB( 500,  50),BG( 500,  50),MSUM( 500)
      COMMON/CNTRL/IPRINT,IRAD,IPERT,KP2X
      COMMON/CPOINT/RZERO,RONE,RTWO,H,KP0,KP1,KP2
C
C     CHANGES INSERTED TO FORCE CONVERGENCE WHEN NU<LL
C     THESE INVOLVE RESETTING AC AND M
      AC=AX
      FNU=FKNU(I)
      LL=LLCH(I)
      M =INT(FNU)+LL+12
      IF(M.LT.  50)GOTO 10
      IF(IPRINT.GT.0)WRITE(6,700)I,M
      GOTO 100
C
   10 DLR = LOG(R)
      R1=1./R
      Y = .5*FNU*R1
C
      F1=1./FNU
      A=FNU*R1-F1
      F2=F1*F1
      B = R*F2+DLR
      C=R1+F2
      E = .5*FNU*FNU*FNU
C
      BG(I,1)=E
      BG(I,2)=F2
      BB(I,1)=FNU
      BB(I,2)=1.
C
      CNTAC=0
      FL = LL
   15 A1=FL-FNU
      A2=FL+FNU+1.
      BN=-2.*FNU-1.
      YN=1.
C
      BET=1.
      GAM=0.
      W1 = C
C
      S=1.
      U=0.
      CX=0.
      CY=0.
      MN=0
C
C
      DO 20 N=3,M
      A1=A1+1.
      A2=A2-1.
      BN=BN+2.
      CN=A1*A2
      DN=BN+F1*CN
      GAM=CN*GAM+DN*BET
      BET=CN*BET
      YN=YN*Y
      U=U+BET*YN
      CY=CY+GAM*YN
      AN=1./(N-2)
      GAM=GAM*AN
      BG(I,N)=GAM
      BET=BET*AN
      BB(I,N)=BET
      S=S+BET*YN
      CX=CX+GAM*YN
      W0 = W1
      W1 = (U*CX-S*CY)*R1 + S*S*C
      IF(ABS((W1-W0)/W1).LT.AC) THEN
        IF(MN.GT.0) GO TO 30
        MN = N+1
      ELSE
        MN = 0
      ENDIF
   20 CONTINUE
C
C HES: PROBLEMS WITH LARGE L IF NU.LT.L+1, START FIX
      IF(FNU.LT.FL+0.01) THEN
       AC=AC*10.
C      M=M-1
       CNTAC=CNTAC+1
       IF(CNTAC.LT.2) GO TO 15
         WRITE(6,"(' SR THETA: NO CONV. AT R,E,W1,W0,I,L,AC='/
     *    1P,4E12.4,2I4,E9.1/' SET RONE.GT.R AND TRY AGAIN')")
     *     R,ETOT,W1,W0,I,LL,AC
         STOP
      ENDIF
C HES END FIX
C  NOT CONVERGED
      IF(IPRINT.GT.0) WRITE(6,610) M-2
      GO TO 100
C
C  SUMMATIONS CONVERGED
C  30 P=EXP(-R*F1)*R**FNU
C
C+++++++++++
C  MODIFICATION TO COPE WITH VAX EXPONENT RANGE
C
C  STATEMENT DELETED -
C  30 P=EXP(-R*F1)*R**FNU
C  RELACED BY 30 P=EXP(-.5*R*F1)*R**(.5*FNU) - REPL.95JUL15WE:
   30 P = EXP((DLR*FNU-R*F1)*.5)
      CFACT=1./P
      BB(I,2)=BB(I,2)*CFACT
      DO 40 J=3,MN
      BB(I,J)=BB(I,J)*CFACT
   40 BG(I,J)=BG(I,J)*CFACT
C
C  END MODIFICATION
C+++++++++++
      T=P*S
      TP = (A*S-R1*U)*P
      TD = (B*S+CX)*E*P
      TDP = ((A*B+C)*S - (B*U+CY)*R1 + A*CX)*E*P
      MSUM(I) = MN
      ICONV=0
      RETURN
C
C  USE SUBROUTINE SC AND SET IPERT = 0
  100 CALL SC(EPS(I),LL,RZERO,AC,FSA,FSPA,FCA,FCPA)
      SINF=SIN(3.141592654*FNU)
      COSF=COS(3.141592654*FNU)
      T=FCA*SINF-FSA*COSF
      TP=FCPA*SINF-FSPA*COSF
CVKL               FSAP WRONG! 19/03/88
      IPERT=0
      ICONV=1
      RETURN
C
  700 FORMAT(//10X,30(1H*)/10X,'SUBROUTINE THETA'/
     1 10X,'FOR I=',I3,' REQUIRE M=',I3/10X,'WHICH IS LARGER',
     *  ' THAN MAXIMUM VALUE OF   50 ALLOWED BY DIMENSIONS'/
     * 10X,'USING SUBROUTINE SC WITH IPERT = 0'/10X,30(1H*))
  610 FORMAT(//10X,30(1H*)//10X,'SUBROUTINE THETA'//
     * 10X,'SUMMATIONS NOT CONVERGED WITH',I4,' TERMS'/
     2 10X,'USING SUBROUTINE SC WITH IPERT = 1'/10X,30(1H*))
      END
C
C***************************************************************
C
      SUBROUTINE SC(E,L,R,AC,S,SP,C,CP)
C
C  COULOMB FUNCTIONS S AND C FROM POWER SERIES EXPANSIONS.
C
C     IMPLICIT REAL*8(A-H,O-Z)
C
C  POWER SERIES FOR F AND G
      CALL COULFG(L,E,R,AC,F,FP,G,GP,K,IERR,ACT)
      IF(IERR.EQ.2)THEN
        WRITE(6,600)L,E
        STOP
      ENDIF
      IF(IERR.EQ.1)WRITE(6,610)L,E,K,ACT
C
C  CASE OF E.GE.0
      IF(E.GE.0)THEN
C  CALCULATE CAP. A
      A=1.
      IF(L.GT.0)THEN
      A1=1.
      A2=-E
      A3=E+E
      DO 10 I=1,L
      A2=A2+A3
      A1=A1+A2
   10 A=A*A1
      ENDIF
C  CALCULATE SCRIPT G AND COULOMB FUNCTION H AND DERIVATIVE
      SG=A*FCHI(E,L,AC)*.318309886
      H=-G-SG*F
      HP=-GP-SG*FP
C  CALCULATE CAP B
      IF(E.GE.0.01) A=A/(1.-EXP(-6.28318531/SQRT(E)))
C
C  CASE OF E.LT.0
      ELSE
      CALL ABG(E,L,AC,A,BG)
      H=-(G+BG*F)
      HP=-(GP+BG*FP)
      ENDIF
C
C  COMPLETE CALCULATION OF S AND C
      C1=1.25331414*SQRT(A)
      S=C1*F
      SP=C1*FP
      C1=C1/A
      C=H*C1
      CP=HP*C1
      RETURN
  600 FORMAT(///10X,'SERIES IN COULFG NOT CONVERGED'/
     + 10X,'SMALL L =',I3,', EPS =',1P,E12.4/)
  610 FORMAT(//5X,'*** FUNCTIONS FROM COULFG INACCURATE ***'/5X,
     * 'SMALL L =',I3,', EPS =',1P,E12.4,', K =',I5,', ACTACC =',E12.4/)
C
      END
C
C**********************************************************
C
      SUBROUTINE ABG(E,L,AC,A,BG)
C
C  COMPUTES FUNCTION G(X,L) TO ACCURACY AC.
C
C     IMPLICIT REAL*8 (A-H,O-Z)
C
      X=1./SQRT(-E)
C
C  CALCULATION OF A AND EAC=E*A*C
      IF(L.GT.0)GOTO 2
      A=1.
      EAC=0.
      GOTO 20
    2 IF(L.GT.1)GOTO 4
      A=1.+E
      EAC=E
      GOTO 20
    4 IF(X.LT.REAL(L+1)) GOTO 12
      C=0.
      A=1.
      A1=1
      A2=-E
      A3=2.*E
      DO 10 I=1,L
      A2=A2+A3
      A1=A1+A2
      A=A*A1
   10 C=C+I/A1
      EAC=E*A*C
      GOTO 20
C  CASE OF X.LT.(L+1)
   12 A=1.
      A1=0.
      DO 16  I=1,L
      A=((I*I)*E+1.)*A
      A2=I
      DO 14 J=1,L
      IF(J.EQ.I) GOTO 14
      A2=((J*J)*E+1.)*A2
   14 CONTINUE
   16 A1=A1+A2
      EAC=E*A1
C
C  COMPUTE A1=PI*BG/A-E*C=1/(2*X)+PSI(X)-LN(X)
   20 A1=0.
C  TEST CONVERGENCE OF ASYMPTOTIC EXPANSION
      XN=(754.*AC)**(-.125)
      IF(X.GT.XN)GOTO 40
C  USE RECURRENCE FORMULAE
      N=XN-X+1
      XN=X+N
      EE=-1./(XN*XN)
      A1=A1-.5*(1./X+1./XN)+LOG(XN/X)
      IF(N.LT.2)GOTO 40
      N=N-1
      DO 30 I=1,N
   30 A1=A1-1./(I+X)
C  USE ASYMPTOTIC EXPANSION
   40 A1=A1+(((1.05*EE+1.)*EE+2.1)*EE+21)*EE*.003968253968
C
C  COMPLETE CALCULATION
      BG=(A*A1+EAC)*.318309886
      RETURN
      END
C
C**********************************************************
C
      SUBROUTINE COULFG(LL,EPS,RHO,ACC,F,FP,G,GP,K,IERR,ACTACC)
C
C  CALCULATES COULOMB FUNCTIONS F AND G AND THEIR DERIVATIVES
C
C  ORIGINAL VERSION PUBLISHED IN COMP. PHYS. COMM.25, 87, 1982.
C  PRESENT VERSION MODIFIED TO AVOID UNDERFLOW AND OVERFLOW
C  CONDITIONS IN THE SUMMATIONS OVER N OF
C         U(N)=A(N)*RHO**(N+L+1)
C     AND V(N)=D(N)*RHO**(N+L+1)
C  U(N) AND V(N) ARE CALCULATED RECURSIVELY.
C
C  A FURTHER MODIFICATION IS THAT THESE SUMMATIONS ARE
C  DONE IN DOUBLE PRECISION EVEN WHEN SINGLE PRECISION
C  IS USED ELSEWHERE.
C
C
C  INPUT -
C        LL=ANGULAR MOMENTUM QUANTUM NUMBER
C        EPS=Z-SCALED ENERGY IN RYDBERGS
C        RHO=Z-SCALED RADIAL VARIABLE IN ATOMIC UNITS
C        ACC=ACCURACY REQUIRED
C
C  OUTPUT -
C        F=REGULAR FUNCTION
C        FP=DERIVATIVE OF F
C        G=IRREGULAR FUNCTION
C        GP=DERIVATIVE OF G
C        K=NUMBER OF TERMS NEEDED IN EXPANSION
C        IERR=ERROR CODE
C        ACTACC=ACCURACY ACTUALLY ACHIEVED
C
C  CONVERGENCE CRITERION -
C        VALUE OF WRONSKIAN CONVERGED TO ACCURACY OF 0.5*ACC
C
C  ERROR CODES -
C        IERR=0, CONVERGED WITH ACTACC.LT.ACC
C        IERR=1, CONVERGED WITH ACTACC.GT.ACC
C        IERR=2, NOT CONVERGED WITH 101 TERMS IN MAIN SUMMATION
C
C  INITIALIZATION
C
C
C     IMPLICIT REAL*8 (A-H,O-Z)
C
      DOUBLE PRECISION DUBU0,DUBU1,DUBU2,DUBV0,DUBV1,DUBV2,
     +  DUBQ1,DUBF,DUBFP,DUBS,DUBSP
C
C  NOTE - IMPLICIT REAL*(A-H,O-Z) IS OPTIONAL BUT IT IS
C         RECOMMENDED THAT THE EXPLICIT DOUBLE PRECISION
C         STATEMENT SHOULD ALWAYS BE RETAINED.
C
      DATA R2PI,PS0/.159154943,-.154431330/
      IERR=0
      LP1=LL+1
      L2=2*LL
      L2P1=L2+1
      FL=LL
      FLP1=LP1
      FL2P1=L2P1
      E2=.5*EPS
      R2=2.*RHO
      ACC2=2.*ACC
C
C     INITIALIZE FA=FACTORIAL(2*LL+1)
C     AND PS=PSI(2*LL+2)+PSI(1)
C
      FA=1.
      PS=PS0
C
C
C  CALCULATE ALPHA(N) AND BETA(N) AND INITIALIZE S AND SP
C  CONTINUE CALCULATION OF FA AND PS
C
C     S AND SP FOR N=0
      X3=-L2
      X2=L2P1
      X1=-2.*R2**(-LP1)
      SP=X3*X1
      X1=R2*X1
      S=X1
C
C     INITIALIZE FOR COEFFICIENTS IN RECURSION FORMULAE
      P1=FL*E2
      P2=P1
      Q1=-E2
C
C     INITIALIZE ALPHA AND BETA
      ALP1=1.
      ALP2=1.+P2
      BET1=0.
      BET2=Q1
C
      IF(LL.EQ.0)GOTO 20
C
C     S AND SP FOR N=1
      X3=X3+2.
      X2=X2-1.
      X1=X1/X2
      SP=SP+X3*X1
      X1=R2*X1
      S=S+X1
C
C     LOOP FOR N=2 TO 2*LL
      DO 10 N=2,L2
C
C     CONTINUE CALCULATION OF FA AND PSI
      FN=N
      FA=FN*FA
      PS=PS+1./FN
C
C     CONTINUE CALCULATION OF S AND SP
      X3=X3+2.
      X2=X2-1.
      X1=X1/(X2*FN)
      SP=SP+X3*X1*ALP2
      X1=R2*X1
      S=S+X1*ALP2
C
C     COMPUTE COEFFICIENTS IN RECURSION FORMULAE
      P1=P1-E2
      P2=P2+P1
      Q1=Q1-E2
C     NOW HAVE P2=-N*(N-2*LL-1)*EPS/4
C     AND Q1=-N*EPS/2
C
C     NEW ALPHA AND BETA
      ALP0=ALP1
      ALP1=ALP2
      ALP2=ALP1+P2*ALP0
      BET0=BET1
      BET1=BET2
   10 BET2=BET1+P2*BET0+Q1*ALP0
C
C     NORMALIZE S AND SP, COMPLETE CALCULATION OF FA AND PS
      S=S*FA
      SP=SP*FA
      FA=FL2P1*FA
      PS=PS+1./FL2P1
C
C     COMPLETE CALCULATION OF ALPHA AND BETA
      P1=P1-E2
      P2=P2+P1
      Q1=Q1-E2
      ALP0=ALP1
      ALP1=ALP2
      BET0=BET1
      BET1=BET2
      BET2=BET1+P2*BET0+Q1*ALP0
C
C     NOW HAVE ALP1=ALPHA(2*LL+1)
C     AND BET1=BETA(2*LL+1), BET2=BETA(2*LL+2)
C
C     VALUE OF A=A(EPS,LL)
   20 A=ALP1
      A4=4.*A
      CL=2.*A*LOG(ABS(R2))
      CLP=2.*A/RHO
C
C  CALCULATE F AND FP AND CONTINUE CALCULATION OF S AND SP
C
C     CALCULATE A0,A1,D0,D1
      A0 = 2.**LP1/FA
      A1=-A0/FLP1
      PS=2.*PS*A
      D0=(BET1-PS)*A0
      D1=(BET2-PS-(2.+1./FLP1)*A)*A1
C
C     INITIALIZE F,FP, CONTINUE CALCULATION OF S,SP
C          -VALUES FOR N=0
C           U0 AND V0
      FNPLP1=FLP1
      C1=RHO**LL
      U0=A0*C1
      V0=D0*C1
      FP=FNPLP1*U0
      SP=SP+FNPLP1*V0
      U0=U0*RHO
      V0=V0*RHO
      F=U0
      S=S+V0
      W1=F*(CLP*F+SP)-FP*S
      NNN=0
C
C          - VALUES FOR N=1
C            U1 AND V1
      FNPLP1=FNPLP1+1.
      C1=C1*RHO
      U1=A1*C1
      V1=D1*C1
      FP=FP+FNPLP1*U1
      SP=SP+FNPLP1*V1
      U1=U1*RHO
      V1=V1*RHO
      F=F+U1
      S=S+V1
      W2=F*(CLP*F+SP)-FP*S
      DW2=ABS(W2-W1)
C
C     INITIALIZE FOR COEFFICIENTS IN RECURSION FORMULAE
      P1=-2.*FLP1
      P2=P1
      Q1=A4+2.*A*FL2P1
      REPS=RHO*EPS
C
C  CONVERT TO DOUBLE
                          DUBU0=U0
                          DUBU1=U1
                          DUBV0=V0
                          DUBV1=V1
                          DUBQ1=Q1
                          DUBS=S
                          DUBSP=SP
                          DUBF=F
                          DUBFP=FP
C     LOOP FOR N=2 TO 100
      DO 40 N=2,100
C
C     COMPUTE COEFFICIENTS IN RECURSION FORMULAE
      P1=P1-2.
      P2=P2+P1
      DUBQ1=DUBQ1+DBLE(A4)
C     NOW HAVE P2=-N*(N+2*LL+1)
C     AND DUBQ1=2*A*(2*N+2*LL+1)
C
C      COMPUTE DUBU2=U(N) AND DUBV2=V(N)
      DUBU2=(2.*DUBU1+DBLE(REPS)*DUBU0)/DBLE(P2)
      DUBV2=(2.*DUBV1+DBLE(REPS)*DUBV0+DUBQ1*DUBU2)/DBLE(P2)
C
C     INCREMENT DUBFP AND DUBSP
      FNPLP1=FNPLP1+1.
      DUBFP=DUBFP+DBLE(FNPLP1)*DUBU2
      DUBSP=DUBSP+DBLE(FNPLP1)*DUBV2
C
C     INCREMENT DUBF AND DUBS
      DUBU2=DUBU2*DBLE(RHO)
      DUBV2=DUBV2*DBLE(RHO)
      DUBF=DUBF+DUBU2
      DUBS=DUBS+DUBV2
C
C     CALCULATE WRONSKIAN
      W1=W2
      DW1=DW2
      W2=DUBF*(CLP*DUBF+DUBSP)-DUBFP*DUBS
      DW2=ABS(W2-W1)
C
C     CONVERGENCE TEST
      K=N+1
      IF(DW1.GT.ACC2)GOTO 30
      IF(DW2.LE.ACC2) GO TO 50
C
C  NEW DUBU0,DUBU1,DUBV0,DUBV1
   30 DUBU0=DUBU1
      DUBU1=DUBU2
      DUBV0=DUBV1
   40 DUBV1=DUBV2
C
C
C  NOT CONVERGED
C
      IERR=2
      ACTACC=ABS(0.25*W2-1.)
      GOTO 60
C
C  CONVERGED
C
   50 ACTACC=ABS(0.25*W2-1.)
      IF(ACTACC.GT.ACC)IERR=1
C
C  COMPLETE CALCULATION OF G AND GP
C
   60                    S=DUBS
                         SP=DUBSP
                         F=DUBF
                         FP=DUBFP
      G=(S+CL*F)*R2PI
      GP=(SP+CL*FP+CLP*F)*R2PI
C
      RETURN
C
      END
C
C**********************************************************
C
      SUBROUTINE COULS(LL,EPS,RHO,S,SP)
C
C  INSERTED 21.11.85
C
C  CALCULATES COULOMB FUNCTION S AND ITS DERIVATIVE SP
C  FROM POWER-SERIES EXPANSION.
C
C  THE STATEMENT "IMPLICIT REAL*8(A-H,O-Z)" CAN BE DELETED FOR
C  SINGLE PRECISION BUT IT IS RECOMMENDED TO KEEP THE
C  EXPLICIT DOUBLE PRECISION STATEMENT.  THE SUMMATION
C  FOR THE COULOMB FUNCTIONS IS THEN CARRIED OUT IN
C  DOUBLE PRECISION.
C     IMPLICIT REAL*8 (A-H,O-Z)
      DOUBLE PRECISION C1,D0,D1,D2,DM,DEPS,DRHO,F,FLP1,FNPLP1,FP,
     +  P1,P2,REPS,U0,U1,U2,UM
C
      CHARACTER*2 Q
      COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,LACC
      COMMON/CNTRL/IPRINT,IRAD,IPERT,KP2X
      LOGICAL T(0:2),TP(0:2)
C
C  INITIALISATIONS
      DRHO=RHO
      DEPS=EPS
      FLP1=LL+1
      ACC10=.1*AC
      DO 10 I=0,1
      T(I)=.FALSE.
   10 TP(I)=.FALSE.
C
C  POWER-SERIES EXPANSION
C  **********************
C
C   VALUES FOR N=0
      FNPLP1=FLP1
      C1=DRHO**LL
      U0=C1
      D0=FNPLP1*U0
      DM=ABS(D0)
      FP=D0
      U0=U0*DRHO
      UM=ABS(U0)
      F=U0
C
C   VALUES FOR N=1
      FNPLP1=FNPLP1+1.
      C1=C1*DRHO
      U1=-C1/FLP1
      D1=FNPLP1*U1
      DM=MAX(ABS(D1),DM)
      FP=FP+D1
      U1=U1*DRHO
      UM=MAX(ABS(U1),UM)
      F=F+U1
C
C   INITIALIZE FOR COEFFICIENTS IN RECURSION FORMULAE
      P1=-2.*FLP1
      P2=P1
      REPS=DRHO*DEPS
C
C   LOOP FOR N=2 TO 100
      DO 40 N=2,100
C     COMPUTE COEFFICIENTS IN RECURSION FORMULAE
      P1=P1-2.
      P2=P2+P1
C     NOW HAVE P2=-N*(N+2*LL+1)
C     COMPUTE U2 AND INCREMENT FP
      FNPLP1=FNPLP1+1.
      U2=(2.*U1+REPS*U0)/P2
      D2=FNPLP1*U2
      DM=MAX(ABS(D2),DM)
      FP=FP+D2
C     MODIFY U2 AND INCREMENT F
      U2=U2*DRHO
      UM=MAX(ABS(U2),UM)
      F=F+U2
C     TEST CONVERGENCE -- 2 FOR 10 LINES, WE'92SEP19:
      T(2)=ABS(U2).LT.ABS(F)*ACC10
      TP(2)=ABS(D2).LT.ABS(FP)*ACC10
      IF(.NOT.T(2))GOTO 21
      IF(.NOT.T(1))GOTO 21
      IF(.NOT.T(0))GOTO 21
      IF(IPRINT.LE.1) GO TO 50
        A=UM/ABS(F)
        Q='UM'
        GO TO 45
   21 IF(.NOT.TP(2))GOTO 23
      IF(.NOT.TP(1))GOTO 23
      IF(.NOT.TP(0))GOTO 23
      IF(IPRINT.LE.1) GO TO 50
        A=DM/ABS(FP)
        Q='DM'
        GO TO 45
C     NEW U0,U1,T AND TP
   23 U0=U1
      U1=U2
      DO 30 I=0,1
      T(I)=T(I+1)
   30 TP(I)=TP(I+1)
   40 CONTINUE
C
C  SERIES NOT CONVERGED
      WRITE(6,600)LL,EPS,RHO
      STOP
C
C  NORMALISATION
C  *************
C
C   NORMALISE FOR FUNCTIONS .5*F AND .5*FP
C
   45 WRITE(6,610) LL,EPS,Q,A
   50 S=F
      SP=FP
      DO 60 K=1,LL
      C=1./(K*(2*K+1))
      S=S*C
   60 SP=SP*C
C
C   CALCULATE CAP B = A AND FUNCTIONS S AND SP
      IF(EPS.GT.0)THEN
        A=1.
        IF(LL.GT.0)THEN
          A1=1.
          A2=-EPS
          A3=EPS+EPS
          DO 70 I=1,LL
          A2=A2+A3
          A1=A1+A2
   70     A=A*A1
        ENDIF
        IF(EPS.GE.0.01) A=A/(1.-EXP(-6.28318531/SQRT(EPS)))
      C=2.50662828*SQRT(A)
      ELSE
      C=2.50662828
      ENDIF
      S=S*C
      SP=SP*C
C
      RETURN
C
  600 FORMAT(//10X,60('*')//10X,'SERIES IN COULS NOT CONVERGED'
     + /10X,' LL =',I3,',  EPS =',1P,E13.5,',  RHO =',
     +  E13.5//10X,60('*')/)
  610 FORMAT(/5X,'SUBROUTINE COULS, LL =',I3,',  EPS =',
     *  1P,E13.5,2X,A2,' =',E10.2)
      END
C
C***************************************************************
C
      SUBROUTINE NUMSC(E,C,R1,H,N,I)
C
C  NUMEROV INTEGRATION OF COULOMB FUNCTIONS S AND C.
C  INTERVAL H, NUMBER OF POINTS N.
C  FUNCTIONS IN ARRAYS FS,FC
C  INWARDS INTEGRATION, FUNCTIONS CALCULATED FOR
C  K=N,(N-1),...,1.
C
C  MODIFIED TO PERFORM BOTH INWARD AND OUTWARD INTEGRATION DEPENDING
C  ON THE SIGN OF H: (H NEGATIVE)=INWARD; (H POSITIVE)=OUTWARD  1.5.86
C
C     IMPLICIT REAL*8(A-H,O-Z)
      COMMON/COULSC/FS( 500,2000),FSP( 500),FC( 500,2000),FCP( 500)
C
      V(X)=EQ+X*(Q2-X*CQ)
C
C  DETERMINE THE DIRECTION OF INTEGRATION AND SET INCREMENT
      IF(H.GT.0.0) THEN
        INC=1
        K=1
      ELSE
        INC=-1
        K=N
      ENDIF
      F1=FS(I,K)
      F1P=FSP(I)
      G1=FC(I,K)
      G1P=FCP(I)
      IF(N.EQ.1)RETURN
C
C  FUNCTIONS AT K=(K+INC)
      Q=H*H
      EQ=E*Q
      Q2=2.*Q
      CQ=C*Q
C
      X1=1./R1
      CX=C*X1
      HX=H*X1
      A1=-2.*HX*HX*H
      U1P=A1*(1.-CX)
      A1=-A1*HX
      U1PP=A1*(2.-3.*CX)
      A1=-A1*HX*6.
      U1PPP=A1*(1.-2.*CX)
      U1=V(X1)
C
      R2=R1+H
      X2=1./R2
      U2=V(X2)
C
      A2=1.+U2*.033333333
      B1=1.+U1*(.025*U1-.46666667)-.13333333*U1P-.025*U1PP
     1 +.0027777778*(4.*U1*U1P-U1PPP)
      C1 = H*(1.+U1*(.0027777778*U1-.13333333)-.05*U1P-.0083333333*U1PP)
C
      F2=(B1*F1+C1*F1P)/A2
      G2=(B1*G1+C1*G1P)/A2
      K=K+INC
      FS(I,K)=F2
      FC(I,K)=G2
C
      U3=U2
      U2=U1
      F3=F2
      F2=F1
      G3=G2
      G2=G1
      IF(N.EQ.2)GOTO 20
C
C  CONTINUE INTEGRATION
      EQ=EQ*.083333333
      Q2=Q2*.083333333
      CQ=CQ*.083333333
      U2=U2*.083333333
      U3=U3*.083333333
      R3=R2
C
      DO 10 M=3,N
      U1=U2
      U2=U3
      F1=F2
      F2=F3
      G1=G2
      G2=G3
      R3=R3+H
      X3=1./R3
      U3=V(X3)
      D3=1./(1.+U3)
      D2=(2.-10.*U2)*D3
      D1=(1.+U1)*D3
      F3=D2*F2-D1*F1
      G3=D2*G2-D1*G1
      K=K+INC
      FS(I,K)=F3
   10 FC(I,K)=G3
C
      U2=12.*U2
      U3=12.*U3
      EQ=12.*EQ
      Q2=12.*Q2
      CQ=12.*CQ
C
C  CALCULATE FINAL DERIVATIVE
   20 IF(INC.EQ.1) RETURN
      H=-H
      CX=C*X3
      HX=H*X3
      A1=-2.*HX*HX*H
      U3P=A1*(1.-CX)
      A1=-A1*HX
      U3PP=A1*(2.-3.*CX)
      A1=-A1*HX*6.
      U3PPP=A1*(1.-2.*CX)
      A2=1.+U2*.033333333
      B3=1.+U3*(.025*U3-.46666667)-.13333333*U3P-.025*U3PP
     1 +.0027777778*(4.*U3*U3P-U3PPP)
      C3 = H*(1.+U3*(.0027777778*U3-.13333333)-.05*U3P-.0083333333*U3PP)
      FSP(I)=(A2*F2-B3*F3)/C3
      FCP(I)=(A2*G2-B3*G3)/C3
      H=-H
C
      RETURN
      END
C***************************************************************
C
      SUBROUTINE NUMT(E,C,R1,HP,N1,N2,I)
C
C  NUMEROV INTEGRATION OF COULOMB FUNCTIONS.
C  THE INTERVAL HP IS POSITIVE.
C  INTEGRATION FROM TABULAR POINT N1 TO TABULAR POINT N2.
C  FUNCTIONS THETA,THETP STORED IN FS,FSP
C  FUNCTIONS THETAD,THETADP STORED IN FC,FCP
C  STARTS WITH FUNCTIONS AND DERIVATIVES AT N1 STORED IN FS, FSP
C  CALCULATES FUNCTIONS AT ALL POINT TO N2 AND DERIVATIVE AT
C  THE POINT N2.
C
C     IMPLICIT REAL*8(A-H,O-Z)
      COMMON/COULSC/FS( 500,2000),FSP( 500),FC( 500,2000),FCP( 500)
      COMMON/CTHET/BB( 500,  50),BG( 500,  50),MSUM( 500)
C
      V(X)=EQ+X*(Q2-X*CQ)
C
      N21=N2-N1
C
C  RENORMALISE FOR CASE OF N2.EQ.N1
      IF(N21.NE.0)GOTO 5
      W=FS(I,N1)*FCP(I)-FSP(I)*FC(I,N1)
      W1=1./SQRT(W)
      FS(I,N1)=FS(I,N1)*W1
      FC(I,N1)=FC(I,N1)*W1
      FSP(I)=FSP(I)*W1
      FCP(I)=FCP(I)*W1
      BB(I,2)=BB(I,2)*W1
      NM=MSUM(I)
      DO 2 M=3,NM
      BB(I,M)=BB(I,M)*W1
    2 BG(I,M)=BG(I,M)*W1
      RETURN
C
C  INTEGRATIONS FOR N2.NE.N1
   5  IP=ABS(N21)
      IS=N21/IP
      H=HP*IS
      K=N1
      IP=IP-1
C
      F1=FS(I,N1)
      F1P=FSP(I)
      G1=FC(I,N1)
      G1P=FCP(I)
C
C  FUNCTIONS AT K=(N1+IS)
      Q=H*H
      EQ=E*Q
      Q2=2.*Q
      CQ=C*Q
C
      X1=1./R1
      CX=C*X1
      HX=H*X1
      A1=-2.*HX*HX*H
      U1P=A1*(1.-CX)
      A1=-A1*HX
      U1PP=A1*(2.-3.*CX)
      A1=-A1*HX*6.
      U1PPP=A1*(1.-2.*CX)
      U1=V(X1)
C
      R2=R1+H
      X2=1./R2
      U2=V(X2)
C
      A2=1.+U2*.033333333
      B1=1.+U1*(.025*U1-.46666667)-.13333333*U1P-.025*U1PP
     1 +.0027777778*(4.*U1*U1P-U1PPP)
      C1=H*(1.+U1*(.0027777778*U1-.13333333)-.05*U1P-.0083333333*U1PP)
      P2=.033333333*Q
C     D1=Q*(-.46666667+.025*U1+.025*U1+.0027777778*U1P) - NOTE '98MAY26?
      D1 = ((.111111111*U1P+U1)*.025-.46666667)*Q
      E1 = (.0027777778*U1-.13333333)*Q*H
C
      F2=(B1*F1+C1*F1P)/A2
      G2=(B1*G1+C1*G1P+D1*F1+E1*F1P-P2*F2)/A2
      K=K+IS
      FS(I,K)=F2
      FC(I,K)=G2
C
C
      U3=U2
      U2=U1
      F3=F2
      F2=F1
      G3=G2
      G2=G1
      IF(IP.EQ.0)GOTO 20
C
C  CONTINUE INTEGRATION
      EQ=EQ*.083333333
      Q2=Q2*.083333333
      CQ=CQ*.083333333
      U2=U2*.083333333
      U3=U3*.083333333
      Q=Q*.083333333333
      R3=R2
C
      DO 10 M=1,IP
      U1=U2
      U2=U3
      F1=F2
      F2=F3
      G1=G2
      G2=G3
      R3=R3+H
      X3=1./R3
      U3=V(X3)
      D3=1./(1.+U3)
      D2=(2.-10.*U2)*D3
      D1=(1.+U1)*D3
      F3=D2*F2-D1*F1
      G3=D2*G2-D1*G1-Q*D3*(F3+10.*F2+F1)
      K=K+IS
      FC(I,K)=G3
   10 FS(I,K)=F3
C
      U2=12.*U2
      U3=12.*U3
      EQ=12.*EQ
      Q2=12.*Q2
      CQ=12.*CQ
      Q=12.*Q
C
C  CALCULATE FINAL DERIVATIVE
   20 H=-H
      CX=C*X3
      HX=H*X3
      A1=-2.*HX*HX*H
      U3P=A1*(1.-CX)
      A1=-A1*HX
      U3PP=A1*(2.-3.*CX)
      A1=-A1*HX*6.
      U3PPP=A1*(1.-2.*CX)
      A2=1.+U2*.033333333
      B3=1.+U3*(.025*U3-.46666667)-.13333333*U3P-.025*U3PP
     1 +.0027777778*(4.*U3*U3P-U3PPP)
      C3=H*(1.+U3*(.0027777778*U3-.13333333)-.05*U3P-.0083333333*U3PP)
      P2=.033333333*Q
      D3 = ((.111111111*U3P+U3)*.025-.46666667)*Q
      E3 = (.0027777778*U3-.13333333)*Q*H
      F3P=(A2*F2-B3*F3)/C3
      G3P=(A2*G2-B3*G3+P2*F2-D3*F3-E3*F3P)/C3
C
C  RE-NORMALISE CLOSED-CHANNEL FUNCTIONS
C ==  MODIFICATION FOR VAX CALCULATION OF
C ==  W1=1./SQRT(F3*G3P-F3P*G3)
      AMAX=1./MAX(ABS(F3),ABS(G3),ABS(F3P),ABS(G3P))
      AF3=F3*AMAX
      AG3=G3*AMAX
      AF3P=F3P*AMAX
      AG3P=G3P*AMAX
      W1=AMAX/SQRT(AF3*AG3P-AF3P*AG3)
C == END MODIFICATION
      FSP(I)=F3P*W1
      FCP(I)=G3P*W1
      IP=IP+2
      DO 30 J=1,IP
      FS(I,J)=FS(I,J)*W1
   30 FC(I,J)=FC(I,J)*W1
C  RE-NORMALISE COEFFICIENTS
      BB(I,2)=BB(I,2)*W1
      NM=MSUM(I)
      DO 40 M=3,NM
      BB(I,M)=BB(I,M)*W1
   40 BG(I,M)=BG(I,M)*W1
C
      RETURN
      END
C***************************************************************
C
      SUBROUTINE NUMS(E,C,R1,HP,N1,N2,F,FP,FST)
C
C  INSERTED 21.11.85
C
C  NUMEROV INTEGRATION OF COULOMB FUNCTION F.
C  INTEGRATION FROM POINT N1 TO N2.
C  INTERVAL HP IS POSITIVE.
C  INPUT FUNCTIONS F,FP AT N1 WHERE FP IS R-DERIVATIVES
C  FUNCTIONS AT INTERMEDIATE POINTS STORED IN FST, FDST
C  OUTPUT FUNCTIONS F,FP,FD,FDP AT N2.
C
C     IMPLICIT REAL*8(A-H,O-Z)
C
      DIMENSION FST(2000)
C
      V(X)=EQ+X*(Q2-X*CQ)
C
      N21=N2-N1
C
      F1=F
      F1P=FP
      FST(N1)=F1
      IF(N21.EQ.0)RETURN
C
      IP=ABS(N21)
      IS=N21/IP
      H=HP*IS
      K=N1
      IP=IP-1
C
C
C  FUNCTIONS AT K=(N1+IS)
      Q=H*H
      EQ=E*Q
      Q2=2.*Q
      CQ=C*Q
C
      X1=1./R1
      CX=C*X1
      HX=H*X1
      A1=-2.*HX*HX*H
      U1P=A1*(1.-CX)
      A1=-A1*HX
      U1PP=A1*(2.-3.*CX)
      A1=-A1*HX*6.
      U1PPP=A1*(1.-2.*CX)
      U1=V(X1)
C
      R2=R1+H
      X2=1./R2
      U2=V(X2)
C
      A2=1.+U2*.033333333
      B1=1.+U1*(.025*U1-.46666667)-.13333333*U1P-.025*U1PP
     1 +.0027777778*(4.*U1*U1P-U1PPP)
      C1=H*(1.+U1*(.0027777778*U1-.13333333)-.05*U1P-.0083333333*U1PP)
      P2=.033333333*Q
C
      F2=(B1*F1+C1*F1P)/A2
      K=K+IS
      FST(K)=F2
C
C
      U3=U2
      U2=U1
      F3=F2
      F2=F1
      IF(IP.EQ.0)GOTO 20
C
C  CONTINUE INTEGRATION
      EQ=EQ*.083333333
      Q2=Q2*.083333333
      CQ=CQ*.083333333
      U2=U2*.083333333
      U3=U3*.083333333
      Q=Q*.083333333333
      R3=R2
C
      DO 10 M=1,IP
      U1=U2
      U2=U3
      F1=F2
      F2=F3
      R3=R3+H
      X3=1./R3
      U3=V(X3)
      D3=1./(1.+U3)
      D2=(2.-10.*U2)*D3
      D1=(1.+U1)*D3
      F3=D2*F2-D1*F1
      K=K+IS
   10 FST(K)=F3
C
      U2=12.*U2
      U3=12.*U3
      EQ=12.*EQ
      Q2=12.*Q2
      CQ=12.*CQ
      Q=12.*Q
C
C  CALCULATE FINAL DERIVATIVE
   20 H=-H
      CX=C*X3
      HX=H*X3
      A1=-2.*HX*HX*H
      U3P=A1*(1.-CX)
      A1=-A1*HX
      U3PP=A1*(2.-3.*CX)
      A1=-A1*HX*6.
      U3PPP=A1*(1.-2.*CX)
      A2=1.+U2*.033333333
      B3=1.+U3*(.025*U3-.46666667)-.13333333*U3P-.025*U3PP
     1 +.0027777778*(4.*U3*U3P-U3PPP)
      C3=H*(1.+U3*(.0027777778*U3-.13333333)-.05*U3P-.0083333333*U3PP)
      P2=.033333333*Q
      F3P=(A2*F2-B3*F3)/C3
C
      F=F3
      FP=F3P
C
      RETURN
      END
C**********************************************************
C
      SUBROUTINE REACT(IOPT1)
C
C  CALCULATION OF REACTANCE MATRIX.  LAST REVISED 06.4.93
C
C     IMPLICIT REAL*8(A-H,O-Y)
      IMPLICIT COMPLEX(Z)
C
      LOGICAL NEWBUT, QDT
C
      COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,LACC
      COMMON/CQDT/R2ST( 500),QDT,NQ
      COMMON/CMESH/EMAX,EMIN,DEOPEN,DQN,QNMAX,EMESH(9000),IMESH
      COMMON/COULSC/FS( 500,2000),FSP( 500),FC( 500,2000),FCP( 500)
C  ***  NOTE CHNAGE OF CC TO CCT IN /CHAN/ ***
      COMMON/CHAN/ECH( 500),EPS( 500),FKNU( 500),CCT( 500)
     1  ,RINF( 500),ITARG( 500),LLCH( 500),NCHF,NCHOP,NCHOP1
      COMMON/CALP/ASS( 500, 500),ASC( 500, 500),ACS( 500, 500)
     1 ,ACC( 500, 500)
      COMMON/CEN/ETOT,MXE,NWT,NZ
      COMMON/CINPUT/ CF( 500, 500,  20),COEFF(3,  25),ENAT(0: 109),
     5 WMAT( 500,7000),VALUE(7000),RA,BSTO,AZAZ,
     1 LAMAX,LRANG2,LRGL2,MNP2,NAST,NCHAN,NELC,NPTY2,NSPN2,NZED,MORE2,
     2 ISAT( 109),LAT( 109),NCONAT( 109)
      COMMON/COMEGA/IE,IRAD0,LOM
      COMMON/CNTRL/IPRINT,IRAD,IPERT,KP2X
      COMMON/CPOINT/RZERO,RONE,RTWO,H,KP0,KP1,KP2
      COMMON/CPOT/BW( 500, 500),LAMP( 500, 500)
      COMMON /CTOP/ TOPS( 500, 500),LRGLAM,LITLAM(5886),NTOP(5886,2),
     * NTCHAN(0: 109),INDM
      COMMON/CDEC/ARAD(5886),ARDEC( 109),IRDEC,IPAR( 109),KCH( 500,3)
      COMMON/CLOGB/NEWBUT
C
      DIMENSION KTST( 109, 109), INOUT( 500), OMEGA(0:5886)
      DIMENSION S( 500),SP( 500),C( 500),CP( 500),SPP( 500),CPP( 500),
     *          A( 500, 500),B( 500, 500),RK( 500, 500),RMAT( 500, 500),
     *          CS( 500, 500),CSP( 500, 500),CC( 500, 500)
     1 ,CCP( 500, 500),CSPP( 500, 500),CCPP( 500, 500)
      DIMENSION DS( 500, 500),DSP( 500, 500),DC( 500, 500),
     1          DCP( 500, 500),DFPP( 500, 500)
      DIMENSION ITST(5886),JTST(5886),OMST(5886),ODIP(0:5886),MDIP(5886)
      DIMENSION Y( 500, 500),P( 500, 500),Q( 500, 500),D( 500, 500)
      DIMENSION ZCHICC(  20,  20),ZCHIOC( 500,  20),ZVAL(  20),ODX(5886)
C
      EQUIVALENCE (Y,CS),(P,CSP),(Q,CC),(D,CCP)
C
C
C  CALCULATE R-MATRIX
C
C    INITIALISE ELEMENTS ABOVE DIAGONAL TO ZERO
      DO 1 J=2,NCHF
       DO 1 I=1,J-1
    1   RMAT(I,J)=0.
C    INITIALISE DIAGONAL ELEMENTS TO BUTTLE CORRECTION
      IF(NEWBUT)THEN
        RA2=RZERO*RZERO
        DO 21 I=1,NCHF
         L=LLCH(I)+1
         NBUT=COEFF(3,L)
   21    RMAT(I,I) = BUT0(NBUT,COEFF(2,L)+RA2*EPS(I)) * COEFF(1,L)
      ELSE
        DO 2 I=1,NCHF
         L=LLCH(I)+1
         E=EPS(I)
    2    RMAT(I,I) = (COEFF(3,L)*E+COEFF(2,L))*E+COEFF(1,L)
      ENDIF
C     COMPLETE FOR UPPER TRIANGLE
      DO 3 K=1,MNP2
       V=1./(VALUE(K)-ETOT)
       DO 3 J=1,NCHF
        VI=V*WMAT(J,K)
        DO 3 I=1,J
    3    RMAT(I,J)=RMAT(I,J)+VI*WMAT(I,K)
C    SYMMETRISE
      DO 4 I=2,NCHF
       DO 4 J=1,I-1
    4   RMAT(I,J)=RMAT(J,I)
C
C  COULOMB FUNCTIONS AT RZERO
C
C       S(I) IS REGULAR FUNCTION FOR OPEN CHANNELS, DECAYING
C       FUNCTION THETA FOR CLOSED CHANNELS.
C       C(I) IS IRREGULAR FUNCTION FOR OPEN CHANNELS,
C       THETAD = ENERGY DERIVATIVE OF DECAYING FUNCTIONS FOR
C       CLOSED CHANNELS.
      CALL COUL(INOUT)
C
C  RESET IPERT AND NCHOP FOR QDT CASE
      IF(QDT) THEN
        IPERTO=IPERT
        NCHOPO=NCHOP
        IPERT=0
        NCHOP=NQ
      ENDIF
      IF(NCHOP.EQ.0) RETURN
C
      AR=1./RZERO
      AR2=AR*AR
      AR3=AR2*AR
      DO 10 I=1,NCHF
       S(I)=FS(I,1)
       SP(I)=FSP(I)-BSTO*S(I)
       HAM=(CCT(I)*AR-2.)*AR-EPS(I)
       IF(I.GE.NCHOP1.AND.IPERT.LE.0) GO TO 10
       C(I)=FC(I,1)
       CP(I)=FCP(I)-BSTO*C(I)
CVKL                     S(I) WRONG --- 19/03/88.
       CPP(I)=HAM*C(I)
   10  SPP(I)=HAM*S(I)
C
C  CALCULATE REACTANCE MATRIX
C
C..CASE OF IPERT.GT.0
      IF(IPERT.GT.0)THEN
      CALL ALPHA
C  CALCULATE PERTURBED FUNCTIONS
C     OPEN-OPEN
      IF(NCHOP.EQ.0)GOTO 240
      DO 210 J=1,NCHOP
       DO 200 I=1,NCHOP
        CS(I,J)=S(I)*ACS(I,J)-C(I)*ASS(I,J)
        CSP(I,J)=SP(I)*ACS(I,J)-CP(I)*ASS(I,J)
        CC(I,J)=S(I)*ACC(I,J)-C(I)*ASC(I,J)
        CCP(I,J)=SP(I)*ACC(I,J)-CP(I)*ASC(I,J)
        CSPP(I,J)=SPP(I)*ACS(I,J)-CPP(I)*ASS(I,J)
        CCPP(I,J)=SPP(I)*ACC(I,J)-CPP(I)*ASC(I,J)
        DS(I,J)=CS(I,J)
        DC(I,J)=CC(I,J)
        DSP(I,J)=CSP(I,J)
  200   DCP(I,J)=CCP(I,J)
       CS(J,J)=CS(J,J)+S(J)
       CSP(J,J)=CSP(J,J)+SP(J)
       CC(J,J)=CC(J,J)+C(J)
       CCP(J,J)=CCP(J,J)+CP(J)
       CSPP(J,J)=CSPP(J,J)+SPP(J)
  210  CCPP(J,J)=CCPP(J,J)+CPP(J)
C     CLOSED-OPEN
      IF(NCHOP.EQ.NCHF)GOTO 270
      DO 220 J=1,NCHOP
       DO 220 I=NCHOP1,NCHF
        CS(I,J)=C(I)*ASS(I,J)-S(I)*ACS(I,J)
        CSP(I,J)=CP(I)*ASS(I,J)-SP(I)*ACS(I,J)
        CC(I,J)=C(I)*ASC(I,J)-S(I)*ACC(I,J)
        CCP(I,J)=CP(I)*ASC(I,J)-SP(I)*ACC(I,J)
        CSPP(I,J)=CPP(I)*ASS(I,J)-SPP(I)*ACS(I,J)
        CCPP(I,J)=CPP(I)*ASC(I,J)-SPP(I)*ACC(I,J)
        DS(I,J)=CS(I,J)
        DC(I,J)=CC(I,J)
        DSP(I,J)=CSP(I,J)
  220   DCP(I,J)=CCP(I,J)
C     OPEN-CLOSED
      DO 230 J=NCHOP1,NCHF
       DO 230 I=1,NCHOP
        CC(I,J)=S(I)*ACS(I,J)-C(I)*ASS(I,J)
        CCP(I,J)=SP(I)*ACS(I,J)-CP(I)*ASS(I,J)
        CCPP(I,J)=SPP(I)*ACS(I,J)-CPP(I)*ASS(I,J)
        DC(I,J)=CC(I,J)
  230   DCP(I,J)=CCP(I,J)
C     CLOSED-   CLOSED
  240 DO 260 J=NCHOP1,NCHF
       DO 250 I=NCHOP1,NCHF
        CC(I,J)=C(I)*ASS(I,J)-S(I)*ACS(I,J)
        CCP(I,J)=CP(I)*ASS(I,J)-SP(I)*ACS(I,J)
        CCPP(I,J)=CPP(I)*ASS(I,J)-SPP(I)*ACS(I,J)
        DC(I,J)=CC(I,J)
  250   DCP(I,J)=CCP(I,J)
       CC(J,J)=CC(J,J)+S(J)
  260  CCP(J,J)=CCP(J,J)+SP(J)
C  CALCULATE MATRICES A AND B
  270 DO 280 J=1,NCHF
       DO 280 I=1,NCHF
        A(I,J)=CC(I,J)
        DO 280 K=1,NCHF
  280    A(I,J)=A(I,J)-RMAT(I,K)*CCP(K,J)
      DO 290 J=1,NCHOP
       DO 290 I=1,NCHF
        B(I,J)=CS(I,J)
        DO 290 K=1,NCHF
  290    B(I,J)=B(I,J)-RMAT(I,K)*CSP(K,J)
C
      ELSE
C..CASE OF IPERT.LE.0
C  CALCULATE MATRICES A AND B
      DO 310 J=1,NCHOP
       DO 300 I=1,NCHF
        A(I,J)=-RMAT(I,J)*CP(J)
  300   B(I,J)=-RMAT(I,J)*SP(J)
       A(J,J)=A(J,J)+C(J)
  310  B(J,J)=B(J,J)+S(J)
      DO 330 J=NCHOP1,NCHF
       DO 320 I=1,NCHF
  320   A(I,J)=-RMAT(I,J)*SP(J)
  330  A(J,J)=A(J,J)+S(J)
      ENDIF
C
C  COMPLETE CALCULATION OF REACTANCE MATRIX
C
      CALL INV(A,NCHF)
      DO 345 J=1,NCHOP
       DO 345 I=1,NCHF
        RKK=0.
        DO 340 K=1,NCHF
  340    RKK=RKK-A(I,K)*B(K,J)
  345   RK(I,J)=RKK
C
C     FOR CASE OF IPERT.GT.0 -- !MJS/PJS-JP'95JUN21/RUB'96FEB27:
C     SYMMETRIZE OPEN-OPEN PART OF REACTANCE MATRIX A-->RK '98MAY-JUNE
      IF(IPERT.GT.0) THEN
       DO 333 I=1,NCHOP-1
        DO 333 J=I+1,NCHOP
         RK(I,J) = (RK(I,J)+RK(J,I))*.5
  333    RK(J,I) = RK(I,J)
      ENDIF
C
C JAT IF (IOPT1.GT.9.AND.IMESH.NE.2) CALL OUTJJM(ETOT,RK,RK,NCHOP,'K  ')
C     MOVED PAST LABEL 910 BECAUSE MISPLACED - WE'93APR6
C     SECTION WRITING RK ONTO 9 (FOR STGBF) RESTORED WE'93APR20
      IF (IOPT1.GT.9) GO TO 395
      IF(IRAD.EQ.0) GOTO 395
C
C      WRITE REACTANCE MATRIX AND FUNCTIONS TO UNIT 9
C   PUT FUNCTIONS IN CS,CSP,CSPP - ORIGINAL CONTENTS DESTROYED
C
      IF(IPERT.EQ.1) THEN
C
C  USE A AS WORKSPACE - ORIGINAL CONTENTS DESTROYED
      DO 385 J=1,NCHOP
       DO 390 I=NCHOP1,NCHF
  390   A(I,J)=S(I)*RK(I,J)
       DO 380 I=1,NCHOP
  380   A(I,J)=C(I)*RK(I,J)
  385  A(J,J)=A(J,J)+S(J)
C
C  INITIALISE DFPP TO V * F (DO 350) ETC
      DO 370 J=1,NCHOP
       DO 370 I=1,NCHF
        VF=0.
        DO 350 K=1,NCHF
         IF(LAMP(I,K).EQ.2) THEN
          VF=VF-BW(I,K)*AR2*A(K,J)
         ELSE IF(LAMP(I,K).EQ.3) THEN
          VF=VF-BW(I,K)*AR3*A(K,J)
         ENDIF
  350    CONTINUE
        DFPP(I,J)=VF
        DO 360 K=1,NCHF
         RKK=RK(K,J)
         CS(I,J)=CS(I,J)+CC(I,K)*RKK
         CSP(I,J)=CSP(I,J)+CCP(I,K)*RKK
         CSPP(I,J)=CSPP(I,J)+CCPP(I,K)*RKK
         DS(I,J)=DS(I,J)+DC(I,K)*RKK
  360    DSP(I,J)=DSP(I,J)+DCP(I,K)*RKK
        DSP(I,J)=DSP(I,J)+BSTO*DS(I,J)
        DFPP(I,J)=DFPP(I,J)+(AR*(AR*CCT(I)-2.)-EPS(I))*DS(I,J)
  370   CSP(I,J)=CSP(I,J)+BSTO*CS(I,J)
C
      ELSE
      DO 450 J=1,NCHOP
       DO 430 I=1,NCHOP
        CS(I,J)=C(I)*RK(I,J)
        CSP(I,J)=(CP(I)+BSTO*C(I))*RK(I,J)
  430   CSPP(I,J)=CPP(I)*RK(I,J)
       DO 440 I=NCHOP1,NCHF
        CS(I,J)=S(I)*RK(I,J)
        CSP(I,J)=(SP(I)+BSTO*S(I))*RK(I,J)
  440   CSPP(I,J)=SPP(I)*RK(I,J)
       CS(J,J)=CS(J,J)+S(J)
       CSP(J,J)=CSP(J,J)+SP(J)+BSTO*S(J)
  450  CSPP(J,J)=CSPP(J,J)+SPP(J)
      ENDIF
C
  480 WRITE(9) QDT
      WRITE(9) NCHOP
      IF(NCHOP.EQ.0) RETURN
      WRITE(9)((RK(I,J),J=1,NCHOP),I=1,NCHF)
      WRITE(9)((CS(I,J),J=1,NCHOP),I=1,NCHF)
      WRITE(9)((CSP(I,J),J=1,NCHOP),I=1,NCHF)
      WRITE(9)((CSPP(I,J),J=1,NCHOP),I=1,NCHF)
C     WRITE(9)((FS(I,KR),KR=1,KP2),I=1,NCHOP)
C     WRITE(9)((FC(I,KR),KR=1,KP2),I=1,NCHOP)
C     WRITE(9)((FS(I,KR),KR=1,KP2),I=NCHOP1,NCHF)
C.....MODIFIED OUTPUT - AND SHORTENED (N1,N2!) '93APR20
      DO 490 I=1,NCHOP
       N1=1
       N2=2
       IF(INOUT(I).EQ.0) GO TO 490
       N1=KP2
       N2=N1-1
  490  WRITE(9)INOUT(I),FS(I,1),FS(I,2),FC(I,N1),FC(I,N2)
      DO 495 I=NCHOP1,NCHF
       N1=1
       N2=2
       IF(INOUT(I).NE.0) GO TO 495
       N1=KP2
       N2=N1-1
  495  WRITE(9)INOUT(I),FS(I,N1),FS(I,N2)
C.....END MODIFICATION
      IF(IPERT.GT.0) THEN
        WRITE(9)((DS(I,J),J=1,NCHOP),I=1,NCHF)
        WRITE(9)((DSP(I,J),J=1,NCHOP),I=1,NCHF)
        WRITE(9)((DFPP(I,J),J=1,NCHOP),I=1,NCHF)
      ENDIF
      IF(IRAD.EQ.2) RETURN
      IF(ITARG(NCHOP).EQ.ITARG(1)) RETURN
C
C
C  WRITE REACTANCE MATRIX
C
  395 IF(IPRINT.GT.0)THEN
      WRITE(6,707)ETOT,(J,J=1,NCHOP)
      DO 400 I=1,NCHF
       IF(I.EQ.1.OR.I.EQ.NCHOP1) WRITE(6,*)
  400  WRITE(6,'(I3,1P,(T4,7E11.3))') I,(RK(I,J),J=1,NCHOP)
      ENDIF
C
C  CALCULATE MATRICES P AND Q,
C    TRANSMISSION MATRIX  IS -2*I*(P+I*Q), I=SQRT(-1)
C
C  A=(1+RK**2)**(-1)
      CALL MATSQ(RK,A,NCHOP)
      DO 800 I=1,NCHOP
  800 A(I,I)=A(I,I)+1.
      CALL INV(A,NCHOP)
C  P=RK*A
      CALL MATMUL(RK,A,P,NCHOP)
C  Q=RK*P
      CALL MATMUL(RK,P,Q,NCHOP)
C..USE VARIATIONAL METHOD FOR IPERT=-1
      IF(IPERT.LT.0.)THEN
      CALL ALPHA
C  Y=-(F/V/F), F HAS K-MATRIX NORMALISATION,
C    V IS PERTURBATION POTENTIAL
      DO 830 I=1,NCHOP
      DO 830 J=1,NCHOP
      Y(I,J)=ASS(I,J)
      DO 810 K=1,NCHOP
      Y(I,J)=Y(I,J)+RK(I,K)*ACS(K,J)+ASC(I,K)*RK(K,J)
      DO 810 L=1,NCHOP
  810 Y(I,J)=Y(I,J)+RK(I,K)*ACC(K,L)*RK(L,J)
      DO 820 K=NCHOP1,NCHF
      Y(I,J)=Y(I,J)+ASS(I,K)*RK(K,J)+RK(K,I)*ASS(K,J)
      DO 820 L=1,NCHOP
  820 Y(I,J)=Y(I,J)+RK(I,L)*ACS(L,K)*RK(K,J)
     +   +RK(K,I)*ASC(K,L)*RK(L,J)
      DO 830 K=NCHOP1,NCHF
      DO 830 L=NCHOP1,NCHF
  830 Y(I,J)=Y(I,J)+RK(K,I)*ASS(K,L)*RK(L,J)
      DO 840 I=2,NCHOP
      DO 840 J=1,I-1
  840 Y(I,J)=Y(J,I)
C  B=Y-RK*Y*RK
      DO 850 J=1,NCHOP
       DO 850 I=1,NCHOP
        B(I,J)=Y(I,J)
        DO 850 K=1,NCHOP
         DO 850 L=1,NCHOP
  850     B(I,J)=B(I,J)-RK(I,K)*Y(K,L)*RK(L,J)
C  D=Y*RK+RK*Y
      DO 860 J=1,NCHOP
       DO 860 I=1,NCHOP
        D(I,J)=0.
        DO 860 K=1,NCHOP
  860    D(I,J)=D(I,J)+Y(I,K)*RK(K,J)+RK(I,K)*Y(K,J)
C  P=P+A*B*A, Q=Q+A*D*A
      DO 870 J=1,NCHOP
       DO 870 I=1,NCHOP
        DO 870 K=1,NCHOP
         DO 870 L=1,NCHOP
          P(I,J)=P(I,J)+A(I,K)*B(K,L)*A(L,J)
  870     Q(I,J)=Q(I,J)+A(I,K)*D(K,L)*A(L,J)
      ENDIF
C...
C
C  QDT CASE  AND IMESH=2, PUNCH CHI MATRICES
C  COMPUTE P AND Q SUCH THAT CHI=P+I*Q
C MJS IF(QDT)THEN -- MOVED BY HES PAST LABEL 920 - WHO FORGOT TO
C  WE              COMPENSATE FOR FACTOR 2.**2 AT LABEL 500! QUB'92SEP22
      DO 890 J=1,NCHOP
       DO 880 I=1,NCHOP
        T=-Q(I,J)*2.
        Q(I,J)=P(I,J)*2.
        P(I,J)=T
  880   Y(I,J)=T
  890  Y(J,J)=2.+Y(J,J)
      IF (IOPT1.GE.14) THEN
C      KAB'94: JBIN FACILITY IN ADAPTED FORM
       CALL OUTJJM(ETOT,P,Q,'BIN')
      ELSE IF(IOPT1.GT.9.AND.IMESH.EQ.2) THEN
       CALL OUTJJM(ETOT,P,Q,'CHI')
      ENDIF
C     REVISING OPTION 10/11 AND CORRECTING K OPTION -- WE'93APR2-3,5-6:
      IF(IMESH.EQ.2) GO TO 910
      IF(IOPT1/2.EQ.5) GO TO 900
      IF(IOPT1/2.NE.6) GO TO 910
      CALL INV(Y,NCHOP)
      CALL MATMUL(Q,Y,RK,NCHOP)
      CALL OUTJJM(ETOT,RK,RK,'K  ')
      GO TO 910
  900 CALL OUTJJM(ETOT,P,Q,'-T ')
C  WE    FOR DEBUGGING STGJJ, NEEDED IN HANNELORE'S JJXITC OF '92JUL28
C        AND '93APR2-3: FOR JAT/UGO - LI-LIKE FE!
      GO TO 910
C     IF(IOPT1.GT.9.AND.IMESH.EQ.2) CALL OUTJJM(ETOT,P,Q,'CHI')'97DEC19:
  910 DO 920 I=1,NCHOP
  920  P(I,I)=1.+P(I,I)
      IF(QDT)THEN
C  REVERT TO ORIGINAL IPERT, NCHOP
        NCHOP=NCHOPO
        IPERT=IPERTO
C  DIAGONALISE CHICC
        NCC=NQ-NCHOP
        IF(NCC.GT.  20)THEN
          WRITE(6,601)NCC
          STOP
        ENDIF
        DO 930 N1=1,NCC
         DO 930 N2=1,NCC
C NOT USED (PHN'94):  Z=CMPLX(P(NCHOP+N1,NCHOP+N2),Q(NCHOP+N1,NCHOP+N2))
  930     ZCHICC(N1,N2)=CMPLX(P(NCHOP+N1,NCHOP+N2),Q(NCHOP+N1,NCHOP+N2))
        CALL CEIGEN(ZCHICC,ZVAL,NCC,AC)
C  CALCULATE ZCHIOC
        DO 940 N1=1,NCHOP
        DO 940 N2=1,NCC
        ZCHIOC(N1,N2)=0.
        DO 940 K=1,NCC
  940   ZCHIOC(N1,N2)=ZCHIOC(N1,N2)+
     +  CMPLX(P(N1,NCHOP+K),Q(N1,NCHOP+K))*ZCHICC(K,N2)
C  RADIATIVE DECAYS
      FDEC=1.
      IF(IRDEC.EQ.1)
     *  FDEC = EXP( ARDEC(ITARG(NCHOP+1)) * FKNU(NCHOP+1)**3 )
C  CALCULATE AVERAGE COLLISION STRENGTH AND STORE IN RK
      DO 970 J=1,NCHOP
       DO 970 I=1,J
        RK(I,J)=P(I,J)*P(I,J)+Q(I,J)*Q(I,J)
        DO 950 K=1,NCC
         VV=ABS(ZVAL(K))
         IF(1.-VV.LT.AC) THEN
          VV=FDEC-1.
          DO 945 M=1,NCHOP
  945     VV=CONJG(ZCHIOC(M,K))*ZCHIOC(M,K)+VV
         ELSE
          VV=FDEC-VV*VV
         ENDIF
  950   RK(I,J)=RK(I,J)+ABS(ZCHIOC(I,K)*ZCHIOC(J,K))**2/VV
        DO 960 K1=1,NCC-1
        DO 960 K2=K1+1,NCC
        IF(1.-ABS(ZVAL(K1)).LT.AC .AND.
     *     1.-ABS(ZVAL(K2)).LT.AC) THEN
          ZVV=ZVAL(K1)*CONJG(ZVAL(K1)-ZVAL(K2))+FDEC-1.
          DO 955 M=1,NCHOP
  955     ZVV=CONJG(ZCHIOC(M,K1))*ZCHIOC(M,K1)+ZVV
        ELSE
          ZVV=FDEC-ZVAL(K1)*CONJG(ZVAL(K2))
        ENDIF
  960   RK(I,J)=RK(I,J)+2.*REAL(ZCHIOC(I,K1)*ZCHIOC(J,K1)*
     +  CONJG(ZCHIOC(I,K2)*ZCHIOC(J,K2))/ZVV)
  970   RK(I,J)=.25*NWT*RK(I,J)
        DO 980 I=2,NCHOP
        DO 980 J=1,I-1
  980   RK(I,J)=RK(J,I)
      ELSE
C
C  CALCULATE COLLISION STRENGTHS FOR NON-QDT CASE AND STORE IN RK
C
        DO 500 J=1,NCHOP
        DO 500 I=1,NCHOP
C MJS   RK(I,J)=(P(I,J)**2+Q(I,J)**2)*NWT --- CONSEQUENTIAL MODIFICATION
  500   RK(I,J)=(P(I,J)*P(I,J)+Q(I,J)*Q(I,J))*NWT*.25
C  WE   N.B.: NWT DECLARED IN SCALE2 AS 4 TIMES THE HALF-INTEGER WEIGHT.
      ENDIF
C
         IF(ITARG(1).EQ.ITARG(NCHOP)) RETURN
C
C    WRITE COLLISION STRENGTH
      IF(IPRINT.GT.0)THEN
C-     DO 520 IW=6,8,2
       IW=6
        WRITE(IW,640)
        DO 519 I=1,NCHOP
         DO 519 J=I,NCHOP
  519     WRITE(IW,650)ETOT,NSPN2,LRGL2,NPTY2,ITARG(I),LLCH(I),
     *    ITARG(J),LLCH(J),RK(I,J)
C-520   CONTINUE
      ENDIF
C
      LMODE = -1
      L=1
      IF(NSPN2.EQ.0) L=2
      IF (ABS(IPERT).EQ.2) THEN
       IF(LRGL2.EQ.LRGLAM-L-L) LMODE=2
       IF(LRGL2.EQ.LRGLAM-L) LMODE=1
       IF(LRGL2.EQ.LRGLAM) LMODE=0
       IF(LMODE.EQ.1) READ(1,REC=IE*3-1) (ODIP(I),I=0,LOM)
      ENDIF
C
C  FOR IABS(IPERT).EQ.2.AND.LRGL2.EQ.LRGLAM, TAKE OUT
C  ALLOWED TRANSITIONS WITH LLCH.GT.LITLAM
      IF (LMODE.EQ.0 .AND. .NOT.QDT) THEN
       DO 1000 I=1,INDM
        IT1=NTOP(I,1)
        IT2=NTOP(I,2)
        IF (NCONAT(IT1)*NCONAT(IT2).EQ.0) GO TO 1000
        K2=NTCHAN(IT2)
C        IF(K2.GT.NCHOP)GOTO 1000
  990    K1=NTCHAN(IT1)
  995     IF(MAX(LLCH(K1),LLCH(K2)).GT.LITLAM(I)) RK(K1,K2)=0.
CC        IF(NSPN2.NE.0) GO TO 1000
          K1=K1-1
          IF(K1.GT.NTCHAN(IT1-1)) GO TO 995
         K2=K2-1
         IF(K2.GT.NTCHAN(IT2-1)) GO TO 990
 1000   CONTINUE
      ENDIF
C
C
C  STORE OMEGA RESULTS IN COMPACT FORM
C
C  5886=( 109*( 109-1))/2
C
C  IT AND JT ARE TARGET STATES
C  OMEGA(IT,JT) FOR JT.GT.IT STORED IN
C  OMST(IJST)=OMEGA(ITST(IJST),JTST(IJST))
C
C  NUMBER OF OPEN TARGET LEVELS
      NTAROP=ITARG(NCHOP)
C
      IJST=0
      I2=0
      DO 1030 IT=1,NTAROP
      IF(NCONAT(IT).EQ.0)GOTO 1030
      I1=I2+1
      I2=I2+NCONAT(IT)
      J2=I2
      DO 1020 JT=IT+1,NTAROP
      IF(NCONAT(JT).EQ.0)GOTO 1020
      J1=J2+1
      J2=J2+NCONAT(JT)
      IJST=IJST+1
      OMT=0.
      KTST(IT,JT)=IJST
      DO 1010 I=I1,I2
      DO 1010 J=J1,J2
 1010 OMT=OMT+RK(I,J)
       OMST(IJST)=OMT
       ITST(IJST)=IT
       JTST(IJST)=JT
       MDIP(IJST) = 0
       IF(LMODE.LT.0) GO TO 1020
       ODX(IJST) = OMT
       DO 1015 I=1,INDM
        IF(NTOP(I,1).NE.IT) GO TO 1015
        IF(NTOP(I,2).NE.JT) GO TO 1015
        MDIP(IJST) = 1
        GO TO 1020
 1015   CONTINUE
 1020 CONTINUE
 1030 CONTINUE
C
C  TOP-UP FOR IABS(IPERT).EQ.2
C     N.B: LRGL2.EQ.LRGLAM ONLY IN BP MODE -
      IF ((ABS(IPERT).EQ.2.AND..NOT.QDT) .AND. LRGL2.GE.LRGLAM-1) THEN
         DO 3000 I=1,INDM
          IT2=NTOP(I,2)
          IF(IT2.GT.NTAROP)GOTO 3000
          IT1=NTOP(I,1)
C         IF(NCONAT(IT1)*NCONAT(IT2).EQ.0) GO TO 3000
          K=KTST(IT1,IT2)
          MDIP(K) = 1
          VV = LITLAM(I)*LITLAM(I)
          DO 2500 I2=NTCHAN(IT2-1)+1,NTCHAN(IT2)
           DO 2400 I1=NTCHAN(IT1-1)+1,NTCHAN(IT1)
C           TOP-UP IN LARGE L
            IF (LRGL2.EQ.LRGLAM) OMST(K) = OMST(K)+RK(I1,I2)*TOPS(I2,I1)
C           TOP-UP IN SMALL L
            IF (TOPS(I1,I2)) 2200,2400,2100
 2100        J=I2
             GO TO 2300
 2200        J=I1
 2300        OMST(K) = (VV*EPS(J)+1)*TOPS(I1,I2)*RK(I1,I2)+OMST(K)
 2400       CONTINUE
 2500      CONTINUE
 3000     CONTINUE
      ENDIF
C
C  WRITE TARGET STATES AND COLLISION STRENGTHS -- CHANGE -2-->-1 '94OCT
      IF (IPRINT) 1070,1060,1050
 1050  WRITE(6,777)
 1060  WRITE(6,6010)ETOT,QDT,IPERT,(ITST(K),JTST(K),OMST(K),K=1,IJST)
C
C  ADD TO OMEGA  THIS OPTION DELETED, HES JAN 92 - RESTORED '93APR20:
 1070 IF (IOPT1.GT.9) RETURN
C
C  ADD TO OMEGA
      READ(1,REC=IE*3) (OMEGA(I),I=0,LOM)
      DO 1040 KK=1,IJST
       IT=ITST(KK)
       JT=JTST(KK)
       K=((JT-2)*(JT-1))/2+IT
       IF(IRAD0.GE.0) K=JT+NAST*(IT-1)-(IT*(IT+1))/2
       OMEGA(K)=OMEGA(K)+OMST(KK)
       IF(LMODE.LE.0) GO TO 1040
        IF(MDIP(KK).EQ.0) THEN
         ODIP(K) = OMEGA(K)
        ELSE
         IF(LMODE.EQ.2) ODIP(K) = -OMEGA(K)
         IF(LMODE.EQ.1) ODIP(K) = ODIP(K)-ODX(KK)
        ENDIF
 1040  CONTINUE
      WRITE(1,REC=IE*3) (OMEGA(I),I=0,LOM)
      ODIP(0) = -OMEGA(0)
      IF(LMODE.GT.0) WRITE(1,REC=IE*3-LMODE) (ODIP(I),I=0,LOM)
C
      RETURN
C
C
  601 FORMAT(//10X,10(1H*),' NUMBER OF DEGENERATE CLOSED CHANNELS',
     * ' NCC =',I3/21X,'LARGER THAN DIMENSION VALUE OF DEG =   20'/)
  640 FORMAT(/7X,'ENERGY',13X,'S  L PI',6X,'TI LI    TJ LJ',7X,'OMEGA'/)
  650 FORMAT(5X,E14.6,5X,3I3,5X,2I3,3X,2I3,5X,E14.6)
  707 FORMAT(//'   REACTANCE MATRIX FOR ETOT =',F11.6//(I9,6I11))
  777 FORMAT(//'   ETOT   QDT  IPERT  ',
     1 'INITIAL AND FINAL TARGET LEVELS, AND COLLISION STRENGTHS')
 6010 FORMAT(/F10.6,L3,I5,1P,(T21,3(2I3,E11.3)))
C
      END
C
C***************************************************************
C
      SUBROUTINE RAD
C
C     IMPLICIT REAL*8(A-H,O-Z)
C
      COMMON/CINPUT/ CF( 500, 500,  20),COEFF(3,  25),ENAT(0: 109),
     5 WMAT( 500,7000),VALUE(7000),RA,BSTO,AZAZ,
     1 LAMAX,LRANG2,LRGL2,MNP2,NAST,NCHAN,NELC,NPTY2,NSPN2,NZED,MORE2,
     2 ISAT( 109),LAT( 109),NCONAT( 109)
C
      COMMON/CDEC/ARAD(5886),ARDEC( 109),IRDEC,IPAR( 109),KCH( 500,3)
C
C
C  RADIATIVE PROBABILITIES
      K=0
      IC0=NCONAT(1)
      DO 40 I=2,NAST
        JC0=0
        DO 30 J=1,I-1
          K=K+1
          IF(ARAD(K).GE.0.) GO TO 30
          DO 20 IC=IC0+1,IC0+NCONAT(I)
          DO 20 JC=JC0+1,JC0+NCONAT(J)
           IF(ABS(CF(IC,JC,1)).LT.1.E-9) THEN
            ARAD(K)=0.
           ELSE
            NT=KCH(JC,1)
            IF(KCH(IC,1).NE.NT) GO TO 20
            IF(NSPN2.NE.0) NT=2*LRGL2
            L1=KCH(IC,2)
            L2=KCH(JC,2)
            I1=KCH(IC,3)
            I2=KCH(JC,3)
            ARAD(K) = (.00729729*(ENAT(I)-ENAT(J)))**3 * CF(IC,JC,1)**2
     +     * (NT+1) / ((I1+1) * (((L1+L2)/2+1)/2) * WSQ(I1,I2,L1,L2,NT))
C            (2J+1)      GI      MAX(LL1,LL2)  W2(LAT1,LAT2,LL1,LL2;1,J)
            GOTO 30
           ENDIF
   20     CONTINUE
   30   JC0=JC0+NCONAT(J)
   40 IC0=IC0+NCONAT(I)
C
C  CHECK COMPLETENESS OF ARAD
      DO 50 K=1,(NAST*(NAST-1))/2
        IF(ARAD(K).GE.0.) GO TO 50
        IRDEC=0
        RETURN
   50 CONTINUE
      IRDEC=1
C
C  CALCULATE ARDEC
      C=6.283185307/((NZED-NELC)**2)
      K=0
      DO 70 I=2,NAST
      A=0.
      DO 60 J=1,I-1
      K=K+1
   60 A=A+ARAD(K)
   70 ARDEC(I)=A*C
C
C  WRITES
      WRITE(6,610)
      K=0
      DO 80 I=2,NAST
      DO 80 J=1,I-1
      K=K+1
      IF(ARAD(K).GT.0.) WRITE(6,620)I,J,ARAD(K),ARAD(K)/2.419E-17
   80 CONTINUE
      WRITE(6,630) (I,ISAT(I),LAT(I),IPAR(I),ARDEC(I),I=1,NAST)
      WRITE(6,640)
      RETURN
C
C  FORMATS
  610 FORMAT(/1X,71(1H+)//' *** GAILITIS-AVERAGE OMEGA',
     +' CALCULATED ALLOWING FOR RADIATIVE DECAYS ***'//8X,
     *'RADIATIVE PROBABILITIES ( UNITS U=2(!!)*RY/HBAR,  NO Z-SCALING )'
     *//14X,'I',9X,'J',12X,'A(I,J)/U',12X,'A(I,J)*SEC'/)
  620 FORMAT(I15,I10,1P,2E21.4)
  630 FORMAT(//12X,'TARGET STATES  --  NUMBERS ARDEC TO UNITS U*Z**2/(2*
     *PI)'//12X,'INDEX     2*S+1',7X,'L',6X,'PARITY P',10X,'ARDEC(I)'/
     *22X,'OR  P   AND J*2'/(I15,3I10,1P,E22.4))
  640 FORMAT(/1X,71(1H+)/)
C
      END
C
C**********************************************************
C
      SUBROUTINE TOP(LLL)
C
C     FOR LS COUPLING SEE V M BURKE AND M J SEATON IN JPB19(1986)L527-33
C     IMPLICIT REAL*8(A-H,O-Z)
C
      CHARACTER*1 LIT
      PARAMETER(NTAR= 109*2)
      LOGICAL PRN
      DIMENSION KTCHAN(0:NTAR)
      DATA KTCHAN(0)/0/
C
      COMMON/CPOT/BW( 500, 500),LAMP( 500, 500)
      COMMON/CEN/ETOT,MXE,NWT,NZ
      COMMON/CINPUT/ CF( 500, 500,  20),COEFF(3,  25),ENAT(0: 109),
     5 WMAT( 500,7000),VALUE(7000),RA,BSTO,AZAZ,
     1 LAMAX,LRANG2,LRGL2,MNP2,NAST,NCHAN,NELC,NPTY2,NSPN2,NZED,MORE2,
     2 ISAT( 109),LAT( 109),NCONAT( 109)
      COMMON/CHAN/ECH( 500),EPS( 500),FKNU( 500),CC( 500)
     1  ,RINF( 500),ITARG( 500),LLCH( 500),NCHF,NCHOP,NCHOP1
      COMMON /CTOP/ TOPS( 500, 500),LRGLAM,LITLAM(5886),NTOP(5886,2),
     * NTCHAN(0: 109),INDM
      COMMON/CDEC/ARAD(5886),ARDEC( 109),IRDEC,IPAR( 109),KCH( 500,3)
      COMMON/CNTRL/IPRINT,IRAD,IPERT,KP2X
C
      LRGLAM=LLL
      JTOT = LRGL2
      MBP=1
      IF(NSPN2.NE.0) THEN
        JTOT = LRGL2*2
        MBP=0
      ENDIF
C
C  FIND ALLOWED TRANSITIONS
C  ************************
C
      DO 20 I1=1,NCHF
      DO 20 I2=I1,NCHF
      TOPS(I1,I2) = 0.
      IF(LAMP(I1,I2).NE.2)GOTO 20
      IT1=ITARG(I1)
      IT2=ITARG(I2)
      IF (LAT(IT2)+LAT(IT1).EQ.0) GO TO 20
      IF (ABS(LAT(IT2)-LAT(IT1)).GT.2) GO TO 20
C     FOR BP CASE.
C     IF(ABS(ENAT(IT2)-ENAT(IT1)).LT.1.E-4) GO TO 20 -- rub'98May15:
      IF(ABS(ENAT(IT2)-ENAT(IT1)).LT.1.E-7) GO TO 20
      DO 10 I=1,INDM
      IF(NTOP(I,1).EQ.IT1.AND.NTOP(I,2).EQ.IT2) GO TO 20
   10 CONTINUE
      INDM=INDM+1
      NTOP(INDM,1)=IT1
      NTOP(INDM,2)=IT2
      LITLAM(INDM) = ABS(MIN(LAT(IT1),LAT(IT2))-MBP) + LRGLAM
      IF(MBP.NE.0) LITLAM(INDM)=LITLAM(INDM)/2
   20 TOPS(I2,I1) = 0.
      IF(LRGL2.LT.LRGLAM-1) RETURN
C     INCIDENTALLY TAKING CARE OF BP CASE LRGL2.NE.LRGLAM
C
C  INITIALISATIONS FOR LRGL2.GE.(LRGLAM-1)
C  ***************************************
C
C     RUB'96MAY1->99MAY10:
      IF(MBP.NE.0) WRITE(6,"(/8X,'N.B. VMB-MJS STYLE TOPPING!')")
C  CALCULATE INVERSE CHANNEL LIST
C     EXPLOITING K QN KCH(N,1)=-1 IN THE CASE OF LS COUPLING:
      NTCHAN(0)=0
      DO 40 I=1,NAST
        M = NTCHAN(I-1)+NCONAT(I)
        KTCHAN(I*2)=M
        KTCHAN(2*I-1)=M
        DO 30 N=NTCHAN(I-1)+1,M
         IF(KCH(N,1).GT.JTOT) GO TO 40
         KTCHAN(2*I-1)=N
   30    CONTINUE
   40   NTCHAN(I)=M
C
C  CASE OF LRGL2.EQ.(LRGLAM-1)  OR  LRGL2.EQ.LRGLAM:  TOP-UP IN SMALL L
C  ************************************************
C
      WRITE(6,620) LRGLAM
      W=1.0
      DO 90 I=1,INDM
       PRN = IPRINT.GT.-3 .OR. I.LE.4 .OR. I.GE.INDM-3
       IF (LRGL2.LT.LRGLAM .AND. LITLAM(I).NE.LRGLAM) GO TO 80
       IT1=NTOP(I,1)
       IT2=NTOP(I,2)
       IF(NCONAT(IT1)*NCONAT(IT2).EQ.0)GOTO 80
       LAM=LITLAM(I)*2
       DO 70 I1=NTCHAN(IT1-1)+1,NTCHAN(IT1)
        LL1=LLCH(I1)*2
        DO 70 I2=NTCHAN(IT2-1)+1,NTCHAN(IT2)
         IF(KCH(I2,1).NE.KCH(I1,1)) GO TO 70
         LL2=LLCH(I2)*2
         IF (LL1.NE.LAM .OR. LL2.NE.LAM-2) GO TO 50
          F = ENAT(IT2)-ENAT(IT1)
          LIT='A'
          GO TO 60
   50    IF (LL2.NE.LAM .OR. LL1.NE.LAM-2) GO TO 70
          F = ENAT(IT1)-ENAT(IT2)
          LIT='B'
   60    IF(LRGL2.EQ.LRGLAM) W=WWW(I1,I2,KCH,0,0,JTOT)
C LS:    IF(LRGL2.EQ.LRGLAM) W=WSQ(LAT(IT1)*2,LAT(IT2)*2,LL1,LL2,JTOT)
         F = 4./((LAM*LAM)*F*W)
         TOPS(I1,I2) = F
         IF(.NOT.PRN) GO TO 90
         WRITE(6,630) I,NTOP(I,1),NTOP(I,2),LITLAM(I),  I1,I2,LIT,F,  W
         GO TO 90
   70    CONTINUE
   80  IF(PRN) WRITE(6,630) I,NTOP(I,1),NTOP(I,2),LITLAM(I)
   90  CONTINUE
C
      IF (LRGL2.NE.LRGLAM) RETURN
C
C  CASE OF LRGL2.EQ.LRGLAM:  TOP-UP IN LARGE L
C  ***********************
C
      WRITE(6,640)
      DO 200 I=1,INDM
       IT1=NTOP(I,1)
       IT2=NTOP(I,2)
       IF (NCONAT(IT1)*NCONAT(IT2).EQ.0) GO TO 200
       PRN = IPRINT.GT.-3 .OR. I.LE.4 .OR. I.GE.INDM-3
       DO 100 K=0,MBP
C
      M=KTCHAN((IT2-1)*2+K)+1
      M2=KTCHAN(2*IT1+K-1)
      DO 110 N=KTCHAN((IT1-1)*2+K)+1,M2
       IF(LLCH(N).NE.LLCH(M)+1) GO TO 110
        M1=N
        GOTO 120
  110  CONTINUE
      GO TO 190
  120 IF(LLCH(M2).GT.LITLAM(I)) M2=M2-1
      IF(M2.LE.M1) GO TO 190
C
      LAM=0
      MAXL=KCH(M,3)+KCH(M,2)-MBP
C LS: MAXL=LAT2+LL2
      DO 160 I1=M1,M2
       IF(MAXL+LAM.LE.JTOT) GO TO 160
        F=0.
        W=1./WWW(M1,M,KCH,0,LAM,JTOT)
        DO 150 L=JTOT+2,MAXL+LAM,2
  150    F=F+W*WWW(M1,M,KCH,L-JTOT,LAM,L)
        I2=I1+M-M1
        IF(PRN) WRITE(6,650)I,I1,I2,F, M,M1,M2,MAXL,LAM
        TOPS(I2,I1)=F
  160   LAM=LAM+4
C
  190 M=KTCHAN((IT1-1)*2+K)+1
      M2=KTCHAN(2*IT2+K-1)
      DO 130 N=KTCHAN((IT2-1)*2+K)+1,M2
       IF(LLCH(N).NE.LLCH(M)+1) GO TO 130
        M1=N
        GOTO 140
  130  CONTINUE
      GO TO 100
  140 IF(LLCH(M2).GT.LITLAM(I)) M2=M2-1
      IF(M2.LE.M1) GO TO 100
C
      LAM=0
      MAXL=KCH(M,3)+KCH(M,2)-MBP
      DO 180 I2=M1,M2
       IF(MAXL+LAM.LE.JTOT) GO TO 180
        F=0.
        W=1./WWW(M,M1,KCH,0,LAM,JTOT)
        DO 170 L=JTOT+2,MAXL+LAM,2
  170    F=F+W*WWW(M,M1,KCH,L-JTOT,LAM,L)
        I1=I2+M-M1
        IF(PRN) WRITE(6,651)I,I1,I2,F, M,M1,M2,MAXL,LAM
        TOPS(I2,I1)=F
  180   LAM=LAM+4
  100 CONTINUE
  200 CONTINUE
C
      RETURN
C
  620 FORMAT(//10X,'TOP-UP FOR (LARGE L).GT.',I2,', (SMALL L).GT.LITLAM'
     *//'  ALLOWED NON-DEGENERATE (DIFF.GE.1.E-4) TRANSITIONS ARE'
     +//'  INDEX  TARGET STATES  LITLAM',7X,'TOP-UP IN SMALL L FOR'
     * /37X,'CHANNELS')
  630 FORMAT(2(I5,I9), I11,I5,6X,'TOP',A1,' =',1P,2E12.4)
  640 FORMAT(/15X,'TOP-UP IN LARGE L FOR'/15X,'INDEX     CHANNELS')
  650 FORMAT(I18,I9,I5,9X,'FTOPA =',1P,E11.4,1X,5I4)
  651 FORMAT(I18,I9,I5,9X,'FTOPB =',1P,E11.4,1X,5I4)
C
      END
C
C**********************************************************
C
      FUNCTION BUT0(NBUT,U)
C
C     IMPLICIT REAL*8(A-H,O-Z)
      LOGICAL POLE
      COMMON/CBUT/FKN(0:  48),UKN(0:  48)
C
      BUT0=0.
C
C  CASE OF U.GT.0.04
      IF(U.GT..04)THEN
        FK=SQRT(U)
        POLE=.FALSE.
        DO 10 N=0,NBUT
        IF(ABS(FK-FKN(N)).GT..3)THEN
          BUT0=BUT0+1./(U-UKN(N))
        ELSE
          POLE=.TRUE.
          D1=FK-FKN(N)
        ENDIF
   10   CONTINUE
        IF(POLE)THEN
          D2=D1*D1
          D=.33333333*D1*(1.+.066666667*D2*(1.+.0952381*D2))
          BUT0=2.*BUT0+(D+1./(2.*FK-D1))/FK
        ELSE
          BUT0=2.*BUT0+TAN(FK)/FK
        ENDIF
C
C  SUM FOR U.LT..04
      ELSE
        DO 20 N=0,NBUT
   20   BUT0=BUT0+1./(U-UKN(N))
C
C  CASE OF U.LT..04 AND U.GT.-.04
        IF(U.GT.-.04)THEN
          BUT0=2.*BUT0+1.+.33333333*U*(1.+.4*U)
C
C  CASE OF U.LT.-.04
        ELSE
          FK=SQRT(-U)
          BUT0=2.*BUT0+TANH(FK)/FK
        ENDIF
C
      ENDIF
C
      RETURN
      END
C
C**********************************************************
C
      SUBROUTINE ALPHA
C
C  CALCULATES ALPHA INTEGRALS
C
C     IMPLICIT REAL*8(A-H,O-Y)
      IMPLICIT COMPLEX(Z)
C
      COMMON/CALP/ASS( 500, 500),ASC( 500, 500),ACS( 500, 500)
     1 ,ACC( 500, 500)
C
      COMMON/COULSC/FS( 500,2000),FSP( 500),FC( 500,2000),FCP( 500)
      COMMON/CPOINT/RZERO,RONE,RTWO,H,KP0,KP1,KP2
      COMMON/CPOT/BW( 500, 500),LAMP( 500, 500)
      COMMON/CBODE/WBODE(2000)
      COMMON/CHAN/ECH( 500),EPS( 500),FKNU( 500),CC( 500)
     1  ,RINF( 500),ITARG( 500),LLCH( 500),NCHF,NCHOP,NCHOP1
      COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,LACC
      COMMON/CNTRL/IPRINT,IRAD,IPERT,KP2X
C
C OUT AIMAG(Z) = (0.,-1.)*Z -- KB'97APR15/17WE: CORR/REPL AIMAG BY  IMAG
      AIMAG(Z) = (0., 1.)*Z
C              = IMAG( Z) -- IF AVAILABLE AS GENERIC FUNCTION '93JUL28WE
C
C  INITIALISE ALPHA TO ZERO
      DO 10 J=1,NCHF
       DO 10 I=1,J
        ASS(I,J)=0.
        ASC(I,J)=0.
        ACS(I,J)=0.
   10   ACC(I,J)=0.
C
C
C  CONTRIBUTION FROM RZERO TO RTWO
C  START LOOP ON RADIAL POINTS
      IF(KP2.EQ.1)GOTO 105
      R=RZERO
      DO 100 K=1,KP2
      W=WBODE(K)
      X=1./R
C  OPEN-OPEN AND OPEN-CLOSED PARTS
      IF(NCHOP.EQ.0)GOTO 80
      DO 50 I=1,NCHOP
      DO 50 J=I,NCHF
      LIJ=LAMP(I,J)
      IF(LIJ.EQ.1)GOTO 50
      WIJ=W*BW(I,J)*X**LIJ
      ASS(I,J)=ASS(I,J)+FS(I,K)*WIJ*FS(J,K)
      ASC(I,J)=ASC(I,J)+FS(I,K)*WIJ*FC(J,K)
      ACS(I,J)=ACS(I,J)+FC(I,K)*WIJ*FS(J,K)
      ACC(I,J)=ACC(I,J)+FC(I,K)*WIJ*FC(J,K)
   50 CONTINUE
C  CLOSED-CLOSED PART
      IF(NCHOP.EQ.NCHF)GOTO 100
   80 DO 90 J=NCHOP1,NCHF
      DO 90 I=NCHOP1,J
      LIJ=LAMP(I,J)
      IF(LIJ.EQ.1)GOTO 90
      WIJ=W*BW(I,J)*X**LIJ
      ASS(I,J)=ASS(I,J)+FS(I,K)*WIJ*FS(J,K)
      ACS(I,J)=ACS(I,J)+FC(I,K)*WIJ*FS(J,K)
      ASC(I,J)=ASC(I,J)+FS(I,K)*WIJ*FC(J,K)
   90 CONTINUE
  100 R=R+H
C
C
C  ASYMPTOTIC INTEGRALS FOR R=RTWO TO INFINITY
  105 IF(IPRINT.GT.1)WRITE(6,750)
C  - OPEN-OPEN PART
      IF(NCHOP.EQ.0)GOTO 180
      DO 150 J=1,NCHOP
      DO 150 I=1,J
      LIJ=LAMP(I,J)
      IF(LIJ.EQ.1) GO TO 150
      NLAG=2*LIJ+LACC
      NLEG=NLAG
      BIJ=BW(I,J)
C  P INTEGRALS
  110 ZP3=ZPLAG(I,J,NLAG,LIJ)*BIJ
C         FORMAT(70X,2I3,' P=',2E12.4)
C  Q INTEGRALS
      ALP = (FKNU(I)-FKNU(J))*RTWO
      IF(ALP.LT.2.) GO TO 141
  120 ZQ3=ZQLAG(I,J,NLAG,LIJ)*BIJ
      GO TO 145
  141 ZQ3=ZQLEG(I,J,NLEG,LIJ)*BIJ
C         FORMAT(70X,2I3,' Q=',2E12.4)
  145 ASS(I,J)=ASS(I,J)+.5*REAL(ZQ3-ZP3)
      ASC(I,J)=ASC(I,J)+.5*AIMAG(ZP3+ZQ3)
      ACS(I,J)=ACS(I,J)+.5*AIMAG(ZP3-ZQ3)
      ACC(I,J)=ACC(I,J)+.5*REAL(ZP3+ZQ3)
      IF(IPRINT.GT.1)WRITE(6,760)I,J,ASS(I,J),ACS(I,J),ASC(I,J),ACC(I,J)
  150 CONTINUE
C  - OPEN-CLOSED PART
      IF(NCHOP.EQ.NCHF)GOTO 300
      DO 170 J=NCHOP1,NCHF
      DO 170 I=1,NCHOP
      LIJ=LAMP(I,J)
      IF(LIJ.EQ.1) GO TO 170
      BIJ=BW(I,J)
      NLAG=2*LIJ+LACC
      CALL ZSLAG(I,J,NLAG,LIJ,ZS3,ZSD3)
      ZS3=ZS3*BIJ
      ZSD3=ZSD3*BIJ
C         FORMAT(70X,2I3,' S=',2E12.4)
C         FORMAT(70X,2I3,' SD=',2E12.4)
      ACC(I,J)=ACC(I,J)+REAL(ZSD3)
      ASC(I,J)=ASC(I,J)+AIMAG(ZSD3)
      ACS(I,J)=ACS(I,J)+REAL(ZS3)
      ASS(I,J)=ASS(I,J)+AIMAG(ZS3)
      IF(IPRINT.GT.1)WRITE(6,760)I,J,ASS(I,J),ACS(I,J)
  170 CONTINUE
C
C  CLOSED-CLOSED PART
  180 IF(NCHOP.EQ.NCHF)GOTO 300
      DO 200 J=NCHOP1,NCHF
      DO 200 I=NCHOP1,J
      LIJ=LAMP(I,J)
      IF(LIJ.EQ.1)GOTO 200
      NLAG=2*LIJ+LACC
      BIJ=BW(I,J)
      CALL TLAG(I,J,NLAG,LIJ,T1,T2,T3)
      ASS(I,J)=ASS(I,J)+T1*BIJ
      ACS(I,J)=ACS(I,J)+T2*BIJ
      ASC(I,J)=ASC(I,J)+T3*BIJ
      IF(IPRINT.GT.1)WRITE(6,760)I,J,ASS(I,J),ACS(I,J),ASC(I,J)
  200 CONTINUE
C
C
C  SYMMETRISE ALPHA
  300 IF(NCHF.EQ.1)RETURN
      DO 310 I=2,NCHF
      DO 310 J=1,I-1
      ASS(I,J)=ASS(J,I)
      ACC(I,J)=ACC(J,I)
      ASC(I,J)=ACS(J,I)
  310 ACS(I,J)=ASC(J,I)
C
C
      RETURN
  750 FORMAT(/' I,J AND ASS(I,J), ACS(I,J), ASC(I,J), ACC(I,J)'/)
  760 FORMAT(2I5,1P,4E14.5)
      END
C
C***************************************************************
C
      BLOCK DATA
C     QUB90'FEB1, WE: BECAUSE OF IBM COMPILERS NOT CALLED AS
CX    SUBROUTINE (BLOCK TO AVOID LINKAGE PROBLEMS WITH LIBRARIES)
C
C  DATA FOR QUADRATURES -
C  LAGUERRE AND LEGENDRE QUADRATURES WITH NUMBERS OF POINTS
C  N = 2, 4, 6, 8 AND 10
C
C     IMPLICIT REAL*8(A-H,O-Z)
      COMMON/CBLK/XLAG(30),WLAG(30),XLEG(15),WLEG(15)
C
      DATA XLAG/
     1 .58578644,3.4142136,
     2 .32254769,1.7457611,4.5366203,9.3950709,
     3 .22284660,1.1889321,2.9927363,5.7751436,9.8374674,
     4 15.982874,
     5 .17027963,.90370178,2.2510866,4.26670017,7.0459054,
     6 10.758516,15.7406786,22.8631317,
     7 .13779347,.72945455,1.8083429,3.4014337,5.5524961,
     8 8.3301527,11.8437858,16.279258,21.996586,29.920697/
      DATA WLAG/
     1 .85355339,.14644661,
     2 .60315410,.35741869,.38887909E-1,.53929471E-3,
     3 .45896467,.41700083,.11337338,.10399197E-1,
     4 .26101720E-3,.89854791E-6,
     5 .36918859,.41878678,.17579499,3.3343492E-2,2.7945362E-3,
     6 9.0765088E-5,8.4857467E-7,1.0480012E-9,
     7 .30844112,.40111993,.21806829,6.2087456E-2,9.5015170E-3,
     8 7.5300839E-4,2.8259233E-5,4.2493140E-7,1.8395648E-9,
     9 9.9118272E-13/
      DATA XLEG,WLEG/.577350269,
     1 .339981044,.861136312,
     2 .238619186,.661209386,.932469514,
     3 .183434642,.525532410,.796666477,.960289856,
     4 .148874339,.433395394,.679409568,.865063367,.973906529,
     5 1.,
     6 .652145159,.347854845,
     7 .467913935,.360761573,.171324492,
     8 .362683783,.313706646,.222381034,.101228536,
     9 .295524225,.269266719,.219086363,.149451349,.066671344/
C
CX    RETURN
      END
C
C***************************************************************
C
      SUBROUTINE BODE
C
C     IMPLICIT REAL*8(A-H,O-Z)
      COMMON/CPOINT/RZERO,RONE,RTWO,H,KP0,KP1,KP2
      COMMON/CBODE/WBODE(2000)
C
      W=.31111111*H
      WBODE(1)=W
      WBODE(KP2)=W
      W=1.4222222*H
      M=KP2-1
      DO 10 K=2,M,2
   10 WBODE(K)=W
      W=.53333333*H
      M=KP2-2
      DO 20 K=3,M,4
   20 WBODE(K)=W
      IF(KP2.EQ.5)RETURN
      W=.62222222*H
      M=KP2-4
      DO 30 K=5,M,4
   30 WBODE(K)=W
C
      RETURN
      END
C
C***************************************************************
C
      SUBROUTINE ZPHI(I,ZR,ZAI,ZPI)
C
C  COMPUTES AMPLITUDE ZAI AND PHASE ZPI OF COULOMB FUNCTION ZPHI
C  FOR COMPLEX RADIAL CO-ORDINATE ZR.
C  USES DATA IN ARRAY D WHICH IS HELD IN COMMON/CJWBK/
C  AND SHOULD HAVE BEEN COMPUTED IN SUBROUTINE INJWBK.
C  THE STRUCTURE OF ZPHI IS SIMILAR TO THAT OF SUBROUTINE JWBK.
C
C
C     IMPLICIT REAL*8(A-H,O-Y)
      IMPLICIT COMPLEX(Z)
      COMMON/CJWBK/D(7500)
      COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,LACC
C
C
      J=(I-1)*15
      E=D(J+1)
      C=D(J+4)
      ZX=1./ZR
C
      ZW=(-C*ZX+2.)*ZX+E
      ZWH=SQRT(ZW)
      Z=ZR*ZWH
C  COMPUTE PHASE
      ZP=Z+D(J+15)
      IF(C.EQ.0)GOTO 30
      ZRMC=ZR-C
      IF(E.EQ.0)GOTO 70
C
C  CASE OF C.GT.0 AND E.GT.0
      FK=D(J+2)
      ZRK=ZR*FK
      ZALP=Z+ZRK
C  LOG TERM TO PHASE
      ZB=FK*ZALP
      IF(ABS(ZB).GT.ACJWBK)GOTO 10
      ZB=-ZB
      ZP = ((((.2*ZB+.25)*ZB+.33333333)*ZB+.5)*ZB+1.)*ZALP + ZP
      GOTO 20
   10 ZP = LOG(ZB+1.)*D(J+3) + ZP
C  ARCTAN TERM
   20 ZA1=1./(ZR*D(J+7))
      ZS=D(J+5)*(Z-FK*ZRMC)*ZA1
      CK=D(J+10)
      ZG=(CK*Z+ZRMC)*ZA1
      ZP =  D(J+6)*(0.,-1.)*LOG((0.,1.)*ZS+ZG) + ZP +
C  CAP PHI TERM
     *   (5.*ZRMC/(Z*Z)-(Z*D(J+9)+ZRK*D(J+8)+CK)/(ZALP*D(J+7)))/(24.*Z)
      ZCC = ((D(J+14)*ZX+D(J+13))*ZX+D(J+12))*ZX+D(J+11)
      GO TO 80
C
   30 IF(E.EQ.0)GOTO 60
C
C  CASE OF C.EQ.0 AND E.GT.0
      FK=D(J+2)
      ZALP=FK*ZR+Z
      ZB=FK*ZALP
      IF(ABS(ZB).GT.ACJWBK)GOTO 40
      ZB=-ZB
      ZP = ((((.2*ZB+.25)*ZB+.33333333)*ZB+.5)*ZB+1.)*ZALP + ZP
      GOTO 50
   40 ZP = D(J+3)*LOG(1.+ZB) + ZP
   50 ZP = 1./(4.*ZALP) + (5.*ZR/(Z*Z)-2.*(Z+ZALP)/ZALP)/(24.*Z) + ZP
      ZCC = -4.*E-3.*ZX
      GO TO 80
C
C  CASE OF C.EQ.0 AND E.EQ.0
   60 ZP = (.093750*ZX+1.)*Z + ZP
      ZCC = -.0234375*ZX
      GO TO 90
C
C  CASE OF E.EQ.0 AND C.GT.0
   70 ZP =  D(J+6) * (0.,-1.) * LOG(ZRMC*ZX+(0.,1.)*D(J+5)*Z*ZX)
     *     -(3.*ZR+C)/((ZRMC+ZR)*Z*24.) + Z + ZP
      ZCC = ((D(J+14)*ZX+D(J+13))*ZX-3.)*ZX
C
C  COMPLETE CALCULATION OF ZPHI
   80 ZCC = (ZX/ZW)**3 * ZCC * .0625
   90 ZAI=(1.-ZCC)/SQRT(ZWH)
      ZPI=ZP
      RETURN
      END
C***************************************************************
C
      SUBROUTINE ZTHETA(I,ZR,ZTA,ZTDA,ZTP)
C
C  CALCULATES THETA AND THETAD FOR CHANNEL I AND COMPLEX ZR
C     THETA = ZTA*CEXP(ZTP)
C     THETAD = ZTDA**CEXP(ZTP)
C     ZTP = FNUI*LOG(ZR) - ZR/FNUI
C
C     IMPLICIT REAL*8(A-H,O-Y)
      IMPLICIT COMPLEX(Z)
C
      COMMON/CTHET/BB( 500,  50),BG( 500,  50),MSUM( 500)
C     E=BG(I,1)
C
      FNUI=BB(I,1)
      Z = ZR/FNUI
      ZY=1./(Z+Z)
      ZAS=1.
      ZS=BB(I,2)
      ZX=0.
      DO 10 M=3,MSUM(I)
       ZAS=ZAS*ZY
       ZX=ZX+BG(I,M)*ZAS
   10  ZS=ZS+BB(I,M)*ZAS
C
      ZDLR=LOG(ZR)
      ZTP = FNUI*ZDLR-Z
      ZTA=ZS
      ZTDA=((BG(I,2)*ZR+ZDLR)*ZS+ZX)*BG(I,1)
C
      RETURN
      END
C***************************************************************
C
      SUBROUTINE ZSLAG(I,J,NLAG,LIJ,ZS3,ZSD3)
C
C  CALCULATES S INTEGRALS USING LAGUERRE QUADRATURE
C
C     IMPLICIT REAL*8(A-H,O-Y)
      IMPLICIT COMPLEX(Z)
C
      COMMON/CBLK/XLAG(30),WLAG(30),XLEG(15),WLEG(15)
      COMMON/CPOINT/RZERO,RONE,RTWO,H,KP0,KP1,KP2
      COMMON/CHAN/ECH( 500),EPS( 500),FKNU( 500),CC( 500)
     1  ,RINF( 500),ITARG( 500),LLCH( 500),NCHF,NCHOP,NCHOP1
C
      ZFK=CMPLX(FKNU(I),1./FKNU(J))
      X=RTWO
      B=SQRT(2.*X)
      ZB=(0.,1.)/B
      ZG=.5*ZFK*B
C
      ZA1=1./ZG
      ZA2=1.+ZG
      ZA2=ZA2*ZA2
      ZA3=ZB*ZG
C
      NS=NLAG/2
      M=NS*(NS-1)
      N1=M+1
      N2=M+NLAG
C
      ZS3=0.
      ZSD3=0.
      DO 10 N=N1,N2
      U=XLAG(N)
      ZET=ZA1*(SQRT(ZA2+ZA3*U)-1.)
      ZMUM=-.5*ZET/(ZB*(1.+ZET*ZG))
      ZET=ZET*ZET
      ZR=ZET*X
      CALL ZPHI(I,ZR,ZAI,ZPI)
      CALL ZTHETA(J,ZR,ZTAJ,ZTDAJ,ZTPJ)
      ZB1 = (ZR**(-LIJ))*WLAG(N)*ZMUM*EXP((0.,1.)*ZPI+ZTPJ+U)*ZAI
      ZS3=ZS3+ZB1*ZTAJ
   10 ZSD3=ZSD3+ZB1*ZTDAJ
C
      RETURN
      END
C***************************************************************
C
      SUBROUTINE TLAG(I,J,NLAG,LIJ,T1,T2,T3)
C
C  CALCULATES T INTEGRALS FOR CHANNELS I,J USING
C  LAGUERRE QUADRATURE WITH NLAG POINTS
C  THE INTEGRALS ARE -
C    T1 FOR (THETAI,THETAJ)
C    T2 FOR (THETADI,THETAJ)
C    T3 FOR (THETAI,THETADJ)
C
C     IMPLICIT REAL*8(A-H,O-Z)
C
      COMMON/CTHET/BB( 500,  50),BG( 500,  50),MSUM( 500)
      COMMON/CBLK/XLAG(30),WLAG(30),XLEG(15),WLEG(15)
      COMMON/CPOINT/RZERO,RONE,RTWO,H,KP0,KP1,KP2
C
C  CAP 1/K
      F1 = 1./(1./BB(I,1)+1./BB(J,1))
C  INITIALISE FOR QUADRATURE
      NS=NLAG/2
      M=NS*(NS-1)
      N1=M+1
      N2=M+NLAG
      T1=0.
      T2=0.
      T3=0.
C  START QUADRATURE
      DO 40 N=N1,N2
       R = XLAG(N)*F1+RTWO
C  CALCULATE THETA FUNCTIONS
       CALL TANDTD(R,I,TI,TDI,TPI)
       CALL TANDTD(R,J,TJ,TDJ,TPJ)
C  ADD TO SUM
C      A1=(R**(-LIJ))*WLAG(N)*EXP(XLAG(N)+TPI+TPJ)
C++ VAX MOD -- REMOVED '97MAR5 WHEN TIDYING UP A1:
       A1 = EXP(XLAG(N)-LIJ*LOG(R)+TPI+TPJ) * WLAG(N)
       T1=T1+TI*A1*TJ
       T2=T2+TDI*A1*TJ
   40  T3=T3+TI*A1*TDJ
      T1=T1*F1
      T2=T2*F1
      T3=T3*F1
C
      RETURN
      END
C
C***************************************************************
C
      SUBROUTINE TANDTD(R,I,TA,TDA,TP)
C
C  CALCULATES THETA AND THETAD =THETA DOT FOR  R REAL
C     THETA = TA*EXP(TP)
C     THETAD = TDA*EXP(TP)
C     TP = FNU*LOG(R) - R/FNU
C
C     IMPLICIT REAL*8(A-H,O-Z)
      COMMON/CTHET/BB( 500,  50),BG( 500,  50),MSUM( 500)
C
      MI=MSUM(I)
      S=BB(I,2)
      CX=0.
      FNU=BB(I,1)
C
      X=2.*R/FNU
      Y=1./X
      AS=1.
      DO 10 L=3,MI
      AS=AS*Y
      S=S+BB(I,L)*AS
   10 CX=CX+BG(I,L)*AS
C
      DLR=LOG(R)
      TP = FNU*DLR-.5*X
      TA=S
      TDA = ((BG(I,2)*R+DLR)*S+CX)*BG(I,1)
C
      RETURN
      END
C***************************************************************
C
      COMPLEX FUNCTION ZPLAG(I,J,NLAG,LIJ)
C
C  CALCULATES P INTEGRALS USING LAGUERRE QUADRATURE
C
C     IMPLICIT REAL*8(A-H,O-Y)
      IMPLICIT COMPLEX(Z)
C
      COMMON/CBLK/XLAG(30),WLAG(30),XLEG(15),WLEG(15)
      COMMON/CPOINT/RZERO,RONE,RTWO,H,KP0,KP1,KP2
      COMMON/CHAN/ECH( 500),EPS( 500),FKNU( 500),CC( 500)
     1  ,RINF( 500),ITARG( 500),LLCH( 500),NCHF,NCHOP,NCHOP1
      COMMON/CACC/AC,ACNUM,ACJWBK,ACZP,LACC
C
      FK=FKNU(I)+FKNU(J)
      B=SQRT(8.*RTWO)
      G=FK*.125*B
      IF(FK.GT.0)THEN
      GM=1./G
      G2=1.+G
      G2=G2*G2
      ENDIF
      ZB=(0.,1.)/B
C
      NS=NLAG/2
      M=NS*(NS-1)
C
      ZPLAG=0.
      DO 30 N=M+1,M+NLAG
      U=XLAG(N)
      A1=FK*U
      IF(A1.LE.ACZP)THEN
      ZET=1.+.5*ZB*U
      ELSE
      ZET=(SQRT(G2+ZB*G*U)-1.)*GM
      ENDIF
      ZMU=-8.*ZB*(G+1./ZET)
      ZET=ZET*ZET
      ZR=ZET*RTWO
      CALL ZPHI(I,ZR,ZAI,ZPI)
      CALL ZPHI(J,ZR,ZAJ,ZPJ)
   30 ZPLAG=ZPLAG+ZAI*ZAJ*EXP((0.,1.)*(ZPI+ZPJ)+U)*
     1 (ZR**(-LIJ))*WLAG(N)/ZMU
C
      RETURN
      END
C***************************************************************
C
      COMPLEX FUNCTION ZQLAG(I,J,NLAG,LIJ)
C
C  CALCULATES Q INTEGRALS USING LAGUERRE QUADRATURE
C
C     IMPLICIT REAL*8(A-H,O-Y)
      IMPLICIT COMPLEX(Z)
C
      COMMON/CBLK/XLAG(30),WLAG(30),XLEG(15),WLEG(15)
      COMMON/CPOINT/RZERO,RONE,RTWO,H,KP0,KP1,KP2
      COMMON/CHAN/ECH( 500),EPS( 500),FKNU( 500),CC( 500)
     1  ,RINF( 500),ITARG( 500),LLCH( 500),NCHF,NCHOP,NCHOP1
C
      ZMUM = (0.,1.)/(FKNU(I)-FKNU(J))
      NS=NLAG/2
      M=NS*(NS-1)
C
      ZQLAG=0.
      DO 30 N=M+1,M+NLAG
       ZR = XLAG(N)*ZMUM+RTWO
       CALL ZPHI(I,ZR,ZAI,ZPI)
       CALL ZPHI(J,ZR,ZAJ,ZPJ)
   30  ZQLAG=ZQLAG+ZAI*ZAJ*EXP((0.,1.)*(ZPI-ZPJ)+XLAG(N))*
     1  (ZR**(-LIJ))*WLAG(N)
      ZQLAG=ZQLAG*ZMUM
C
      RETURN
      END
C***************************************************************
C
      COMPLEX FUNCTION ZQLEG(I,J,NLEG,LIJ)
C
C  CALCULATES Q INTEGRALS USING LEGENDRE QUADRATURE
C
C     IMPLICIT REAL*8(A-H,O-Y)
      IMPLICIT COMPLEX(Z)
C
      COMMON/CBLK/XLAG(30),WLAG(30),XLEG(15),WLEG(15)
      COMMON/CPOINT/RZERO,RONE,RTWO,H,KP0,KP1,KP2
      COMMON/CHAN/ECH( 500),EPS( 500),FKNU( 500),CC( 500)
     1  ,RINF( 500),ITARG( 500),LLCH( 500),NCHF,NCHOP,NCHOP1
C
      FK=FKNU(I)-FKNU(J)
      X=RTWO
      ZA=1.
      IF(FK.NE.0.) ZA=(0.,1.)/(FK*X+1.)
C
      NS=NLEG/2
      J1=(NS*(NS-1))/2
      J2=J1+NS
      J1=J1+1
C
      ZQLEG=0.
      DO 10 JJ=J1,J2
      U=1.+XLEG(JJ)
      W=1.-XLEG(JJ)
      ZR1 = (ZA*W/U+1.)*X
      ZR2 = (ZA*U/W+1.)*X
      CALL ZPHI(I,ZR1,ZAI,ZPI)
      CALL ZPHI(J,ZR1,ZAJ,ZPJ)
      ZF = EXP((0.,1.)*(ZPI-ZPJ))*(ZR1**(-LIJ))*ZAI*ZAJ/(U*U)
      CALL ZPHI(I,ZR2,ZAI,ZPI)
      CALL ZPHI(J,ZR2,ZAJ,ZPJ)
      ZF = EXP((0.,1.)*(ZPI-ZPJ))*(ZR2**(-LIJ))*ZAI*ZAJ/(W*W)+ZF
   10 ZQLEG=ZQLEG+ZF*WLEG(JJ)
      ZQLEG=2.*X*ZA*ZQLEG
C
      RETURN
      END
C***************************************************************
C
      SUBROUTINE INV(A,N)
C
C  MATRIX INVERSION WITH FULL PIVOTING
C
C     IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION IPIVOT( 500),A( 500, 500),INDEX( 500,2),PIVOT( 500)
C
C
C  INITIALIZATION
      DO 20 J=1,N
   20  IPIVOT(J)=0
      DO 550 I=1,N
C  SEARCH FOR PIVOT ELEMENT
       AMAX=0.
       DO 105 J=1,N
        IF(IPIVOT(J).EQ.1) GO TO 105
        DO 100 K=1,N
         IF(IPIVOT(K)-1) 80,100,740
   80    T=ABS(A(J,K))
         IF(T.LE.AMAX) GO TO 100
         IROW=J
         ICOL=K
         AMAX=T
  100    CONTINUE
  105   CONTINUE
       IPIVOT(ICOL)=IPIVOT(ICOL)+1
C  INTERCHANGE ROWS TO PUT PIVOT ELEMENT ON DIAGONAL
       IF(IROW.EQ.ICOL) GO TO 260
       DO 200 L=1,N
        SWAP=A(IROW,L)
        A(IROW,L)=A(ICOL,L)
  200   A(ICOL,L)=SWAP
  260  INDEX(I,1)=IROW
       INDEX(I,2)=ICOL
       PIVOT(I)=A(ICOL,ICOL)
C  DIVIDE PIVOT ROW BY PIVOT ELEMENT
       A(ICOL,ICOL)=1.
       DO 350 L=1,N
  350   A(ICOL,L)=A(ICOL,L)/PIVOT(I)
C  REDUCE NON-PIVOT ROWS
       DO 550 K=1,N
       IF(K.EQ.ICOL) GO TO 550
       T=A(K,ICOL)
       A(K,ICOL)=0.
       DO 450 L=1,N
  450   A(K,L)=A(K,L)-A(ICOL,L)*T
  550  CONTINUE
C  INTERCHANGE COLUMNS
      DO 710 L=N,1,-1
       IF(INDEX(L,1).EQ.INDEX(L,2)) GO TO 710
       JROW=INDEX(L,1)
       JCOL=INDEX(L,2)
       DO 700 K=1,N
        SWAP=A(K,JROW)
        A(K,JROW)=A(K,JCOL)
  700   A(K,JCOL)=SWAP
  710  CONTINUE
  740 RETURN
      END
C
C**************************************************************
C
      SUBROUTINE MATSQ(A,B,N)
C
C  CALCULATES B=A*A
C
C     IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A( 500, 500),B( 500, 500)
C
      DO 2 I=1,N
      DO 2 J=1,N
      S=0.
      DO 1 K=1,N
    1 S=S+A(I,K)*A(K,J)
    2 B(I,J)=S
C
      RETURN
      END
C
C***************************************************************
C
      SUBROUTINE CEIGEN(A,VAL,N,DELTA)
C
C  DIAGONALISATION OF COMPLEX SYMMETRIC MATRIX (N<3 ERROR 1998-99May20)
C  USING METHOD DESCRIBED BY M.J.SEATON, COMPUTER JOURNAL 12 (1969) 156
C
C     IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX A(  20,  20),X(  20,  20),VAL(  20),C,S,H,P,Q,CIM
      DATA CIM/(0.,1.)/, IROTM/20/
C                        ===== SET
C
      SQ(H)=H*CONJG(H)
C
C
      C = (1.,0.)
      H = (0.,0.)
      IF(N-2) 710,720,500
C
C INITIAL D1,D2 AND X
C
  500 D1=0.
      D2=0.
      DO 520 I=1,N
       X(I,I)=(1.,0.)
       DO 510 J=1,I-1
        X(I,J)=(0.,0.)
        X(J,I)=(0.,0.)
  510   D2=D2+SQ(A(J,I))
       D1=D1+SQ(A(I,I))
  520  H=H+A(I,I)
C
      D1 = (N-1) * (D1-SQ(H)/N + D2+D2) * DELTA*DELTA/(N+N+2.)
C
C BEGIN ROTATIONS
C
      DO 1010 IROT=1,IROTM
      DO 1000 IP=1,N-1
      DO 1000 IQ=IP+1,N
C
C ROTATION CONSTANTS
C
      Q=A(IP,IQ)
      P=(A(IP,IP)-A(IQ,IQ))*.5
      FL=ABS(P*P+Q*Q)
      BETA=0.5*LOG(SQ(P-CIM*Q)/FL)
      T=(CONJG(P)*Q+P*CONJG(Q))/FL
      D=(SQ(P)-SQ(Q))/FL
      U=-0.25*ATAN2(T,D)
C
      T=0.
      D=0.
      DO 560 I=1,N
      IF(I.EQ.IP.OR.I.EQ.IQ) GO TO 560
      D=D+SQ(A(IP,I))+SQ(A(IQ,I))
      T=T+SQ(A(IP,I)+CIM*A(IQ,I))
  560 CONTINUE
      T=D-T
      FN= SQRT(D*D-T*T)
      GAMMA=LOG((D+T)/FN)
C
C ITERATION FOR V
C
      V0=-0.5*(BETA+GAMMA)
      DO 570 ITERV=1,100
      V=V0-(FL*SINH(2.*(V0+BETA))+FN*SINH(V0+GAMMA))/
     1  (2.*FL*COSH(2.*(V0+BETA))+FN*COSH(V0+GAMMA))
      IF(ABS(V-V0).LT.DELTA*0.1) GO TO 580
  570 V0=V
      WRITE(6,5010) DELTA
C
C NEW A,X AND D2
C
  580 V=0.5 0*V
      S= SIN(U)* COSH(V)+CIM* COS(U)* SINH(V)
      C= COS(U)* COSH(V)-CIM* SIN(U)* SINH(V)
      H=(S*P+C*Q)*S*2.
      A(IP,IP)=A(IP,IP)-H
      A(IQ,IQ)=A(IQ,IQ)+H
      A(IP,IQ)=(C*C-S*S)*Q+C*S*P*2.
      D2=D2-SQ(A(IQ,IP))+SQ(A(IP,IQ))
      A(IQ,IP)=A(IP,IQ)
C
      DO 590 I=1,N
      H=X(I,IP)
      X(I,IP)=H*C-X(I,IQ)*S
      X(I,IQ)=H*S+X(I,IQ)*C
      IF(I.EQ. IP.OR.I.EQ.IQ)GO TO 590
      H=A(IP,I)
      A(IP,I)=C*H-A(IQ,I)*S
      A(IQ,I)=S*H+A(IQ,I)*C
      D2=D2-SQ(A(I,IP))-SQ(A(I,IQ))+SQ(A(IP,I))+SQ(A(IQ,I))
      A(I,IP)=A(IP,I)
      A(I,IQ)=A(IQ,I)
  590 CONTINUE
C
C TEST CONVERGENCE
      IF(D2.LT.D1) GO TO 610
C
C END ROTATIONS
C
 1000 CONTINUE
C
C  RECALCULATE D2
C
      D2=0.
      DO 1001 J=2,N
       DO 1001 I=1,J-1
 1001   D2=D2+SQ(A(I,J))
 1010 CONTINUE
      WRITE(6,5020)IROTM,DELTA
C
C EIGENVALUES AND EIGENVECTORS
C
  610 DO 620 I=1,N
       VAL(I)=A(I,I)
       DO 620 J=1,N
  620   A(J,I)=X(J,I)
      RETURN
C
C N=1 AND N=2
C
  720 Q=A(1,2)
      P = (A(1,1)-A(2,2))*.5
      FL= ABS(P*P+Q*Q)
      V=-0.25*LOG(SQ(P-CIM*Q)/FL)
      T= CONJG(P)*Q+P* CONJG(Q)
      U= -0.25 * ATAN2(T/FL,(SQ(P)-SQ(Q))/FL)
      S= SIN(U)* COSH(V)+CIM* COS(U)* SINH(V)
      C= COS(U)* COSH(V)-CIM* SIN(U)* SINH(V)
      H=(S*P+C*Q)*S*2.
      VAL(2)=A(2,2)+H
      A(1,2)=S
      A(2,1)=-S
      A(2,2)=C
  710 VAL(1)=A(1,1)-H
      A(1,1)=C
C
 5010 FORMAT(//10X,'*** SUBROUTINE CEIGEN ***'
     + //' WARNING - NO CONVERGENCE IN ITERATIONS FOR V'/
     +  ' ACCURACY PARAMETER DELTA =',1P,E9.2/
     +  ' NEXT VALUE OF AVERAGED OMEGA MAY BE INCORRECT'//)
 5020 FORMAT(//11X,'*** SUBROUTINE CEIGEN ***'
     * //' NO CONVERGENCE FOR IROTM =',I3/' DELTA =',E11.3/
     +  ' NEXT VALUE OF AVERAGED OMEGA MAY BE INCORRECT'//)
C
      RETURN
      END
C
C**********************************************************
      FUNCTION WWW(IC1,IC2,KCH,KADD,LADD,JTOT)
C     PAIR COUPLING EQUIVALENT TO WSQ IN LS COUPLING - RUB'96APR-M,97OCT
C     IMPLICIT REAL*8(C,W)
      DIMENSION KCH( 500,3)
C     TARGET ORBITAL QN:
      I1 = KCH(IC1,3)
      I2 = KCH(IC2,3)
C     CHANNEL ORBITAL QN:
      L1 = KCH(IC1,2)+LADD
      L2 = KCH(IC2,2)+LADD
C     CHANNEL PAIR COUPLING QN:
      KK = KCH(IC1,1)
      IF (KK.LT.0) THEN
        KK=JTOT
        C = 1.0
      ELSE
        WWW = 0.
        IF (KCH(IC2,1).NE.KK) RETURN
        KK=KADD+KK
        C = (JTOT+1)*.5/(KK+1)
      ENDIF
      WWW = WSQ(I1,I2,L1,L2,KK) * C
      RETURN
      END
C
      FUNCTION WSQ(A,B,C,D,F)
C
C  CALCULATES 3*(2*F+1)*W(A,B,C,D,1,F)**2 WHERE W IS A RACAH COEFFICIENT
C                REFORMULATED FOR TWICE THE VALUES OF A-F -- RUB'96MAR24
      IMPLICIT INTEGER(A-F)
C     IMPLICIT REAL*8(W)
C
      IF(B-A)10,20,30
C
   10 IF(D-C)11,12,13
   20 IF(D-C)21,22,23
   30 IF(D-C)31,32,33
C
   11 WSQ = ((F+1)*(B+D+F+4)*(B+D+F+6)*(B+D-F+2)*(B+D-F+4)*3)/
     *    REAL((B+1)*(B+2)*(B+3)*(D+1)*(D+2)*(D+3)*16)
      RETURN
   12 WSQ = ((F+1)*(B+D+F+4)*(-B+D+F)*(B-D+F+2)*(B+D-F+2)*3)/
     *    REAL((B+1)*(B+2)*(B+3)*D*(D+1)*(D+2)*8)
      RETURN
   13 WSQ = ((F+1)*(-B+D+F-2)*(-B+D+F)*(B-D+F+2)*(B-D+F+4)*3)/
     *    REAL((B+1)*(B+2)*(B+3)*(D-1)*D*(D+1)*16)
      RETURN
C
   21 WSQ = ((F+1)*(B+D+F+4)*(-B+D+F+2)*(B-D+F)*(B+D-F+2)*3)/
     *    REAL(B*(B+1)*(B+2)*(D+1)*(D+2)*(D+3)*8)
      RETURN
   22 WSQ = ((B*(B+2)+D*(D+2)-F*(F+2)))**2 * (F+1)*3/
     *    REAL(B*(B+1)*(B+2)*D*(D+1)*(D+2)*4)
C      ((2)F+1)**2 CORR'96APR19-20.
      RETURN
   23 WSQ = ((F+1)*(B+D+F+2)*(-B+D+F)*(B-D+F+2)*(B+D-F)*3)/
     *    REAL(B*(B+1)*(B+2)*(D-1)*D*(D+1)*8)
      RETURN
C
   31 WSQ = ((F+1)*(-B+D+F+2)*(-B+D+F+4)*(B-D+F-2)*(B-D+F)*3)/
     *    REAL((B-1)*B*(B+1)*(D+1)*(D+2)*(D+3)*16)
      RETURN
   32 WSQ = ((F+1)*(B+D+F+2)*(-B+D+F+2)*(B-D+F)*(B+D-F)*3)/
     *    REAL((B-1)*B*(B+1)*D*(D+1)*(D+2)*8)
      RETURN
   33 WSQ = ((F+1)*(B+D+F)*(B+D+F+2)*(B+D-F-2)*(B+D-F)*3)/
     *    REAL((B-1)*B*(B+1)*(D-1)*D*(D+1)*16)
      RETURN
C
      END
C***************************************************************
C
      SUBROUTINE MATMUL(A,B,C,N)
C
C  CALCULATES C=A*B
C
C     IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A( 500, 500),B( 500, 500),C( 500, 500)
C
      DO 2 I=1,N
      DO 2 J=1,N
      S=0.
      DO 1 K=1,N
    1 S=S+A(I,K)*B(K,J)
    2 C(I,J)=S
C
      RETURN
      END
C*******************************************************************
      SUBROUTINE OUTJJM(ETOT,A,B,LIT)
C
C    TO OUTPUT COLLISION MATRICES FOR JJOM ON UNIT 11 (UNFORMATTED)
C instructions by KAB Feb 91.
C in MAIN, insert after C     IOPT1 = 1 FOR ALL SLPI CASES
C     =10,11 FOR JJOM WITH ALGEBRAIC RECOUPLING WITHOUT,WITH TOP-UP
C      10,11 REDUNDANT, NOW SPECIFYING T OUTPUT. WE/JAT'93APR2-3.
C     =12,13 FOR JJOM WITH TERM COUPLING (COUPLING COEFFS READ IN JJOM)
C     =14,15 KAB'S BINARY MODE CODED HERE RUB'97DEC BY WE.
C in MAIN, insert after statement READ(5,*)IOPT1:
C     IF(IOPT1.GE.10) CALL OUTJJ(-1)
C in MAIN, insert after statement labelled 3100:
C     IF(IOPT1.GE.10) CALL OUTJJ(IOPT1)
C in REACT, for K insert before  IF(ITARG(NCHOP).EQ.ITARG(1)) RETURN
C     IF(IOPT1.GE.10) CALL OUTJJM(ETOT,RK,RK,NCHOP,'K  ')
C     while CHI matrix output requires IMESH.EQ.2
C new subroutines OUTJJM and OUTJJ:
C    ETOT * RY = ELECTRON ENERGY
C    RK(I,J) = K-MATRIX
C    NCHF (NCHOP) = NUMBER OF (OPEN) CHANNELS
C    /CINPUT/ CONTAINS BASIC COLLISION DATA.
C
C    STORE PARTIAL WAVES (ILT(ISLM),IST(ISLM)) AND ENERGIES (EK(NEN))
C    IN /STOJJ/ WITH NO DUPLICATION, FOR WRITING OUT TO JJOM IN
C    SUBROUTINE OUTJJ. KEEP A RUNNING KOUNT TOTAL OF THE NUMBER OF
C    OUTPUT MATRIX ELEMENTS.
C
C     IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(IEX=0)
      CHARACTER*3 LIT
      COMMON/CHAN/ECH( 500),EPS( 500),FKNU( 500),CC( 500)
     1  ,RINF( 500),ITARG( 500),LLCH( 500),NCHF,NCHOP,NCHOP1
      COMMON/CINPUT/ CF( 500, 500,  20),COEFF(3,  25),ENAT(0: 109),
     5 WMAT( 500,7000),VALUE(7000),RA,BSTO,AZAZ,
     1 LAMAX,LRANG2,LRGL2,MNP2,NAST,NCHAN,NELC,NPTY2,NSPN2,NZED,MORE2,
     2 ISAT( 109),LAT( 109),NCONAT( 109)
      COMMON/STOJJ/EK(5900),QNMAX,NEN,ILT(100),IST(100),ISLM,NAOP,KOUNT
      DIMENSION A( 500, 500),B( 500, 500)
C
C    CONVERT ENERGIES TO UNSCALED RYDBERG RELATIVE TO GROUND-STATE
C
      E=ETOT*AZAZ
C
C    STORE OR LOCATE PARTIAL-WAVE AND ENERGY IN /STOJJ/
C
      IF(ISLM.EQ.0) GO TO 4
      DO 3 I=1,ISLM
      IF(LRGL2.EQ.ILT(I).AND.NSPN2.EQ.IST(I)) GO TO 5
    3 CONTINUE
      IF(ISLM.LT.100) GO TO 4
        WRITE(6,*)' STOP. MAKE ARRAYS ISX LONGER THAN 100!'
        RETURN
    4 ISLM=ISLM+1
      ILT(ISLM)=LRGL2
      IST(ISLM)=NSPN2
C
    5 IF(NEN.EQ.0) GO TO 9
      DO 6 I=1,NEN
C      IF(ABS(E-EK(I)).LT.TINY) GO TO 10 -- HES'93MAR10:
       IE = I
       IF(EK(I).EQ.E) GO TO 10
    6  CONTINUE
      IF(NEN-5900) 9,7,8
    7   WRITE(6,*)' STGF LIMITED TO 5900 ENERGIES, MAKE EST BIGGER'
        NEN=NEN+1
    8   RETURN
    9 NEN=NEN+1
      IE = NEN
      EK(NEN)=E
C
C
C    WRITE OUT MATRIX K, -T OR CHI
C
   10 M = (1-2*NPTY2)*(NSPN2*1000+LRGL2)
      N = NCHOP
      IF (LIT.EQ.'CHI') N=NCHF
      DO 15 J=1,N
       NAOP=MAX(ITARG(J),NAOP)
   15  CONTINUE
      IF(LIT.EQ.'BIN') THEN
       KNUM = ((NCHOP+1)*NCHOP)/2
       WRITE(11) IE, M, KNUM
       KOUNT = KNUM + KOUNT
       WRITE(11) ((A(I,J),I=1,J),J=1,N), ((B(I,J),I=1,J),J=1,N)
      ELSE
       BB=0.
       DO 18 J=1,N
        IF(J.GT.NCHOP.AND.FKNU(J)+.2E-4.GT.MAX(QNMAX,LLCH(J)+.5))GOTO 20
        KAB1=ITARG(J)*1000+LLCH(J)
        DO 16 I=1,J
          KAB2=ITARG(I)*1000+LLCH(I)
          IF(LIT.NE.'K  ') BB = B(I,J)
          KOUNT=KOUNT+1
          K=MOD(KOUNT,100000)
   16     WRITE(11,1013) M,KAB1,KAB2, A(I,J),BB, NZED,IEX,E,NELC,LIT, K
C  !!     BUT HES'93MAR10   BB=B(J,I);A(J,I),=SYMM! FOR UNKNOWN REASONS
   18   CONTINUE
      ENDIF
C1013 FORMAT(I5,2I6,E16.8,E17.8,I3,I2,F12.8,I3,' K   ',I5) - WE'92SEP18:
 1013 FORMAT(I5,2I6,1P,2E16.8,I4,I2,0P,F12.7,I3,1X,A3,I6)
C1013 FORMAT(I5,2I6,1P,2E16.8,I4,I2,0P,F12.8,I3,1X,A3,I6) CATAST'93MAY1!
   20 RETURN
      END
C*******************************************************************
      SUBROUTINE OUTJJ(IOPT1,NRANG2)
C
C    TO SUPPLY NON-K-MATRIX DATA FOR JJOM ON UNIT 12 (FORMATTED).
C
C    IOPT1=-1 FOR INITIALISATION
C     =10,11 FOR JJOM WITH ALGEBRAIC RECOUPLING WITHOUT,(WITH TOP-UP)
C      10,11 REDUNDANT, NOW SPECIFYING T MATRIX OUTPUT. WE/JAT'93APR2-3.
C     =12,13 FOR JJOM WITH TERM COUPLING (COUPLING COEFFS READ IN STGF)
C     =14,15 KAB'S JBIN BINARY FORM; 11, 13 AND 15 SETTING ITOP=2!
C     /STOJJ/ CONTAINS DATA FOR JJDAT PASSED FROM SR OUTJJM
C
C    /CINPUT/ CONTAINS BASIC COLLISION DATA.
C
C     IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(JOUT=12,NOX=0,IFACT=105,ntc= 109* 109)
      PARAMETER(IREFL=1,IPART=0, IB=7,IW=8)
      CHARACTER*1 PARITY(-1:1),LSPECT(0:7), LIT*3
      COMMON/CDEC/ARAD(5886),ARDEC( 109),IRDEC,IPAR( 109),KCH( 500,3)
      COMMON/CINPUT/ CF( 500, 500,  20),COEFF(3,  25),ENAT(0: 109),
     5 WMAT( 500,7000),VALUE(7000),RA,BSTO,AZAZ,
     1 LAMAX,LRANG2,LRGL2,MNP2,NAST,NCHAN,NELC,NPTY2,NSPN2,NZED,MORE2,
     2 ISAT( 109),LAT( 109),NCONAT( 109)
      COMMON/CMESH/EMAX,EMIN,DEOPEN,DQN,QNMAX,EMESH(9000),IMESH
      COMMON/CNTRL/IPRINT,IRAD,IPERT,KP2X
      COMMON/STOJJ/EK(5900),QNMXX,NEN,ILT(100),IST(100),ISLM,NAOP,KCOUNT
      DIMENSION FJ(  25), IDF(ntc),ITF(ntc),FCF(ntc)
      DATA LSPECT/'S','P','D','F','G','H','I','X'/, PARITY/'?','e','o'/
      DATA J2F,N2F/2*0/
C
      IF(IOPT1.EQ.-1) THEN
      NEN=0
      ISLM=0
      NAOP=0
      KCOUNT=0
      QNMXX=QNMAX
      RETURN
      ENDIF
C
      JCHAN=0
      IF(NAOP.EQ.0) GO TO 42
      IF(NEN.EQ.0.OR.ISLM.EQ.0.OR.KCOUNT.EQ.0) GO TO 42
      NEK=NEN
      IF(NEN.GT.5900) NEK=5900
C
C
C    WRITE OUT JJOM INPUT CARDS A1-A2
C
      ITOP=0
      IF(MOD(IOPT1,2).NE.0) ITOP=2
      MAT=IMESH
      IF(IOPT1.EQ.10.OR.IOPT1.EQ.11) MAT=2
      IONQ=NZED-NELC
      WRITE(JOUT,'(2X,I2,2I3,10I5)')
     *      NZED,NELC,NAOP,IREFL,IW,IB,ITOP,MAT,IONQ,IPART,IRDEC,IMESH
C     FIND TO NAOP = INDEX OF HIGHEST TARGET STATE INVOLVED:
C         MNSPN,MXSPN = MINIMUM,MAXIMUM TARGET STATE 2S+1;
C         JCHAN= NUMBER OF CHANNELS IN INTERMEDIATE COUPLING.
      MNSPN=999
      MXSPN=0
      IOJM=0
      DO 8 I=1,NAOP
       MNSPN=MIN(MNSPN,ISAT(I))
       MXSPN=MAX(MXSPN,ISAT(I))
       K=2*LAT(I)+1
       JHIGH=K+ISAT(I)-1
       IOJM=MAX(JHIGH/2+1,IOJM)
       JLOW=ABS(K-ISAT(I))+1
       JCHAN=JCHAN + ((JHIGH+JLOW)/2) * (1+(JHIGH-JLOW)/2)
       K=ISAT(I)
       IF(IPAR(I).EQ.1) K=-K
    8  WRITE(JOUT,1002) I, K,LAT(I),ENAT(I)*AZAZ,
     *    ISAT(I),LSPECT(MIN(LAT(I),7)),PARITY(IPAR(I)),ARDEC(I)
C1002 FORMAT(3I4,F20.8,I3,2A1,E16.8)
 1002 FORMAT(3I4,F20.8,I3,2A1,1P,E16.6)
C
C    CARDS B1-B4
C
      LIT='T/K'
      IF (IMESH.EQ.2) LIT='CHI'
C     WRITE(JOUT,1004) NZED,NELC,LIT, LRANG2,ISLM,NOX, (I,I=0,LRANG2-1)
C1004 FORMAT(' Z=',I2,' N=',I2,
C    *' F.S. COLLISION STRENGTHS FROM ',A3,'-MATRICES'/3I5/(20I3))
      WRITE(JOUT,1004) IMESH,DEOPEN,QNMAX,DQN,IPERT,LIT,IOPT1,
     *LRANG2,ISLM,NOX,(I,I=0,LRANG2-1)
 1004 FORMAT(' IMESH=',I2,' DE=',F8.6,' QNMAX=',F7.3,' DQN=',F8.6,
     * ' IPERT=',I1,1X,A3,'-MTRX OPT=',I2/3I5/(20I3))
      WRITE(JOUT,'(18I4)') (IST(I),ILT(I),I=1,ISLM)
C
C    CARDS C1-C4
C
      WRITE(JOUT,2004) NZED,NELC,RA,NRANG2,LAMAX,KCOUNT,
     *      NEK,NEN,   (K,EK(K),K=1,NEN)
 2004 FORMAT(' Z=',I2,' N=',I2,' RA=',F8.4,' NRANG2.LE.',I2.2,
     *   ' LAMAX=',I2,' COUNT=',I8/ 2I5/(4(I5,F13.7)))
C-OLD*  2I5/(4(I3,F12.8)))
      WRITE(JOUT,'(1P,E13.6,4HEEEE)') EK(NEK)
C
C    FIND ALLOWED TOTAL J VALUES (FJ(JJFSL))
C
      JJFSL=0
      J2=MOD(MNSPN,2)
   22 IZ=MNSPN
      IF(IZ.EQ.2) IZ=0
   23 L2=ABS(J2-IZ)
   24 DO 25 I=1,ISLM
      IF(IZ.EQ.IST(I)-1 .AND.L2.EQ.2*ILT(I)) GO TO 26
   25 CONTINUE
      GO TO 27
   26 L2=L2+2
      IF(L2.LE.J2+IZ) GO TO 24
      IZ=IZ+2
      IF(IZ.LE.MXSPN) GO TO 23
      JJFSL=JJFSL+1
      FJ(JJFSL)=0.5*J2
      J2=J2+2
      GO TO 22
C
C    CARDS E1-E2
C
   27 IFIT=JJFSL+1
C     JPUN=NAOP+1
      JPUN=0
      WRITE(JOUT,'(5H F.S.,2I5,I3)') JJFSL,IFIT,JPUN
      IF(JJFSL.EQ.0) THEN
      WRITE(6,'(/42H NO J VALUES -- SYMMETRIES SLPI INCOMPLETE)')
      GO TO 41
      ENDIF
      WRITE(JOUT,'(7X,13F5.1)') (FJ(J),J=1,JJFSL)
C
C    CARDS F1-F2: TERM COUPLING COEFFICIENTS TO BE READ BY JJOM
C
C-    IF(IOPT1.LT.12)GOTO32 - DELETED MAKING ROOM FOR JAT-T'93APR3-4
C     " A RATHER REDUNDANT CHOICE, BECAUSE OF TERMINATING BLANK!! W.
C      READ TERMCOUPLING COEFFICIENTS, PUT THEM INTO JAJOM INPUT
       WRITE(6,*)
C      BUT SKIP TERM LIST
       READ(5,*,END=32) NTER
       IF(NTER.GT.0) READ(5,'(I4)') (K,I=1,NTER)
   31   READ(5,1106)J2F,N2F
   32   WRITE(JOUT,1106) J2F,N2F
        IF (N2F.EQ.0) GO TO 33
        READ(5,1107) (ITF(I),IDF(I),FCF(I),I=1,N2F)
        WRITE(JOUT,1107) (ITF(I),IDF(I),FCF(I),I=1,N2F)
        WRITE(6,1108)N2F,J2F
        GO TO 31
   33 WRITE(JOUT,1106) NOX,NOX
   41 WRITE(JOUT,'(7HENDDATA)')
 1106 FORMAT(7X,I3,I5)
 1107 FORMAT(5(I3,I2,F9.6))
 1108 FORMAT(I5,' TERM COUPLING COEFFICIENTS PASSED ON FOR 2J =',I3)
C
C     FINALISE OUTPUT TO JJOM, WRITE PROCESSOR DATA FOR JJOM ONTO 13
C
      NEST= NAOP*(NAOP+1)*2*NEK
      IF (NEST.LT.KCOUNT) NEST =KCOUNT
      NEST = NEST/100 +1
      IF (NEST.GT. 9990) NEST = 9990
      WRITE(13,1313) NEST,NEK,JCHAN,NAOP,IOJM,IFACT,NAST
 1313 FORMAT('   11   12   13   -1'/6I5,'    R  CDC',I5,' JOMSPECS')
      WRITE(6,1313) NEST,NEK,JCHAN,NAOP,IOJM,IFACT,NAST
C     NOTE PARAMETRICALLY FIXED IOJM AND IFACT.
   42 WRITE(6,1099) NEN,NAOP,ISLM,KCOUNT,LIT,JCHAN
      RETURN
C
 1099 FORMAT(/I9,' ENERGIES',I6,' STATES',I6,' SL-SYMMETRIES'/I9,1X,
     * A3,'-MATRIX ELEMENTS (TOTAL)',I8,' PAIR-COUPLING CHANNELS'/)
      END