      SUBROUTINE FANOFS(IRHO,ISIG,IRHOP,ISIGP)
C
C     -- adapted for 2-body FS from Coulomb FANO of STG2 by CGX, & WE'04
C
C{-R8}IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER(LL53=  21*2+3,LL75=  21+2, NL1=LL75,NL2=NL1+2,NL3=NL1+5)
c     PARAMETER(KFL1=&L76,     KFL2=&L74, !stg2 AND recupd:
      PARAMETER(KFL1=3*  21+19,KFL2=  21+4, !cor'03Nov12: +2 --> +4
     +          KFL6=120,KFL7=150,KFL8=120,KFL9=40,KFLW=20,
     +          KFLS=12,KFLN=10,KFLV=40)
      PARAMETER(HALF=.5, EPS=1.E-5, ALP4=1.3313E-5)
C
c err DIMENSION RMEDIR(99),RMEEX(99),NBAR(21)
      DIMENSION NBAR(21) !JRELOP(4:5); DATA 2*1 
C
      LOGICAL FAILSD,FAILSE,FAILAD,FAILAE,FREE, J23REL, LET
C
      LOGICAL SMVRSD,SMVRSE,SMVRAD,SMVRAE
      DIMENSION K6SD(KFL6),K7SD(KFL7),K8SD(KFL8),K9SD(KFL9),KWSD(6,KFLW)
     +          ,LDELSD(KFLW,2),SMVRSD(KFL1)
      DIMENSION K6SE(KFL6),K7SE(KFL7),K8SE(KFL8),K9SE(KFL9),KWSE(6,KFLW)
     +          ,LDELSE(KFLW,2),SMVRSE(KFL1)
      DIMENSION K6AD(KFL6),K7AD(KFL7),K8AD(KFL8),K9AD(KFL9),KWAD(6,KFLW)
     +          ,LDELAD(KFLW,2),SMVRAD(KFL1)
      DIMENSION K6AE(KFL6),K7AE(KFL7),K8AE(KFL8),K9AE(KFL9),KWAE(6,KFLW)
     +          ,LDELAE(KFLW,2),SMVRAE(KFL1)
C
      DIMENSION J6PSD(KFLV),J7PSD(KFLV),J8PSD(KFLV),J9PSD(KFLV),
     + JWRDSD(6,KFLW),
     + NBJSD(KFLN),NB6JSD(KFLN),K6CPSD(KFLN),K7CPSD(KFLN),K8CPSD(KFLN),
     + K9CPSD(KFLN),JSM6SD(KFLS),JSM4SD(KFLS,KFLW),JSM5SD(KFLS,KFLW),
     + IN6JSD(KFLW)
      DIMENSION J6PSE(KFLV),J7PSE(KFLV),J8PSE(KFLV),J9PSE(KFLV),
     + JWRDSE(6,KFLW),
     + NBJSE(KFLN),NB6JSE(KFLN),K6CPSE(KFLN),K7CPSE(KFLN),K8CPSE(KFLN),
     + K9CPSE(KFLN),JSM6SE(KFLS),JSM4SE(KFLS,KFLW),JSM5SE(KFLS,KFLW),
     + IN6JSE(KFLW)
      DIMENSION J6PAD(KFLV),J7PAD(KFLV),J8PAD(KFLV),J9PAD(KFLV),
     + JWRDAD(6,KFLW),
     + NBJAD(KFLN),NB6JAD(KFLN),K6CPAD(KFLN),K7CPAD(KFLN),K8CPAD(KFLN),
     + K9CPAD(KFLN),JSM6AD(KFLS),JSM4AD(KFLS,KFLW),JSM5AD(KFLS,KFLW),
     + IN6JAD(KFLW)
      DIMENSION J6PAE(KFLV),J7PAE(KFLV),J8PAE(KFLV),J9PAE(KFLV),
     + JWRDAE(6,KFLW),
     + NBJAE(KFLN),NB6JAE(KFLN),K6CPAE(KFLN),K7CPAE(KFLN),K8CPAE(KFLN),
     + K9CPAE(KFLN),JSM6AE(KFLS),JSM4AE(KFLS,KFLW),JSM5AE(KFLS,KFLW),
     + IN6JAE(KFLW)
C
      COMMON/INFORM/IREAD,IWRITE,IPUNCH
      COMMON/DEBUG/IBUG1,IBUG2,IBUG3,IBUG4,IBUG5,IBUG6,IBUG7,IBUG8,IBUG9
c     COMMON/KRON/IDEL(LL75,LL75)
      COMMON/KRON/IDEL(  21+2,  21+2)
      COMMON/TERMS/ NROWS,ITAB(18),JTAB(18),NTAB(189)
      COMMON/MEDEFN/IHSH,NJ(LL75),LJ(LL75),NOSH1(LL75),NOSH2(LL75),
     1              J1QN1(LL53,3),J1QN2(LL53,3),IJFUL(LL75) ! L53<-L43!
      COMMON/MVALUE/M1,M2,M3,M4,M5,M6,M7,M8,M9,M10,M11,M12,M13,M14,M15,
     1     M16,M17,M18,M19,M20
      COMMON/NJLJ/NRHO,LRHO,NSIG,LSIG,NRHOP,LRHOP,NSIGP,LSIGP
      COMMON/INTERM/J1BAR1(LL75,3),J1BAR2(LL75,3),J1TLD1(3),J1TLD2(3)
      COMMON/COUPLE/NJ1S,NJ23S,J1(KFL1),J2(KFL2,3),J3(KFL2,3),FREE(KFL1)
C
      LOGICAL SUMVAR ! sr.SETDIM,SETTAB,SPRATE,TRDEL,TRIANG,VAR,SQUARE-
      COMMON/ARGU/J6C,J7C,J8C,J9C,JWC,J6(KFL6),J7(KFL7),J8(KFL8),
     + J9(KFL9),KW(6,KFLW),JDEL,LDEL(KFLW,2),SUMVAR(KFL1),MP
C
      COMMON/SUMARG/J6P(KFLV),J7P(KFLV),J8P(KFLV),J9P(KFLV),
     + JWORD(6,KFLW),NLSUM,
     + NBJ(KFLN),NB6J(KFLN),K6CP(KFLN),K7CP(KFLN),K8CP(KFLN),K9CP(KFLN),
     + JSUM6(KFLS),JSUM4(KFLS,KFLW),JSUM5(KFLS,KFLW),INV6J(KFLW)
C     COMMON/XATION/ AMULT(99),BMULT(99),KD1,KD2,KE1,KE2,MULTD,MULTE
      COMMON/DEMULT/DMULT0(NL2,3),DMULT1(NL2,3),EMULT0(NL2,3),
     *       EMULT1(NL2,3),MULDSO,MULESO, KD1,KD2,KE1,KE2
      COMMON/DESSNK/D00SNK(NL2),D11SNK(NL2),E01SNK(NL2),E10SNK(NL2),
     :     MULDSS,MULDSP,MULESS,MULESP
      COMMON/RMEPRD/RMEDIR(99),RMEEX(99)
c
      COMMON/SSKLIM/KDS1,KDS2,KDS3,KDS4,KES1,KES2,KES3,KES4
      COMMON/SSMULT/DMULS0(NL2),DMULS1(NL2),EMULS0(NL2),EMULS1(NL2)
C
C === DETERMINES POTENTIAL MATRIX ELEMENT FOR GIVEN RHO,SIG,RHOP,SIGP
C
C
  301 FORMAT(21H NO POSSIBLE K-VALUES)
  302 FORMAT(66H SPECTATOR QUANTUM NUMBERS NOT DIAGONAL FOR NON-INTERACT
     1ING SHELLS)
  306 FORMAT(37H DIRECT AND EXCHANGE SPIN INTEGRALS =,2F10.6)
  307 FORMAT(7H NJ,LJ:,2(I8,I4),2X,2(I8,I4))
* 308 FORMAT(9H    KD1 =,I4,8H   KD2 =,I4,9H    KE1 =,I4,8H   KE2 =,I4)
  308 FORMAT(18H    KD1-2, KE1-2 =,4I4/
     *       18H   KDS1-2,KES1-2 =,4I4/18H   KDS3-4,KES3-4 =,4I4)
  309 FORMAT(55H ROWS OF TERM TABLE CONTAINING PARENTS ARE RESPECTIVELY,
     1 2(I6,I3))
C
C --- SET USEFUL CONSTANTS
C
      M1=ISIG-IRHO
      M2=ISIGP-IRHOP
      M19=IRHO-IRHOP
      M20=ISIG-ISIGP
      I2HSH=IHSH+IHSH-1
C
C     JSNDIR,JANGDI=0  IMPLIES APPROPRIATE J2,J3 ARRAYS FOR CALL OF
C     NJGRAF HAVE NOT BEEN SET
C
      JSNDIR=0
      JANGDI=0
      ISPIND=0
      ISPINE=0
      IANGD=0
      IANGE=0
      ICALL=0
C
C     SET THE FAIL PARAMETERS .FALSE. INITIALLY SO THAT NJGRAF CAN
C     BE CALLED.
C
      FAILSD=.FALSE.
      FAILSE=.FALSE.
      FAILAD=.FALSE.
      FAILAE=.FALSE.
C
C --- SET N,L VALUES OF INTERACTING SHELLS
C
      NRHO=NJ(IRHO)
      LRHO=LJ(IRHO)
      NSIG=NJ(ISIG)
      LSIG=LJ(ISIG)
      NRHOP=NJ(IRHOP)
      LRHOP=LJ(IRHOP)
      NSIGP=NJ(ISIGP)
      LSIGP=LJ(ISIGP)
*     IF(JRELOP(4).NE.0.OR.JRELOP(5).NE.0)
*     IREL=1
      IF(IBUG1.NE.0 .OR. IBUG2.GT.2)
     * WRITE(IWRITE,307) NRHO,LRHO,NSIG,LSIG,NRHOP,LRHOP,NSIGP,LSIGP
C
C --- EVALUATE MULTIPLICATIVE FACTORS
C
      IL=IDEL(IRHO,ISIG)
      IR=IDEL(IRHOP,ISIGP)
      Q=NOSH1(IRHO)*(NOSH1(ISIG)-IL)*NOSH2(IRHOP)*(NOSH2(ISIGP)-IR)
      IEXCHG=2-IL-IR
       XMULT1=SQRT(Q/((LRHO+LRHO+1)*(LSIG+LSIG+1)))*0.5
       MULDSO=0
       MULESO=0
       MULDSS=0
       MULDSP=0
       MULESS=0
       MULESP=0
      DO 13 J=1,IHSH
   13  NBAR(J)=NOSH1(J)-IDEL(J,IRHO)-IDEL(J,ISIG)
      IDELP=0
      IF(M1.EQ.0) GO TO 15
      DO 14 J=IRHO+1,ISIG
   14  IDELP=IDELP+NBAR(J)
   15 IF(M2.EQ.0) GO TO 18
      DO 16 J=IRHOP+1,ISIGP
   16  IDELP=IDELP+NBAR(J)
   18 K = 1-MOD(IDELP,2)*2
      XX = XMULT1
**---
C not XMULT1 = XMULT1 * K * SQRT(REAL(J1QN2(I2HSH,2)*J1QN2(I2HSH,3)))
      XMULT1 = XMULT1 * K * SQRT(REAL(J1QN1(I2HSH,2)*J1QN1(I2HSH,3)))
CCC                                     2L^+1          2S^S+1
C    *                *ALP4  !ALPHA**2/4 -- before '03Nov18/19
**---
      if(IBUG1.gt.1) then
       print"(' FANOFS: XX,XMULT1 =',F10.6,
     *  1P,E12.5E1,3I4/22X,24HI2HSH,J1QNi(-,2:3) =  '',2I4)", XX,
     *  XMULT1, i2hsh, J1QN1(I2HSH,2),J1QN1(I2HSH,3)
     *               , J1QN2(I2HSH,2),J1QN2(I2HSH,3)
      else
       if(IBUG2.gt.3) print"(' FANOFS',52X,'X,MLT=',2F7.3)",XX,XMULT1
      endif
C
C --- DETERMINE RANGES OF K FOR DIRECT AND EXCHANGE INTEGRALS
C     TRIANGULAR RELATIONS LIMIT (K+1) VALUES TO LIE BETWEEN KD1 AND KD2
C     FOR =DIRECT= INTEGRALS,  KE1 AND KE2 FOR =EXCHANGE= INTEGRALS
C
      K1=ABS(LSIG-LSIGP)
      K2=LSIG+LSIGP
      K3=ABS(LRHO-LRHOP)
      K4=LRHO+LRHOP
      KD1=MAX(K1,K3)+1
      KD2=MIN(K2,K4)+1
      K1=ABS(LRHOP-LSIG)
      K2=LRHOP+LSIG
      K3=ABS(LRHO-LSIGP)
      K4=LRHO+LSIGP
      KE1=MAX(K1,K3)+1
      KE2=MIN(K2,K4)+1
*     IF(JRELOP(5).NE.0)
C     CALL SKLIM(LET) -- in-line:
       LA=2
       CALL KLIM(LRHO,LRHOP,LSIG,LSIGP,LA,K1,K2)
       KDS1=K1+1
       KDS2=K2+1
       CALL KLIM(LSIG,LSIGP,LRHO,LRHOP,LA,K1,K2)
       KDS3=K1+1
       KDS4=K2+1
       CALL KLIM(LRHO,LSIGP,LSIG,LRHOP,LA,K1,K2)
       KES1=K1+1
       KES2=K2+1
       CALL KLIM(LSIG,LRHOP,LRHO,LSIGP,LA,K1,K2)
       KES3=K1+1
       KES4=K2+1
      IF(IBUG1.GT.1)  WRITE(IWRITE,308) KD1,KD2,KE1,KE2,
     *   KDS1,KDS2,KES1,KES2, KDS3,KDS4,KES3,KES4
C
C --- THE LOWEST VALUES  KD1  AND  KD2  ARE ALWAYS ALLOWED (UNLESS THEY
C     ARE GREATER THAN KD2 AND KE2 RESPECTIVELY). ALLOWED K-VALUES THEN
C     STEP BY 2 TO SATISFY THE EVEN CONDITION OF THE REDUCED MATRIX
C     ELEMENTS, WHICH ARE THEN CALCULATED AND STORED
      LET = .TRUE.
      IF(KD1.LE.KD2) THEN
       IF(KD2.GE.99) GO TO 999
       DO 23 K=KD1,KD2,2
   23   RMEDIR(K)=RME(LRHO,LRHOP,K-1)*RME(LSIG,LSIGP,K-1)
       LET = .FALSE.
      ENDIF
      IF(KE1.LE.KE2) THEN
       IF(KE2.GE.99) GO TO 999
       DO 24 K=KE1,KE2,2
   24   RMEEX(K)=RME(LRHO,LSIGP,K-1)*RME(LSIG,LRHOP,K-1)
       LET = .FALSE.
      ENDIF
      IF (LET.and.IBUG1.GT.1) WRITE(IWRITE,301)
      IF ( KDS1.GT.KDS2.AND.KDS3.GT.KDS4.AND.
     *     KES1.GT.KES2.AND.KES3.GT.KES4 )  RETURN != SKLIM-LET
cc    IF(LET) RETURN
C
*     IF(JRELOP(4).NE.0)
      CALL SOORED
*     IF(JRELOP(5).NE.0)
      CALL SSRED
C
C --- SET SENIORITY, S AND L VALUES OF SPECTATOR SHELLS FOR SR MUMDAD
C                                            AND SETJ1  VIA /INTERM/
      DO 26 J=1,IHSH
       IF(IRHO.EQ.J) GO TO 29
       IF(ISIG.EQ.J) GO TO 29
       DO 28 K=1,3
   28   J1BAR1(J,K)=J1QN1(J,K)
   29  IF(IRHOP.EQ.J) GO TO 26
       IF(ISIGP.EQ.J) GO TO 26
       DO 30 K=1,3
   30   J1BAR2(J,K)=J1QN2(J,K)
       IF(IRHO.EQ.J) GO TO 26
       IF(ISIG.EQ.J) GO TO 26
C      CHECK THAT SPECTATOR SHELLS HAVE SAME RESPECTIVE QUANTUM NUMBERS
       DO K=1,3
        IF(J1BAR1(J,K).NE.J1BAR2(J,K)) GO TO 27
       ENDDO
   26  CONTINUE
      GO TO 32
   27 IF(IBUG1.GT.1) WRITE(IWRITE,302)
      GO TO 527
C
C --- FROM WHICH ROWS OF NTAB DO WE FIND THE QUANTUM NUMBERS WITH BARS
C     OR TILDES
C
   32 NELCTS=NOSH1(ISIG)
      K2=NTAB1(NELCTS,LSIG+1)
      NELCTS=NOSH1(IRHO)
      IF(M1.EQ.0) NELCTS=NOSH1(ISIG)-1
      K1=NTAB1(NELCTS,LRHO+1)
      NELCTS=NOSH2(ISIGP)
      K4=NTAB1(NELCTS,LSIGP+1)
      NELCTS=NOSH2(IRHOP)
      IF(M2.EQ.0) NELCTS=NOSH2(ISIGP)-1
      K3=NTAB1(NELCTS,LRHOP+1)
      IF(IBUG1.GT.1)  WRITE(IWRITE,309) K1,K2,K3,K4
      KK1=ITAB(K1)
      KK2=ITAB(K2)
      KK3=ITAB(K3)
      KK4=ITAB(K4)
C
C
C === SUM OVER QUANTUM NUMBERS WITH BARS OR TILDES
C
      DO 151 JJ2=1,KK2
C
C --- TEST TO SEE WHICH PARENT TERMS ARE ALLOWABLE.  ONLY TEST THIS ON
C     L AND S  VALUES AT THIS STAGE, BY MEANS OF TRIANGULAR CONDITIONS
C     FOR TWICE THE QUANTUM NUMBERS -- RELYING ON INTEGER QUANTITIES!
C
      IN3=2*LSIG
      IJK2=3*(JJ2-1)+JTAB(K2)
      DO 140 K=2,3
       IN1=NTAB(IJK2+K)-1
       IN2=J1QN1(ISIG,K)-1
       IF(IN1.GT.IN2+IN3) GO TO 151
       IF(IN1.LT.ABS(IN2-IN3)) GO TO 151
  140  IN3=1
      DO 152 JJ1=1,KK1
      IN3=2*LRHO
      IJK1=3*(JJ1-1)+JTAB(K1)
  162 DO 145 K=2,3
      IN1=NTAB(IJK1+K)-1
      IF(M1.EQ.0) GO TO 142
      IN2=J1QN1(IRHO,K)-1
      GO TO 143
  142 IN2=NTAB(IJK2+K)-1
  143 IF(IN1.GT.IN2+IN3) GO TO 152
      IF(IN1.LT.ABS(IN2-IN3)) GO TO 152
  145 IN3=1
      DO 153 JJ4=1,KK4
      IN3=2*LSIGP
      IJK4=3*(JJ4-1)+JTAB(K4)
      DO 33 K=2,3
       IN1=NTAB(IJK4+K)-1
       IN2=J1QN2(ISIGP,K)-1
       IF(IN1.GT.IN2+IN3) GO TO 153
       IF(IN1.LT.ABS(IN2-IN3)) GO TO 153
   33  IN3=1
C
C
      DO 154 JJ3=1,KK3
      IN3=2*LRHOP
      IJK3=3*(JJ3-1)+JTAB(K3)
      DO 34 K=2,3
       IN1=NTAB(IJK3+K)-1
       IN2=NTAB(IJK4+K)-1
       IF(M2.NE.0) IN2=J1QN2(IRHOP,K)-1
       IF(IN1.GT.IN2+IN3) GO TO 154
       IF(IN1.LT.ABS(IN2-IN3)) GO TO 154
   34  IN3=1
      if(IBUG1.gt.1)  print "(/' FANOFS: DO JJ3 =',I3,
     * '  at JJ2,JJ1,JJ4 =',3I3)", JJ3, JJ2,JJ1,JJ4
C
C     SUMMATIONS NOW PERFORMED OVER ALLOWED QUANTUM NUMBERS
C     THE TILDES CORRESPOND TO  IRHO=ISIG  AND/OR IRHOP=ISIGP
C
C --- SET THE REMAINING QUANTUM NUMBERS WITH BARS OR TILDES
C
      DO 35 K=1,3
       J1BAR1(IRHO,K)=NTAB(IJK1+K)
       J1BAR2(IRHOP,K)=NTAB(IJK3+K)
       IF(M1.NE.0) THEN
        J1BAR1(ISIG,K)=NTAB(IJK2+K)
       ELSE
        J1TLD1(K)=NTAB(IJK2+K)
       ENDIF
       IF(M2.NE.0) THEN
        J1BAR2(ISIGP,K)=NTAB(IJK4+K)
       ELSE
        J1TLD2(K)=NTAB(IJK4+K)
       ENDIF
   35  CONTINUE
C
C --- IS POTENTIAL DIAG. IN BARRED QU. NOS. FOR INTERACTING SHELLS
C
      I5=0
      I=ISIG
      GO TO 50
   42 I=IRHO
      IF(M1.NE.0) GO TO 50
      I5=I5+1
   45 I=ISIGP
      GO TO 50
   46 I=IRHOP
      IF(M2.EQ.0) GO TO 48
   50 I5=I5+1
      DO K=1,3
       IF(J1BAR1(I,K).NE.J1BAR2(I,K)) GO TO 154
      ENDDO
      GO TO (42,45,46,48),I5
   48 PICFP=1.0
C
C --- EVALUATE FRACTIONAL PARENTAGE COEFFICIENTS
C
      I=1
      CALL MUMDAD (I,ISIG,IRHO,M1,CFPLHS)
      PICFP=PICFP*CFPLHS
      IF(ABS(PICFP).LT.EPS) GO TO 154
      I=2
      CALL MUMDAD(I,ISIGP,IRHOP,M2,CFPRHS)
      PICFP=PICFP*CFPRHS
      IF(ABS(PICFP).LT.EPS) GO TO 154
C
C === SET UP J1,J2,J3 AND EVALUATE RECOUPLING COEFFICIENTS
C
C --- FIRST OF ALL, THE SPIN COEFFICIENTS
C
      I=3
      CALL SETJ1(I,IRHO,ISIG,IRHOP,ISIGP,ISPIND,ISPINE,KK1,KK2,KK3,KK4)
      IF (ISPIND.NE.0) GO TO 90
         CALL J23SPN(IRHO,ISIG,IRHOP,ISIGP,JSNDIR)
         ISPIND=1
C
C --- DIRECT SPIN INTEGRAL
C
   90 SPINDT=0.0
      IF(KD1.GT.KD2) GO TO 78
      IF (FAILSD) GO TO 78
        IF (ISPIND.NE.2) THEN
          CALL NJGRAF(SPINDT,FAILSD)
          ISPIND=2
          IF(FAILSD) GO TO 78
          CALL KNJ(J6CSD,J7CSD,J8CSD,J9CSD,JWCSD,K6SD,K7SD,K8SD,K9SD,
     +             KWSD,JDELSD,LDELSD,SMVRSD,MPSD,
     +             J6PSD,J7PSD,J8PSD,J9PSD,JWRDSD,NLSMSD,NBJSD,NB6JSD,
     +             K6CPSD,K7CPSD,K8CPSD,K9CPSD,JSM4SD,JSM5SD,JSM6SD,
     +             IN6JSD)
        ELSE
          CALL GENSUM(J6CSD,J7CSD,J8CSD,J9CSD,JWCSD,K6SD,K7SD,K8SD,K9SD,
     +               KWSD,JDELSD,LDELSD,SMVRSD,MPSD,
     +               J6PSD,J7PSD,J8PSD,J9PSD,JWRDSD,NLSMSD,NBJSD,NB6JSD,
     +               K6CPSD,K7CPSD,K8CPSD,K9CPSD,JSM4SD,JSM5SD,JSM6SD,
     +               IN6JSD,SPINDT)
C
        ENDIF
C
C     IEXCHG IS ZERO WHENEVER  M1=0=M2 , IN WHICH CASE THE EXCHANGE
C     INTEGRAL HAS ZERO COEFFICIENT.  THERE IS THEN NO POINT IN
C     CALCULATING THIS INTEGRAL, AND SPINEX IS SET ZERO (AT LABEL 93)
C     AS A MARKER OF THIS SITUATION
C
   78 SPINEX=0.0
      IF(IEXCHG.EQ.0) GO TO 94
C
C --- MODIFY J2 AND J3 TO CALCULATE THE EXCHANGE SPIN INTEGRAL
C
      IF (ISPINE.EQ.0) THEN
         I=1
         CALL MODJ23(I)
         ISPINE=1
      ENDIF
C
C --- EXCHANGE SPIN INTEGRAL
C
      IF(KE1.GT.KE2) GO TO 94
      IF (FAILSE) GO TO 94
        IF (ISPINE.NE.2) THEN
          CALL NJGRAF(SPINEX,FAILSE)
          ISPINE=2
          IF (FAILSE) GO TO 94
          CALL KNJ(J6CSE,J7CSE,J8CSE,J9CSE,JWCSE,K6SE,K7SE,K8SE,K9SE,
     +             KWSE,JDELSE,LDELSE,SMVRSE,MPSE,
     +             J6PSE,J7PSE,J8PSE,J9PSE,JWRDSE,NLSMSE,NBJSE,NB6JSE,
     +             K6CPSE,K7CPSE,K8CPSE,K9CPSE,JSM4SE,JSM5SE,JSM6SE,
     +             IN6JSE)
        ELSE
          CALL GENSUM(J6CSE,J7CSE,J8CSE,J9CSE,JWCSE,K6SE,K7SE,K8SE,K9SE,
     +               KWSE,JDELSE,LDELSE,SMVRSE,MPSE,
     +               J6PSE,J7PSE,J8PSE,J9PSE,JWRDSE,NLSMSE,NBJSE,NB6JSE,
     +               K6CPSE,K7CPSE,K8CPSE,K9CPSE,JSM4SE,JSM5SE,JSM6SE,
     +               IN6JSE,SPINEX)
        ENDIF
   94 IF(IBUG1.GT.1)  WRITE(IWRITE,306) SPINDT,SPINEX
C
C
C --- THE ANGULAR RECOUPLING COEFFICIENTS
C     SET J1,J2,J3   (COMPARE SPIN INTEGRALS)
C
C     IF BOTH SPIN INTEGRALS ARE ZERO,  THERE IS NO PURPOSE IN
C     CALCULATING THE ANGULAR INTEGRALS
C
      IF(ABS(SPINDT).LT.EPS.AND.ABS(SPINEX).LT.EPS) GO TO 55
C
      I=2
      CALL SETJ1(I,IRHO,ISIG,IRHOP,ISIGP,IANGD,IANGE,KK1,KK2,KK3,KK4)
      IF (IANGD.EQ.0) CALL J23ANG(IRHO,ISIG,IRHOP,ISIGP,JANGDI)
C
C>    %%%%%%%%%%%%%%%%%%%%%%%
*>    if(ICALL.le.0) go to 55 ! eliminating 2/r_ij case '04Jun04
C>    %%%%%%%%%%%%%%%%%%%%%%%
C>
*> 55 IF(IREL.EQ.1)  CALL FANO2B(IRHO,ISIG,IRHOP,ISIGP,
*>   :     PICFP,ICALL,KK1,KK2,KK3,KK4)
c>
*>    SUBROUTINE FANO2B(IRHO,ISIG,IRHOP,ISIGP,PICFP,ICALL,KK1,KK2,...
*
   55 M21=M17+1
*
* --- SET UP J1, J2, AND J3 ARRAYS
*
*  FIRST OF ALL, THE  RECOUPLING COEFFICIENT
*
      IF (ICALL.EQ.0) THEN
         I = 2
         CALL SETJ1(I,IRHO,ISIG,IRHOP,ISIGP,  0,0,KK1,KK2,KK3,KK4)
         MLIMIT=M16
         NJ1S=M8+11
         NJ23S=M18+1
         FAILAD = .FALSE.
         FAILAE = .FALSE.
C
         FREE(M17) = .TRUE.
         FREE(M21) = .TRUE.
         FREE(NJ1S) = .TRUE.
		 J1(M17) = 1
		 J1(M21) = 1
		 J1(NJ1S) = 1
         J23REL=.FALSE.
         CALL J23SP1(J23REL)
         CALL NJGRAF(ANG,FAILAD)
          CALL KNJ(J6CAD,J7CAD,J8CAD,J9CAD,JWCAD,K6AD,K7AD,K8AD,K9AD,
     +             KWAD,JDELAD,LDELAD,SMVRAD,MPAD,
     +             J6PAD,J7PAD,J8PAD,J9PAD,JWRDAD,NLSMAD,NBJAD,NB6JAD,
     +             K6CPAD,K7CPAD,K8CPAD,K9CPAD,JSM4AD,JSM5AD,JSM6AD,
     +             IN6JAD)
         CALL J23SP1(J23REL)
         CALL MODREL
         CALL NJGRAF(ANG,FAILAE)
          CALL KNJ(J6CAE,J7CAE,J8CAE,J9CAE,JWCAE,K6AE,K7AE,K8AE,K9AE,
     +             KWAE,JDELAE,LDELAE,SMVRAE,MPAE,
     +             J6PAE,J7PAE,J8PAE,J9PAE,JWRDAE,NLSMAE,NBJAE,NB6JAE,
     +             K6CPAE,K7CPAE,K8CPAE,K9CPAE,JSM4AE,JSM5AE,JSM6AE,
     +             IN6JAE)
      ENDIF
*
      I=3
      CALL SETJ1(I,IRHO,ISIG,IRHOP,ISIGP,0,0,KK1,KK2,KK3,KK4)
      if(IBUG1.gt.1) print "(' FANOBP  at  I = 3: KK1-4,ANG,FAILAE =',
     *   4I4,F10.5,L4)", KK1,KK2,KK3,KK4, ANG,FAILAE ! ### .gt.2
      J1(M15)=2
      J1(M16)=2
      MLIMIT=M16
      NJ1S = M8+11
      NJ23S = M18+1
      IF (ICALL.EQ.0) THEN
         ICALL=1
             FAILSD = .FALSE.
             FAILSE = .FALSE.
             FREE(M17) = .TRUE.
             FREE(M21) = .TRUE.
             FREE(NJ1S) = .TRUE.
		     J1(M17) = 1
		     J1(M21) = 1
		     J1(NJ1S) = 1
             J23REL=.FALSE.
             CALL J23SP1(J23REL)
         CALL NJGRAF(SPNDT1,FAILSD)
         if(IBUG1.gt.1) print"(' FANO2B exit NJGRAF 1: J23REL,SPNDT1 =',
     *      L3,F10.5)", J23REL,spndt1
          CALL KNJ(J6CSD,J7CSD,J8CSD,J9CSD,JWCSD,K6SD,K7SD,K8SD,K9SD,
     +             KWSD,JDELSD,LDELSD,SMVRSD,MPSD,
     +             J6PSD,J7PSD,J8PSD,J9PSD,JWRDSD,NLSMSD,NBJSD,NB6JSD,
     +             K6CPSD,K7CPSD,K8CPSD,K9CPSD,JSM4SD,JSM5SD,JSM6SD,
     +             IN6JSD)
         CALL J23SP1(J23REL)
         CALL MODREL
         CALL NJGRAF(SPNDT1,FAILSE)
          CALL KNJ(J6CSE,J7CSE,J8CSE,J9CSE,JWCSE,K6SE,K7SE,K8SE,K9SE,
     +             KWSE,JDELSE,LDELSE,SMVRSE,MPSE,
     +             J6PSE,J7PSE,J8PSE,J9PSE,JWRDSE,NLSMSE,NBJSE,NB6JSE,
     +             K6CPSE,K7CPSE,K8CPSE,K9CPSE,JSM4SE,JSM5SE,JSM6SE,
     +             IN6JSE)
      ENDIF
*     IF(JRELOP(4).EQ.0) GO TO 100
*
* --- EVALUATION OF THE DIRECT SPIN RECOUPLING COEFFICIENTS
*
      SPNDT1= 0.0
      SPNDT2= 0.0
      IF (FAILSD) GO TO 3
      IF(KD1.gt.KD2) GO TO 3
      J1(M17)=3
      J1(M21)=1
      J1(NJ1S)=3
C ICI GENSUMSD
          CALL GENSUM(J6CSD,J7CSD,J8CSD,J9CSD,JWCSD,K6SD,K7SD,K8SD,K9SD,
     +               KWSD,JDELSD,LDELSD,SMVRSD,MPSD,
     +               J6PSD,J7PSD,J8PSD,J9PSD,JWRDSD,NLSMSD,NBJSD,NB6JSD,
     +               K6CPSD,K7CPSD,K8CPSD,K9CPSD,JSM4SD,JSM5SD,JSM6SD,
     +               IN6JSD,SPNDT1)
      J1(M17)=1
      J1(M21)=3
C ICI GENSUMSD
          CALL GENSUM(J6CSD,J7CSD,J8CSD,J9CSD,JWCSD,K6SD,K7SD,K8SD,K9SD,
     +               KWSD,JDELSD,LDELSD,SMVRSD,MPSD,
     +               J6PSD,J7PSD,J8PSD,J9PSD,JWRDSD,NLSMSD,NBJSD,NB6JSD,
     +               K6CPSD,K7CPSD,K8CPSD,K9CPSD,JSM4SD,JSM5SD,JSM6SD,
     +               IN6JSD,SPNDT2)
*
* --- MULTIPLY SPIN RECOUPLING COEFFICIENTS BY PRODUCT OF
*     FRACTIONAL PARENTAGE COEFFICIENTS
*
  3   SPIND0=(SPNDT1+SPNDT2+SPNDT2)*PICFP
      SPIND1=(SPNDT1+SPNDT1+SPNDT2)*PICFP
      if(IBUG1.gt.0) print"(9x,'KD1,KD2,SPNDTi,PICFP =',2I3,2X,3F10.5)",
     * KD1,KD2,SPNDT1,SPNDT2,PICFP
*
* --- EVALUATION OF THE EXCHANGE SPIN RECOUPLING COEFFICIENTS
*
*
* --- IF M1=0=M2 THE EXCHANGE SPIN RECOUPLING COEFFICIENT HAS ZERO
*     COEFFICIENT.  THERE IS THEN NO POINT IN CALCULATING THIS
*     COEFFICIENT AND SPNEX1 AND SPNEX2 ARE SET ZERO
*
      SPNEX1= 0.0
      SPNEX2= 0.0
      IF (FAILSE) GO TO 5
      IF ((M1.EQ.0.AND.M2.EQ.0).OR.KE1.GT.KE2) GO TO 5
*
* --- MODIFY J2 AND J3 ARRAYS TO CALCULATE THE EXCHANGE SPIN
*     RECOUPLING COEFFICIENT
*
      J1(M17)=3
      J1(M21)=1
      J1(NJ1S)=3
C ICI GENSUMSE
          CALL GENSUM(J6CSE,J7CSE,J8CSE,J9CSE,JWCSE,K6SE,K7SE,K8SE,K9SE,
     +               KWSE,JDELSE,LDELSE,SMVRSE,MPSE,
     +               J6PSE,J7PSE,J8PSE,J9PSE,JWRDSE,NLSMSE,NBJSE,NB6JSE,
     +               K6CPSE,K7CPSE,K8CPSE,K9CPSE,JSM4SE,JSM5SE,JSM6SE,
     +               IN6JSE,SPNEX1)
      J1(M17)=1
      J1(M21)=3
C ICI GENSUMSE
          CALL GENSUM(J6CSE,J7CSE,J8CSE,J9CSE,JWCSE,K6SE,K7SE,K8SE,K9SE,
     +               KWSE,JDELSE,LDELSE,SMVRSE,MPSE,
     +               J6PSE,J7PSE,J8PSE,J9PSE,JWRDSE,NLSMSE,NBJSE,NB6JSE,
     +               K6CPSE,K7CPSE,K8CPSE,K9CPSE,JSM4SE,JSM5SE,JSM6SE,
     +               IN6JSE,SPNEX2)
*
* --- MULTIPLY SPIN RECOUPLING COEFFICIENTS BY PRODUCT OF
*     FRACTIONAL PARENTAGE COEFFICIENTS
*
  5   if(IBUG1.gt.0) print "(9X,'M1,M2,KE1,KE2,FAILSE,SPNEX1,SPNEX2 =',
     *           4I3,L3,2F9.5)", M1,M2,KE1,KE2,FAILSE,SPNEX1,SPNEX2
      SPINE0=(SPNEX1+SPNEX2+SPNEX2)*PICFP
      SPINE1=(SPNEX2+SPNEX1+SPNEX1)*PICFP
*
* --- IF THE DIRECT AND EXCHANGE SPIN RECOUPLING COEFFICIENTS ARE
*     ZERO,  WE NEED NOT CALCULATE THE CORRESPONDING ORBITAL RECOUPLING
*     COEFFICIENTS
*
      IF(ABS(SPIND0).LT.EPS.AND.ABS(SPIND1).LT.EPS.AND.ABS(SPINE0).LT
     :     .EPS.AND.ABS(SPINE1).LT.EPS) GO TO 100
      I=2
      CALL SETJ1(I,IRHO,ISIG,IRHOP,ISIGP,2,2,KK1,KK2,KK3,KK4)
      NJ1S=M8+11
      NJ23S=M18+1
*
* --- EVALUATION OF THE DIRECT ORBITAL RECOUPLING COEFFICIENTS
*
*
* --- IF THE DIRECT SPIN RECOUPLING COEFFICIENT IS ZERO, WE NEED NOT
*     CALCULATE THE CORRESPONDING ORBITAL RECOUPLING COEFFICIENT
*
   17 IF(ABS(SPIND0).LT.EPS.AND.ABS(SPIND1).LT.EPS) GO TO 19
      IF (FAILAD) GO TO 19
      L=0
      DO 20 J=KD1,KD2,2
*
* --- CONSIDER ALL ALLOWED K-VALUES
*
      N=0
      L=L+1
      K=J-1
      JI=J+2
      DO 21 JJ=J,JI
      N=N+1
      I=JJ-2
      IF(I.LE.0.AND.K.EQ.0) GO TO 21
      J1(M17)=I+I+1
      J1(M21)=K+K+1
          CALL GENSUM(J6CAD,J7CAD,J8CAD,J9CAD,JWCAD,K6AD,K7AD,K8AD,K9AD,
     +               KWAD,JDELAD,LDELAD,SMVRAD,MPAD,
     +               J6PAD,J7PAD,J8PAD,J9PAD,JWRDAD,NLSMAD,NBJAD,NB6JAD,
     +               K6CPAD,K7CPAD,K8CPAD,K9CPAD,JSM4AD,JSM5AD,JSM6AD,
     +               IN6JAD,ANGDIR)
      DMULT0(L,N)=DMULT0(L,N)+ANGDIR*SPIND0*(1-MOD(I+K+1,2)*2)
      IF(I.EQ.K) GO TO 22
      J1(M17)=K+K+1
      J1(M21)=I+I+1
          CALL GENSUM(J6CAD,J7CAD,J8CAD,J9CAD,JWCAD,K6AD,K7AD,K8AD,K9AD,
     +               KWAD,JDELAD,LDELAD,SMVRAD,MPAD,
     +               J6PAD,J7PAD,J8PAD,J9PAD,JWRDAD,NLSMAD,NBJAD,NB6JAD,
     +               K6CPAD,K7CPAD,K8CPAD,K9CPAD,JSM4AD,JSM5AD,JSM6AD,
     +               IN6JAD,ANGDIR)
   22 DMULT1(L,N)=DMULT1(L,N)+ANGDIR*SPIND1
*
* --- MULDSO=1 WHEN A DIRECT INTEGRAL COEFFICIENT HAS BEEN CALCULATED
*     FOR USE, SEE RADWTS
*
      MULDSO=1
      if(IBUG1.gt.0) print"(9X,'J,JJ,I,K,M17,M21,ANGDIR  = ',6I4,2F9.5)"
     *  ,J,JJ,I,K,M17,M21,ANGDIR
   21 CONTINUE
   20 CONTINUE
*
* --- EVALUATION OF THE EXCHANGE ORBITAL RECOUPLING COEFFICIENTS
*
*
* --- IF THE EXCHANGE SPIN RECOUPLING COEFFICIENT IS ZERO, WE NEED NOT
*     CALCULATE THE CORRESPONDING ORBITAL RECOUPLING COEFFICIENT
*
   19 if(IBUG1.gt.0) print "(9X,'>>> SPINE0-E1,M1,M2,FAILAE =',
     *  2F9.4,2I3,L3)", SPINE0,SPINE1,M1,M2,FAILAE
      IF(ABS(SPINE0).LT.EPS.AND.ABS(SPINE1).LT.EPS) GO TO 99
      IF(M1.EQ.0.AND.M2.EQ.0) GO TO 99
*
* --- CONSIDER ALL ALLLOWED K-VALUES
*
      IF(KE1.GT.KE2 .OR. FAILAE) GO TO 99
      L=0
      DO J=KE1,KE2,2
      N=0
      L=L+1
      K=J-1
      JI=J+2
      DO 31 JJ=J,JI
      N=N+1
      I=JJ-2
      IF(I.LE.0.AND.K.EQ.0) GO TO 31
*
* --- MODIFY J2 AND J3 ARRAYS TO CALCULATE THE EXCHANGE ORBITAL
*     RECOUPLING COEFFICIENT
*
      J1(M17)=I+I+1
      J1(M21)=K+K+1
C ICI GENSUMAE
          CALL GENSUM(J6CAE,J7CAE,J8CAE,J9CAE,JWCAE,K6AE,K7AE,K8AE,K9AE,
     +               KWAE,JDELAE,LDELAE,SMVRAE,MPAE,
     +               J6PAE,J7PAE,J8PAE,J9PAE,JWRDAE,NLSMAE,NBJAE,NB6JAE,
     +               K6CPAE,K7CPAE,K8CPAE,K9CPAE,JSM4AE,JSM5AE,JSM6AE,
     +               IN6JAE,ANGEX)
      EMULT0(L,N)=EMULT0(L,N)-ANGEX*SPINE0*(1-MOD(I+K+1,2)*2)
      IF(I.EQ.K) GO TO 132
      J1(M17)=K+K+1
      J1(M21)=I+I+1
C ICI GENSUMAE
          CALL GENSUM(J6CAE,J7CAE,J8CAE,J9CAE,JWCAE,K6AE,K7AE,K8AE,K9AE,
     +               KWAE,JDELAE,LDELAE,SMVRAE,MPAE,
     +               J6PAE,J7PAE,J8PAE,J9PAE,JWRDAE,NLSMAE,NBJAE,NB6JAE,
     +               K6CPAE,K7CPAE,K8CPAE,K9CPAE,JSM4AE,JSM5AE,JSM6AE,
     +               IN6JAE,ANGEX)
  132 EMULT1(L,N)=EMULT1(L,N)-ANGEX*SPINE1
*
* --- MULESO=1 WHEN AN EXCHANGE INTEGRAL COEFFICIENT HAS BEEN CALCULATED
*     FOR USE,  SEE RADWTS
*
      MULESO=1
      if(IBUG1.gt.0)print"(9X,'J,JJ,K,ANGEX = ',3I4,2F9.5)",J,JJ,K,ANGEX
   31 CONTINUE
      ENDDO ! end J=KE1,KE2,2
C
*  99 IF(JRELOP(5).EQ.0) RETURN
   99 I=3
      CALL SETJ1(I,IRHO,ISIG,IRHOP,ISIGP,2,2,KK1,KK2,KK3,KK4)
      J1(M15)=2
      J1(M16)=2
      NJ1S=M8+11
      NJ23S=M18+1
*
* --- EVALUATION OF THE DIRECT  SPIN RECOUPLING COEFFICIENT
*
* 100 IF(JRELOP(5).EQ.0) RETURN
  100 SPINDT= 0.0
      IF(KDS1.GT.KDS2.AND.KDS3.GT.KDS4) GO TO 102
      IF (FAILSD) GO TO 102
      J1(M17)=3
      J1(M21)=3
      J1(NJ1S)=5
C ICI GENSUMSD
          CALL GENSUM(J6CSD,J7CSD,J8CSD,J9CSD,JWCSD,K6SD,K7SD,K8SD,K9SD,
     +               KWSD,JDELSD,LDELSD,SMVRSD,MPSD,
     +               J6PSD,J7PSD,J8PSD,J9PSD,JWRDSD,NLSMSD,NBJSD,NB6JSD,
     +               K6CPSD,K7CPSD,K8CPSD,K9CPSD,JSM4SD,JSM5SD,JSM6SD,
     +               IN6JSD,SPINDT)
  102 if(IBUG1.gt.0) print "(' FANO2B 102 RC: FAILSD, SPINDT*PICFP =',
     *     L4,2F10.5)", FAILSD,SPINDT,PICFP
*
* --- MULTIPLY SPIN RECOUPLING COEFFICIENT BY PRODUCT OF FRACTIONAL
*     PARENTAGE COEFFICIENTS
*
      SPINDT=SPINDT*PICFP
*
* --- EVALUATION OF THE EXCHANGE SPIN RECOUPLING COEFFICIENT
*
*
* --- IF M1=0=M2 THE EXCHANGE SPIN RECOUPLING COEFFICIENT HAS ZERO
*     COEFFICIENT.  THERE IS THEN NO POINT IN CALCULATING THIS
*     COEFFICIENT AND SPINEX IS SET ZERO
*
      SPINEX= 0.0
      IF(M1.EQ.0.AND.M2.EQ.0) GO TO 104
      IF(KES1.GT.KES2.AND.KES3.GT.KES4) GO TO 104
*
* --- MODIFY J2 AND J3 ARRAYS TO CALCULATE THE EXCHANGE SPIN RECOUPLING
*     COEFFICIENT
*
      IF (FAILSE) GO TO 104
      J1(M17)=3
      J1(M21)=3
      J1(NJ1S)=5
C ICI GENSUMSE
          CALL GENSUM(J6CSE,J7CSE,J8CSE,J9CSE,JWCSE,K6SE,K7SE,K8SE,K9SE,
     +               KWSE,JDELSE,LDELSE,SMVRSE,MPSE,
     +               J6PSE,J7PSE,J8PSE,J9PSE,JWRDSE,NLSMSE,NBJSE,NB6JSE,
     +               K6CPSE,K7CPSE,K8CPSE,K9CPSE,JSM4SE,JSM5SE,JSM6SE,
     +               IN6JSE,SPINEX)
  104 if(IBUG1.gt.0) print "(' FANO2B 104 RC: FAILSE, SPINEX*PICFP =',
     *     L4,2F10.5)", FAILSE,SPINEX,PICFP
*
* --- MULTIPLY SPIN RECOUPLING COEFFICIENT BY PRODUCTS OF FRACTIONAL
*     PARENTAGE COEFFICIENTS
*
      SPINEX=SPINEX*PICFP
*
* --- IF THE DIRCET AND EXCHANGE SPIN RECOUPLING COEFFICIENTS ARE
*     ZERO, WE NEED NOT CALCULATE THE CORRESPONDING ORBITAL RECOUPLING
*     COEFFICIENTS
*
      IF(ABS(SPINDT).LT.EPS.AND.ABS(SPINEX).LT.EPS) GO TO 154
      I=2
      CALL SETJ1(I,IRHO,ISIG,IRHOP,ISIGP,2,2,KK1,KK2,KK3,KK4)
      NJ1S=M8+11
      NJ23S=M18+1
*
* --- EVALUATION OF THE DIRECT RECOUPLING COEFFICIENTS
*
*
* --- IF THE DIRECT SPIN RECOUPLING COEFFICIENT IS ZERO, WE NEED NOT
*     CALCULATE THE CORRESPONDING ORBITAL RECOUPLING COEFFICIENT
*
      IF(ABS(SPINDT).LT.EPS) GO TO 59
      IF(KDS1.GT.KDS2 .OR. FAILAD) GO TO 39
*
* --- CONSIDER ALL ALLOWED K-VALUES
*
      L=0
      DO 36 J=KDS1,KDS2,2
       K=J-1
       L=L+1
       J1(M17)=K+K+5
       J1(M21)=K+K+1
C ICI  GENSUMAD
          CALL GENSUM(J6CAD,J7CAD,J8CAD,J9CAD,JWCAD,K6AD,K7AD,K8AD,K9AD,
     +               KWAD,JDELAD,LDELAD,SMVRAD,MPAD,
     +               J6PAD,J7PAD,J8PAD,J9PAD,JWRDAD,NLSMAD,NBJAD,NB6JAD,
     +               K6CPAD,K7CPAD,K8CPAD,K9CPAD,JSM4AD,JSM5AD,JSM6AD,
     +               IN6JAD,ANGDIR)
*
       MULDSS=1 ! WHEN A DIRECT RECOUPLING COEFFICIENT HAS BEEN
*       CALCULATED FOR USE, SEE READ2B
*
   36  DMULS0(L)=DMULS0(L)+ANGDIR*SPINDT
*
* --- CONSIDER ALL ALLOWED K-VALUES
*
   39 IF(KDS3.GT.KDS4) GO TO 59
      IF(M1.EQ.0.OR.M2.EQ.0 .OR. FAILAD) GO TO 59
      L=0
      DO J=KDS3,KDS4,2
      K=J-1
      L=L+1
      J1(M17)=K+K+1
      J1(M21)=K+K+5
C ICI GENSUMAD
          CALL GENSUM(J6CAD,J7CAD,J8CAD,J9CAD,JWCAD,K6AD,K7AD,K8AD,K9AD,
     +               KWAD,JDELAD,LDELAD,SMVRAD,MPAD,
     +               J6PAD,J7PAD,J8PAD,J9PAD,JWRDAD,NLSMAD,NBJAD,NB6JAD,
     +               K6CPAD,K7CPAD,K8CPAD,K9CPAD,JSM4AD,JSM5AD,JSM6AD,
     +               IN6JAD,ANGDIR)
      DMULS1(L)=DMULS1(L)+ANGDIR*SPINDT
*
      MULDSP=1 !WHEN A DIRECT RECOUPLING COEFFICIENT HAS BEEN CALCULATED
*
      ENDDO
*
* --- EVALUATION OF THE EXCHANGE ORBITAL RECOUPLING COEFFICIENTS
*
*
* --- IF THE EXCHANGE SPIN RECOUPLING COEFFICIENT IS ZERO, WE NEED NOT
*     CALCULATE THE CORRESPONDING ORBITAL RECOUPLING COEFFICIENT
*
   59 IF(ABS(SPINEX).LT.EPS) GO TO 154
c#### IF(ABS(SPINEX).ge.EPS) GO TO 154 ! ########### bypass test 04Jul31
      IF(KES1.GT.KES2.OR.M2.EQ.0 .OR. FAILAE) GO TO 60
*
* --- CONSIDER ALL ALLOWED K-VALUES
*
      L=0
      DO 51 J=KES1,KES2,2
       K=J-1
       L=L+1
*
* ---  MODIFY J2 AND J3 ARRAYS TO CALCULATE THE EXCHANGE ORBITAL
*      RECOUPLING COEFFICIENT
*
       J1(M17)=K+K+5
       J1(M21)=K+K+1
C ICI  GENSUMAE
          CALL GENSUM(J6CAE,J7CAE,J8CAE,J9CAE,JWCAE,K6AE,K7AE,K8AE,K9AE,
     +               KWAE,JDELAE,LDELAE,SMVRAE,MPAE,
     +               J6PAE,J7PAE,J8PAE,J9PAE,JWRDAE,NLSMAE,NBJAE,NB6JAE,
     +               K6CPAE,K7CPAE,K8CPAE,K9CPAE,JSM4AE,JSM5AE,JSM6AE,
     +               IN6JAE,ANGEX)
       EMULS0(L)=EMULS0(L)-ANGEX*SPINEX
*
* ---  MULESS=1 WHEN AN EXCHANGE RECOUPLING COEFFICIENT HAS BEEN
*      CALCULATED FOR USE, SEE RADWTS
*
   51  MULESS=1
      if(ibug1+IBUG2.gt.4) print "(' FANOFS past label 59: EMULSO =',
     *  5F9.4)",(EMULS0(J),J=1,L)
   60 IF(KES3.GT.KES4.OR.M1.EQ.0 .OR. FAILAE) GO TO 154
      L=0
*
* --- CONSIDER ALL ALLOWED K-VALUES
*
      DO 61 J=KES3,KES4,2
       K=J-1
       L=L+1
*
* ---  MODIFY J2 AND J3 ARRAYS TO CALCULATE THE EXCHANGE ORBITAL
*      RECOUPLING COEFFICIENT
*
       J1(M17)=K+K+1
       J1(M21)=K+K+5
C ICI  GENSUMAE
          CALL GENSUM(J6CAE,J7CAE,J8CAE,J9CAE,JWCAE,K6AE,K7AE,K8AE,K9AE,
     +               KWAE,JDELAE,LDELAE,SMVRAE,MPAE,
     +               J6PAE,J7PAE,J8PAE,J9PAE,JWRDAE,NLSMAE,NBJAE,NB6JAE,
     +               K6CPAE,K7CPAE,K8CPAE,K9CPAE,JSM4AE,JSM5AE,JSM6AE,
     +               IN6JAE,ANGEX)
*
       MULESP=1 ! WHEN AN EXCHANGE RECOUPLING COEFFICIENT HAS BEEN
*      CALCULATED FOR USE, SEE READ2B
*
   61  EMULS1(L)=EMULS1(L)-ANGEX*SPINEX
c >>>>>>>>>>>>
c >>>>>>>>>>>>
C
  154 CONTINUE
  153 CONTINUE
  152 CONTINUE
  151 CONTINUE
C
C
C
*     IF(JRELOP(4).NE.0)
  527 CALL SOOPAR(XMULT1)
*     IF(JRELOP(5).NE.0)
      CALL SSPAR(XMULT1)
      RETURN
C
C *** DEFINITION OF DIMENSION LIST
C
C     RMEDIR(K),K=KD1,KD2,2  -  DIRECT REDUCED MATRIX ELEMENT PRODUCT
C     RMEEX(K),K=KE1,KE2,2  -  EXCHANGE REDUCED MATRIX ELEMENT PRODUCT
C     KD1,KE1   ARE ALWAYS .GE. 1
C     KD2,KE2   ARE .LE.  1+2*MAX(L-VALUE)
C     NBAR(I), I=1,IHSH  -  NUMBER OF SPECTATOR ELECTRONS IN EACH SHELL
C      THE OTHER ARRAYS ARE DEFINED IN NJGRAF.
C
  999 PRINT *,' FANOFS: FIXED LENGTH 99 TOO SHORT IN former /XATION/.'
      STOP
C
      END
C***********************************************************************
C
C***********************************************************************
      SUBROUTINE SOOPAR(XMULT1)
C{-R8}IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER(LL53=  21*2+3,LL75=  21+2, NL1=LL75,NL2=NL1+2,NL3=NL1+5)
      COMMON/DEBUG/IBUG1,IBUG2,IBUG3,IBUG4,IBUG5,IBUG6,IBUG7,IBUG8,IBUG9
C     COMMON/DEMULT/DMULT0(12,3),DMULT1(12,3),EMULT0(12,3),EMULT1(12,3),
C    1     MULDSO,MULESO
C     COMMON/IMAGNT/CONST,CONSOO,CONSS,ICM1,ISPORB,ISOORB,ISPSPN,
C    1     ITENPR,IELST,IMASS,IDAR
C
C     COMMON/XATION/AMULT(25),BMULT(25),KD1,KD2,KE1,KE2,MULTD,MULTE
      COMMON/DEMULT/DMULT0(NL2,3),DMULT1(NL2,3),EMULT0(NL2,3),
     *       EMULT1(NL2,3),MULDSO,MULESO, KD1,KD2,KE1,KE2
      COMMON/DENKVK/D00N2(NL2),D00NK(NL2),D00VK(NL2),D11N2(NL2),
     1D11NK(NL2),D11VK(NL2),E01N2(NL2),E01NK(NL2),E01VK(NL2),
     2E10N2(NL2),E10NK(NL2),E10VK(NL2)
      COMMON/RMEPRD/RMEDIR(99),RMEEX(99)
      COMMON/MEDEFN/IHSH,NJ(LL75),LJ(LL75),NOSH1(LL75),NOSH2(LL75),
     1              J1QN1(LL53,3),J1QN2(LL53,3),IJFUL(LL75)
C
      COMMON/INFORM/IREAD,IWRITE,IPUNCH
      COMMON/MVALUE/M1,M2,M3,M4,M5,M6,M7,M8,M9,M10,M11,M12,M13,M14,M15,
     1     M16,M17,M18,M19,M20
C     COMMON/RMEPRD/RMEDIR(25),RMEEX(25)
C     COMMON/SOORME/RMEDR0(12,3),RMEDR1(12,3),RMEEX0(12,3),RMEEX1(12,3)
      COMMON/SOORME/RMEDR0(NL2,3),RMEDR1(NL2,3),RMEEX0(NL2,3),
     1RMEEX1(NL2,3)
C
C --- EVALUATE THE COEFFICIENTS
C
C --- DIRECT INTEGRAL
C
      XC=XMULT1*2.0
      IF(IBUG1.GT.0) WRITE(IWRITE,"(8X,12HSOOPAR: XC =,F12.6)") XC
      IF(MULDSO.LT.1) GO TO 2
      L=0
      DO 3 J=KD1,KD2,2
      L=L+1
      K=J-1
C
C --- INCLUDE MULTIPLICATIVE FACTORS COMMON TO ALL TERMS WITHIN THE
C     FOUR-FOLD SUMMATION
C
      FLOAT1=K+1
      FLOAT3=K+K+1
      IF(K.EQ.0) GO TO 50
      AN21=FLOAT3*SQRT(REAL(K+K-1))
   50 AN22=-FLOAT1*SQRT(FLOAT3)
      ANK2=K*SQRT(FLOAT3)
      ANK3=-FLOAT3*SQRT(REAL(K+K+3))
      AVK=SQRT(K*FLOAT1*FLOAT3)
      IF(K) 3,6,7
    7 D00N2(L)=(AN21*DMULT0(L,1)*RMEDR0(L,1)+AN22*DMULT0(L,2)*RMEDR0(L,2
     1))*XC
    6 D00NK(L)=(ANK2*DMULT0(L,2)*RMEDR0(L,2)+ANK3*DMULT0(L,3)*RMEDR0(L,3
     1))*XC
      IF(K.EQ.0) GO TO 5
      D00VK(L)=AVK*DMULT0(L,2)*RMEDIR(J)*XC
    5 IF(M1.EQ.0.OR.M2.EQ.0) GO TO 3
      IF(K) 3,8,9
    9 D11N2(L)=(AN21*DMULT1(L,1)*RMEDR1(L,1)+AN22*DMULT1(L,2)*RMEDR1(L,2
     1))*XC
    8 D11NK(L)=(ANK2*DMULT1(L,2)*RMEDR1(L,2)+ANK3*DMULT1(L,3)*RMEDR1(L,3
     1))*XC
      IF(K.EQ.0) GO TO 3
      D11VK(L)=AVK*DMULT1(L,2)*RMEDIR(J)*XC
    3 CONTINUE
C
C --- EXCHANGE INTEGRAL
C
    2 IF(MULESO.LT.1) GO TO 10
      L=0
      DO 12 J=KE1,KE2,2
      L=L+1
      K=J-1
      FLOAT1=K+1
      FLOAT3=K+K+1
      IF(K.EQ.0) GO TO 60
      AN21=FLOAT3*SQRT(REAL(K+K-1))
   60 AN22=-FLOAT1*SQRT(FLOAT3)
      ANK2=K*SQRT(FLOAT3)
      ANK3=-FLOAT3*SQRT(REAL(K+K+3))
      AVK=SQRT(K*FLOAT1*FLOAT3)
      IF(M2.EQ.0) GO TO 14
      IF(K) 12,15,16
   16 E01N2(L)=(AN21*EMULT0(L,1)*RMEEX0(L,1)+AN22*EMULT0(L,2)*RMEEX0(L,2
     1))*XC
   15 E01NK(L)=(ANK2*EMULT0(L,2)*RMEEX0(L,2)+ANK3*EMULT0(L,3)*RMEEX0(L,3
     1))*XC
      IF(K.EQ.0) GO TO 14
      E01VK(L)=AVK*EMULT0(L,2)*RMEEX(J)*XC
   14 IF(M1.EQ.0) GO TO 12
      IF(K) 12,17,18
   18 E10N2(L)=(AN21*EMULT1(L,1)*RMEEX1(L,1)+AN22*EMULT1(L,2)*RMEEX1(L,2
     1))*XC
   17 E10NK(L)=(ANK2*EMULT1(L,2)*RMEEX1(L,2)+ANK3*EMULT1(L,3)*RMEEX1(L,3
     1))*XC
      IF(K.EQ.0) GO TO 12
      E10VK(L)=AVK*EMULT1(L,2)*RMEEX(J)*XC
   12 CONTINUE
   10 RETURN
      END
C***********************************************************************
      SUBROUTINE SOORED
C{-R8}IMPLICIT REAL*8 (A-H,O-Z)
c     PARAMETER(NL1=&L75,NL2=NL1+2,NL3=NL1+5)
      PARAMETER(NL1=  21+2,NL2=NL1+2,NL3=NL1+5)
      COMMON/DEBUG/IBUG1,IBUG2,IBUG3,IBUG4,IBUG5,IBUG6,IBUG7,IBUG8,IBUG9
C     COMMON/XATION/AMULT(25),BMULT(25),KD1,KD2,KE1,KE2,MULTD,MULTE
      COMMON/DEMULT/DMULT0(NL2,3),DMULT1(NL2,3),EMULT0(NL2,3),
     *       EMULT1(NL2,3),MULDSO,MULESO, KD1,KD2,KE1,KE2
      COMMON/INFORM/IREAD,IWRITE,IPUNCH
      COMMON/MVALUE/M1,M2,M3,M4,M5,M6,M7,M8,M9,M10,M11,M12,M13,M14,M15,
     1     M16,M17,M18,M19,M20
      COMMON/NJLJ/NRHO,LRHO,NSIG,LSIG,NRHOP,LRHOP,NSIGP,LSIGP
C     COMMON/RMEPRD/RMEDIR(25),RMEEX(25)
      COMMON/RMEPRD/RMEDIR(99),RMEEX(99)
C     COMMON/SOORME/RMEDR0(12,3),RMEDR1(12,3),RMEEX0(12,3),RMEEX1(12,3)
      COMMON/SOORME/RMEDR0(NL2,3),RMEDR1(NL2,3),RMEEX0(NL2,3),
     1RMEEX1(NL2,3)
C
C --- DIRECT INTEGRALS
C
  301 FORMAT(4X,9HRMEDRN = ,6F11.5)
  302 FORMAT(4X,9HRMEEXN = ,6F11.5)
      IF(KD1.GT.KD2) GO TO 2
      L=0
      DO 3 J=KD1,KD2,2
      RMEP=RMEDIR(J)
      N=0
      L=L+1
      K=J-1
      J2=J+2
      DO 4 JJ=J,J2
      N=N+1
      I=JJ-2
      DMULT0(L,N)= 0.0
      DMULT1(L,N)= 0.0
      RMEDR0(L,N)= 0.0
      RMEDR1(L,N)= 0.0
      IF(LRHOP.EQ.0.AND.LSIGP.EQ.0) GO TO 4
      IF(I.LE.0.AND.K.EQ.0) GO TO 4
C
C --- EPSILON=ZERO, EPSILONP=ZERO,  /EPSILON-EPSILONP/=ZERO
C
      IF(LRHOP.EQ.0) GO TO 5
      CALL RMECLL(LRHO,LRHOP,K,I,RMECL)
      RMEDR0(L,N)=RMECL*RMEP
    5 IF(M1.EQ.0.OR.M2.EQ.0) GO TO 4
C
C --- EPSILON=ONE, EPSILONP=ONE ,  /EPSILON-EPSILONP/=ZERO
C
      IF(LSIGP.EQ.0) GO TO 4
      CALL RMECLL(LSIG,LSIGP,K,I,RMECL)
      RMEDR1(L,N)=RMECL*RMEP
    4 CONTINUE
      IF(IBUG1.GT.1) WRITE(IWRITE,301) (RMEDR0(L,N),N=1,3),
     1     (RMEDR1(L,N),N=1,3)
    3 CONTINUE
C
C --- EXCHANGE INTEGTALS
C
    2 IF(M1.EQ.0.AND.M2.EQ.0) RETURN
      L=0
      IF(KE1.GT.KE2) GO TO 7
      DO 8 J=KE1,KE2,2
      RMEP=RMEEX(J)
      N=0
      L=L+1
      K=J-1
      J2=J+2
      DO 9 JJ=J,J2
      N=N+1
      I=JJ-2
      EMULT0(L,N)=0.0
      EMULT1(L,N)=0.0
      RMEEX0(L,N)=0.0
      RMEEX1(L,N)=0.0
      IF(LRHOP.EQ.0.AND.LSIGP.EQ.0) GO TO 9
      IF(I.LE.0.AND.K.EQ.0) GO TO 9
      IF(M2.EQ.0.OR.LSIGP.EQ.0) GO TO 10
C
C --- EPSILON=ZERO, EPSILONP=ONE,  /EPSILON-EPSILONP/=ONE
C
      CALL RMECLL(LRHO,LSIGP,K,I,RMECL)
      RMEEX0(L,N)=RMECL*RMEP
   10 IF(M1.EQ.0.OR.LRHOP.EQ.0) GO TO 9
C
C --- EPSILON=ONE, EPSILONP=ZERO, /EPSILON-EPSILONP/=ONE
C
      CALL RMECLL(LSIG,LRHOP,K,I,RMECL)
      RMEEX1(L,N)=RMECL*RMEP
    9 CONTINUE
      IF(IBUG1.GT.1) WRITE(IWRITE,302) (RMEEX0(L,N),N=1,3),
     1     (RMEEX1(L,N),N=1,3)
    8 CONTINUE
    7 RETURN
      END
C***********************************************************************
      SUBROUTINE SSPAR(XMULT1)
C{-R8}IMPLICIT REAL*8 (A-H,O-Z)
c     PARAMETER(NL1=&L75,NL2=NL1+2,NL3=NL1+5,THREE=3.0)
      PARAMETER(NL1=  21+2,NL2=NL1+2,NL3=NL1+5,THREE=3.0,FIVE=5.0)
      COMMON/DEBUG/IBUG1,IBUG2,IBUG3,IBUG4,IBUG5,IBUG6,IBUG7,IBUG8,IBUG9
      COMMON/DESSNK/D00SNK(NL2),D11SNK(NL2),E01SNK(NL2),E10SNK(NL2),
     :     MULDSS,MULDSP,MULESS,MULESP
      COMMON/SSKLIM/KDS1,KDS2,KDS3,KDS4,KES1,KES2,KES3,KES4
      COMMON/SSMULT/DMULS0(NL2),DMULS1(NL2),EMULS0(NL2),EMULS1(NL2)
C     COMMON/SSRME/DRRME0(12),DRRME1(12),EXRME0(12),EXRME1(12)
      COMMON/SSRME/DRRME0(NL2),DRRME1(NL2),EXRME0(NL2),EXRME1(NL2)
C
C --- EVALUATE THE COEFFICIENTS
C
C
C --- DIRECT INTEGRAL
C
      XC=XMULT1*SQRT(1.0/FIVE)*THREE
      IF(IBUG1.GT.0) WRITE(6,"(9X,11HSSPAR: XC =,F12.6)") XC
      IF(MULDSS.LT.1) GO TO 2
      L=0
      DO 3 J=KDS1,KDS2,2
       L=L+1
       K=J-1
C ---  INCLUDE MULTIPLICATIVE FACTORS COMMON TO ALL TERMS WITHIN THE
C      FOUR-FOLD SUMMATION
       DKJ=XC*SQRT(REAL((K+K+1)*(K+K+2)*(K+K+3)*(K+K+4)*(K+K+5)))
    3  D00SNK(L)=DKJ*DMULS0(L)*DRRME0(L)
    2 IF(MULDSP.LE.0) GO TO 4
      L=0
      DO 6 J=KDS3,KDS4,2
       L=L+1
       K=J-1
       DKJ=XC*SQRT(REAL((K+K+1)*(K+K+2)*(K+K+3)*(K+K+4)*(K+K+5)))
    6  D11SNK(L)=DKJ*DMULS1(L)*DRRME1(L)
C
C --- EXCHANGE INTEGRAL
C
    4 IF(MULESS.LE.0) GO TO 7
      L=0
      DO 8 J=KES1,KES2,2
       L=L+1
       K=J-1
       DKJ=XC*SQRT(REAL((K+K+1)*(K+K+2)*(K+K+3)*(K+K+4)*(K+K+5)))
    8  E01SNK(L)=DKJ*EMULS0(L)*EXRME0(L)
    7 IF(MULESP.EQ.0) RETURN
      L=0
      DO 9 J=KES3,KES4,2
       L=L+1
       K=J-1
       DKJ=XC*SQRT(REAL((K+K+1)*(K+K+2)*(K+K+3)*(K+K+4)*(K+K+5)))
    9  E10SNK(L)=DKJ*EMULS1(L)*EXRME1(L)
      RETURN
      END
C***********************************************************************
      SUBROUTINE SSRED
C{-R8}IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER(NL1=  21+2,NL2=NL1+2,NL3=NL1+5)
      COMMON/DEBUG/IBUG1,IBUG2,IBUG3,IBUG4,IBUG5,IBUG6,IBUG7,IBUG8,IBUG9
      COMMON/MVALUE/M1,M2,M3,M4,M5,M6,M7,M8,M9,M10,M11,M12,M13,M14,M15,
     1     M16,M17,M18,M19,M20
      COMMON/NJLJ/NRHO,LRHO,NSIG,LSIG,NRHOP,LRHOP,NSIGP,LSIGP
      COMMON/SSKLIM/KDS1,KDS2,KDS3,KDS4,KES1,KES2,KES3,KES4
C
      COMMON/SSMULT/DMULS0(NL2),DMULS1(NL2),EMULS0(NL2),EMULS1(NL2)
      COMMON/SSRME/DRRME0(NL2),DRRME1(NL2),EXRME0(NL2),EXRME1(NL2)
      IF(KDS1.GT.KDS2) GO TO 1
      L=0
      DO 5 J=KDS1,KDS2,2
       L=L+1
       K=J-1
       K2=K+2
       DMULS0(L)= 0.0
    5  DRRME0(L)=RME(LRHO,LRHOP,K2)*RME(LSIG,LSIGP,K)
    1 IF(KDS3.GT.KDS4) GO TO 2
      IF(M1.EQ.0.OR.M2.EQ.0) GO TO 2
      L=0
      DO 6 J=KDS3,KDS4,2
       L=L+1
       K=J-1
       K2=K+2
       DMULS1(L)= 0.0
    6  DRRME1(L)=RME(LSIG,LSIGP,K2)*RME(LRHO,LRHOP,K)
    2 IF(KES1.GT.KES2.OR.M2.EQ.0) GO TO 3
      L=0
      DO 7 J=KES1,KES2,2
       L=L+1
       K=J-1
       K2=K+2
       EMULS0(L)= 0.0
    7  EXRME0(L)=RME(LRHO,LSIGP,K2)*RME(LRHOP,LSIG,K)
    3 IF(KES3.GT.KES4.OR.M1.EQ.0) RETURN
      L=0
      DO 8 J=KES3,KES4,2
       L=L+1
       K=J-1
       K2=K+2
       EMULS1(L)=0.0
    8  EXRME1(L)=RME(LSIG,LRHOP,K2)*RME(LRHO,LSIGP,K)
      RETURN
      END
C     **********************************************************
C
C     **********************************************************
      SUBROUTINE MATANS(IMATX)
C
C     -- SETS UP THE CALL TO MKWTS FOR THE BOUND-BOUND
C     (IMATX=1), BOUND-CONTINUUM (IMATX=2), AND CONTINUUM-CONTINUUM
C     DIAGONAL AND OFF-DIAGONAL (IMATX=3) MATRIX ELEMENTS.
C
C     IMPLICIT REAL*8 (R)
      PARAMETER (LL53=  21*2+3, LL75=  21+2) ! cor Oct24: +20->+2
      COMMON /ANSWER/ SOO(  60,  60),SS(  60,  60)
      COMMON/BASICS/LRANG2,NRANG2
      COMMON/C2BODY/ K2BDY,LRANG3, LRGL,NSPN ,LLRGL,NNSPN
      COMMON/MEDEFN/IHSH,NJ(LL75),LJ(LL75),NOSH1(LL75),NOSH2(LL75),
     *              J1QN1(LL53,3),J1QN2(LL53,3),IJFUL(LL75)
      COMMON/MKMATX/KSAME,ILIMIT,JLIMIT
C
      IF(IMATX.LT.0) GO TO 8
      DO 4 J=1,NRANG2
       DO 4 I=1,NRANG2
        SOO(I,J) = 0.
    4    SS(I,J) = 0.
      IF(IMATX.EQ.0) RETURN ! for testing purposes '04Jun17/Jul13
    8 IF(ABS(IMATX)-2) 1,2,3
C
C      BOUND-BOUND
C
    1 KSAME=1
      ILIMIT=1
      JLIMIT=1
      CALL MKWTS
      RETURN
C
C      BOUND-CONTINUUM
C
    2 KSAME=2
      ILIMIT= NRANG2 !! 1
      JLIMIT=1 !! NRANG2 !NDCT(2)
      CALL MKWTS
      RETURN
C
C      CONTINUUM-CONTINUUM. OFF DIAGONAL I.E. KI.NE.KJ
C     SPNORB: (LRGL,NSPN) ON THE L.H.S. AND (LLRGL,NNSPN) ON THE R.H.S.
C
    3 ILIMIT=NRANG2 !NDCT(1)
      JLIMIT=NRANG2 !NDCT(2)
      IF(LJ(IHSH).EQ.LJ(IHSH-1)) GO TO 5
C
C     CONTINUUM-CONTINUUM. DIAGONAL I.E. KI= KJ IF L.NE.LP (I.E.LJ(IHSH
C      -1).NE.LJ(IHSH)). CALL MATRX AGAIN WITH SAME COUPLING SCHEMES
C
      KSAME=0
      CALL MKWTS
      RETURN
C
C     CONTINUUM-CONTINUUM. DIAGONAL I.E. KI=KJ BUT NOW L=LP. HENCE
C      CONTRACTION OF COUPLING SCHEMES SINCE THE TWO CONTINUUM ORBITALS
C      ARE IDENTICAL
C
    5 KSAME=4
      CALL MKWTS
      KSAME=3
c     ILIMIT = NRANG2 ! ??? ILIMIT=MIN(NDCT(1),NDCT(2))
c     JLIMIT=ILIMIT
      I10=IHSH
      IHSH=IHSH-1
      I11=2*IHSH-2
      I12=I11+1
      NOSH2(IHSH)=1
      J1QN2(IHSH,1)=1
      J1QN2(IHSH,2)=2*LJ(IHSH)+1
      J1QN2(IHSH,3)=2
      IF(I11.LT.I10) I11=I10
      DO 7 K=1,3
       DO 6 J=I10,I11
        J1QN1(J,K)=J1QN1(J+1,K)
    6   J1QN2(J,K)=J1QN2(J+1,K)
       J1QN1(I12,K)= J1QN1(I12+2,K) ! restored '06Jul11: KSAME=3 ok
    7  J1QN2(I12,K)= J1QN2(I12+2,K) ! 0, 2*LLGRL+1, NNSPN+1 on LHS?
C
C      NOW CALL MKWTS WITH CONTRACTED COUPLING SCHEMES
C
      CALL MKWTS
      RETURN
      END
C***********************************************************************
      SUBROUTINE MKWTS
C     -- CONTROLLING 2-BODY FS SET-UP
C
      PARAMETER (LL53=  21*2+3, LL75=  21+2)
      COMMON/DEBUG/IBUG1,IBUG2,IBUG3,IBUG4,IBUG5,IBUG6,IBUG7,IBUG8,IBUG9
      COMMON/INFORM/IREAD,IWRITE,IPUNCH
      COMMON/INSTO7/NBUG,INK1,INK4,NCNT,NMIN(0:   5)
      COMMON/MEDEFN/IHSH,NJ(LL75),LJ(LL75),NOSH1(LL75),NOSH2(LL75),
     *              J1QN1(LL53,3),J1QN2(LL53,3),IJFUL(LL75)
C
C     THE =INTERACTING= SHELLS ARE DESIGNATED  IRHO,ISIG,IRHOP,ISIGP.
C     THE FIRST TWO REFER TO THE L.H.S. OF (PSI/V/PSIP),
C     WHILE  THE SECOND TWO REFER TO THE R.H.S.  FOR DIAGONAL AND
C     CERTAIN OFF-DIAGONAL MATRIX ELEMENTS; THESE MAY NOT BE UNIQUE,
C     AND EACH  POSSIBILITY MUST BE CONSIDERED: IN TURN THE CONDITION
C     =IRHO.LE.ISIG , IRHOP.LE.ISIGP=  MUST BE SATISFIED
C
   61 FORMAT(/9X,5HRHO =,I3,4X,5HSIG =,I3,5X,6HRHOP =,I3,3X,6HSIGP =,I3)
C
C
C === DETERMINE THE INTERACTING SHELLS AS FAR AS POSSIBLE BY
C     CONSIDERING THE DIFFERENCES BETWEEN PSI AND PSIP
C
      IF(IBUG2.LE.2) go to 2
      print "(' MKWTS: NOSH1 =',15I3)", (nosh1(i),i=1,ihsh) !**
      print "('        NOSH2 =',15I3)", (nosh2(i),i=1,ihsh) !**
    2 IX=0
      IRHO=0
      ISIG=0
      IRHOP=0
      ISIGP=0
      DO 4 J=1,IHSH
      N=NOSH1(J)-NOSH2(J)
      IF(ABS(N).GT.2) GO TO 1
      IF(N) 7,4,6
    6 IF(N.GT.1) GO TO 9
C
C --- TO SATISFY =IRHO.LE.ISIG=  IRHO IS SET FIRST,  ETC.
C
      IF(IRHO.GT.0) GO TO 11
      IRHO = J
      GO TO 17
   11 ISIG=J
      GO TO 17
    9 IRHO=J
      GO TO 15
    7 IF(N+1) 13,14,4
   14 IF(IRHOP.GT.0) GO TO 16
      IRHOP = J
      GO TO 17
   16 ISIGP=J
   17 IX=IX+1
      GO TO 4
   13 IRHOP=J
   15 IX=IX+2
    4 CONTINUE
c
      if(IBUG1+IBUG2.gt.3) print "(/' MKWTS: IHSH,IX0 ='2I4)", IHSH,IX
C
C     IX MEASURES THE TOTAL NUMBER OF ELECTRONS IN EITHER CONFIGURATION
C     WHICH DO NOT OCCUR IN THE OTHER.  THEN IF  IX  IS GREATER THAN 4,
C     ORTHOGONALITY OF THE ORBITALS PREVENTS A NON-ZERO MATRIX ELEMENT.
C     IF  IX  IS LESS THAN 4, THEN WE DIVIDE IX BY 2 AND NOW IX MEASURES
C     THE NUMBER OF ELECTRONS WHICH HAVE BEEN CHANGED IN GOING FROM PSI
C     TO PSIP.  IF NOW IX=0, WE HAVE A DIAGONAL MATRIX ELEMENT.  RHO AND
C     SIG MAY TAKE ON ANY VALUES LESS THAN IHSH.  IF IX=1, ONE INTER-
C     ACTING SHELL ON EACH SIDE IS FIXED, WHILE THE OTHER MAY VARY.  IF
C     IX=2, ALL INTERACTING SHELLS ARE DETERMINED
C
      IF(IX.GT.4) GO TO 1
      IX=IX/2
      IF(IX-1) 19,20,21
C
C
C === UNIQUE SPECIFICATION OF INTERACTING SHELLS
C     -- FANOFS COMPUTES ANGULAR CORFFICIENTS AND ARGUMENTS OF MAGNETIC
C        INTEGRALS, READ2 ACCUMULATES THE PRODUCTS IN SOO(,) AND SS(,)
C
   21 IF(ISIG.EQ.0) ISIG=IRHO
      IF(ISIGP.EQ.0) ISIGP = IRHOP
      IF(IBUG1+IBUG2.GT.3) WRITE(IWRITE,61) IRHO,ISIG,IRHOP,ISIGP
      IF  (LJ(IRHO)+LJ(ISIG)+LJ(IRHOP)+LJ(ISIGP)-
     * MAX(LJ(IRHO),LJ(ISIG),LJ(IRHOP),LJ(ISIGP)).EQ.0) GO TO 1
      IF (NJ(IRHO).LT.NMIN(LJ(IRHO)) .OR.
     *    NJ(IRHOP).LT.NMIN(LJ(IRHOP)) .OR.
     *    NJ(ISIG).LT.NMIN(LJ(ISIG)) .OR.
     *    NJ(ISIGP).LT.NMIN(LJ(ISIGP))) GO TO 1
      CALL REDUCE(IRHO,ISIG,IRHOP,ISIGP,LESSEN)
      CALL SETM
      CALL FANOFS(IRHO,ISIG,IRHOP,ISIGP)
      IF(LESSEN.NE.0) CALL MEREST(IRHO,ISIG,IRHOP,ISIGP)
      CALL READ2B
      GO TO 1
C
C
C === ONE INTERACTING SHELL SPECIFIED ON EACH SIDE: SUMMATION OVER OTHER
C
   20 IRSTO=IRHO
      IRPSTO=IRHOP
      DO 125 K1=1,IHSH
       IF(NOSH1(K1).EQ.0) GO TO 125
       ISIG=K1
       IF(NOSH2(K1).EQ.0) GO TO 125
       ISIGP = K1
       IRHO=IRSTO
       IRHOP=IRPSTO
C
C      ORTHOGONALITY OF THE ORBITALS REQUIRES THAT THE VARYING INTER-
C       ACTING SHELL BE THE SAME ON BOTH SIDES OF THE MATRIX ELEMENT
C
C ---  IRHO.LE.ISIG,   IRHOP.LE.ISIGP
C
       IF(IRHO-ISIG) 27,26,25
   25  ISTO=IRHO
       IRHO=ISIG
       ISIG = ISTO
       GO TO 27
   26  IF(NOSH1(ISIG).LE.1) GO TO 125
   27  IF(IRHOP-ISIGP) 30,29,28
   28  ISTO=IRHOP
       IRHOP = ISIGP
       ISIGP = ISTO
       GO TO 30
   29  IF(NOSH2(ISIGP).LE.1) GO TO 125
   30  IF(IBUG1+IBUG2.GT.3) WRITE(IWRITE,61) IRHO,ISIG,IRHOP,ISIGP
       IF  (LJ(IRHO)+LJ(ISIG)+LJ(IRHOP)+LJ(ISIGP)-
     *  MAX(LJ(IRHO),LJ(ISIG),LJ(IRHOP),LJ(ISIGP)).EQ.0) GO TO 125
       IF (NJ(IRHO).LT.NMIN(LJ(IRHO)) .OR.
     *     NJ(IRHOP).LT.NMIN(LJ(IRHOP)) .OR.
     *     NJ(ISIG).LT.NMIN(LJ(ISIG)) .OR.
     *     NJ(ISIGP).LT.NMIN(LJ(ISIGP))) GO TO 125
       CALL REDUCE(IRHO,ISIG,IRHOP,ISIGP,LESSEN)
       CALL SETM
       CALL FANOFS(IRHO,ISIG,IRHOP,ISIGP)
       IF(LESSEN.NE.0) CALL MEREST(IRHO,ISIG,IRHOP,ISIGP)
       CALL READ2B
  125  CONTINUE
      GO TO 1
C
C
C === NO INTERACTING SHELLS SPECIFIED
C      SUMMATION OVER ALL POSSIBLE COMBINATIONS
C     IN THIS CASE ORTHOGONALITY OF ORBITALS PRECLUDES ALL CASES
C      EXCEPT  IRHO=IRHOP   AND   ISIG=ISIGP
C
   19 DO 32 K1=1,IHSH
       IF(NOSH1(K1).EQ.0) GO TO 32
       ISIG = K1
       ISIGP = K1
       DO 33 K2=1,K1
        IF(NOSH1(K2).EQ.0) GO TO 33
        IRHO = K2
        IRHOP = K2
        IF(K2.NE.K1) GO TO 50
        IF(NOSH1(K1).LE.1) GO TO 33
   50   IF(IBUG1+IBUG2.GT.3) WRITE(IWRITE,61) K2,K1,K2,K1
C ???   IF(LJ(K2).GT.4) GO TO 74 ! ???
C ???   IF(LJ(K1).GT.4) GO TO 74 ! ???
CX      IF(ICHOP(K1).EQ.1.OR.ICHOP(K2).EQ.1)GOTO 33;/REMOVE/ICHOP(LL75):
        IF (LJ(K2).EQ.0 .AND. LJ(K1).EQ.0) GO TO 33
        IF (NJ(K2).LT.NMIN(LJ(K2)) .OR. NJ(K1).LT.NMIN(LJ(K1))) GO TO 33
   74   CALL REDUCE(IRHO,ISIG,IRHOP,ISIGP,LESSEN)
        CALL SETM
        CALL FANOFS(IRHO,ISIG,IRHOP,ISIGP)
        IF(LESSEN.NE.0) CALL MEREST(IRHO,ISIG,IRHOP,ISIGP)
        CALL READ2B
   33   CONTINUE
   32  CONTINUE
C
    1 RETURN
      END
C***********************************************************************
      SUBROUTINE RMECLL(L,LP,K,K1,RMECL)
C
C     -- EVALUATES THE REDUCED MATRIX ELEMENT  (L//(C(K1)*L)(K)//LP)
C
C{-R8}IMPLICIT REAL*8 (R,S)
      RMECL=0.0
      IF(LP.EQ.0) RETURN
      L2=L+L
      LP2=LP+LP
      K12=K1+K1
      IL2=2
      K2=K+K
      CALL DRACAH(LP2,LP2,K12,K2,IL2,L2,RACAH)
      RMECL = SQRT(REAL(K12+1)*LP*(LP+1)*(LP2+1))*
     C        (1-MOD(ABS(LP+K-L),2)*2) * RACAH
      RETURN
      END
C***********************************************************************
      SUBROUTINE MODREL
      LOGICAL FREE
      PARAMETER(LL76=3*  21+19,LL74=  21+4) ! '03Oct25: +4, not +2
      COMMON/COUPLE/NJ1S,NJ23S,J1(LL76),J2(LL74,3),J3(LL74,3),FREE(LL76)
      COMMON/DEBUG/IBUG1,IBUG2,IBUG3,IBUG4,IBUG5,IBUG6,IBUG7,IBUG8,IBUG9
      COMMON/INFORM/IREAD,IWRITE,IPUNCH
      COMMON/MVALUE/M1,M2,M3,M4,M5,M6,M7,M8,M9,M10,M11,M12,M13,M14,M15,
     1     M16,M17,M18,M19,M20
      J2(1,1)=M16
      J2(2,1)=M15
      IF(IBUG1.GT.1.AND.IBUG3.NE.1)
     *  WRITE(IWRITE,"(3H J2,18X,2HJ3,15X,'in MODREL:'/(3I5,I10,2I5))")
     *    ((J2(J,K),K=1,3),(J3(J,K),K=1,3),J=1,M18)
      RETURN
      END
C***********************************************************************
      SUBROUTINE J23SP1(J23REL)
      PARAMETER(LL76=3*  21+19,LL74=  21+4) ! '03Oct25: +4, not +2
      LOGICAL J23REL, FREE
      DIMENSION J2REL(69),J3REL(69)
      COMMON/COUPLE/NJ1S,NJ23S,J1(LL76),J2(LL74,3),J3(LL74,3),FREE(LL76)
      COMMON/DEBUG/IBUG1,IBUG2,IBUG3,IBUG4,IBUG5,IBUG6,IBUG7,IBUG8,IBUG9
      COMMON/HOLD/ J2SPIN(LL76),J3SPIN(LL76),J2ANG(LL76),J3ANG(LL76)
      COMMON/INFORM/IREAD,IWRITE,IPUNCH
      COMMON/MVALUE/M1,M2,M3,M4,M5,M6,M7,M8,M9,M10,M11,M12,M13,M14,M15,
     1     M16,M17,M18,M19,M20
  304 FORMAT(3H J2,18X,2HJ3,15X,'in J23SP1:'/(3I5,I10,2I5))
C
      SAVE J2REL,J3REL
C
      IF(J23REL) GO TO 10
C
C --- ROWS 1 TO M4 OF SPIN J3(J23SPN) ARE THE SAME AS ROWS 1 TO M4 OF
C     SPIN J3(J23SP1); AND J2(") SAME AS ROWS 3 TO M4+2 OF SPIN J2(..P1)
C
      I1=0
      DO 1 J=1,M4
       DO 2 K=1,3
        I1=I1+1
        J3(J,K)=J3SPIN(I1)
    2   J2(J+2,K)=J2SPIN(I1)
    1  CONTINUE
C
C --- SET FIRST TWO ROWS OF J2, CORRESPONDING TO THE COUPLING OF
C     INTERACTING ELECTRONS WITHIN THEIR SHELLS
C
      J2(1,1)=M15
      J2(1,2)=M17
      J2(1,3)=M13
      J2(2,1)=M16
      J2(2,2)=M17+1
      J2(2,3)=M14
C
C --- RESET THE FOLLOWING ELEMENTS IN THE J3 ARRAY
C
      J3(1,2)=M15
      J3(2,2)=M16
      J3(M4,3)=M8
C
C --- SET ROWS M5 AND M5+1 OF J3
C
      J3(M5,1)=M17
      J3(M5,2)=M17+1
      J3(M5,3)=NJ1S
      J3(M18,1)=M8
      J3(M18,2)=NJ1S
      J3(M18,3)=M6
C
C --- STORE THE ARRAYS J2 AND J3 FOR FUTURE RESETTING OF J2 AND J3
C
      I1=0
      DO 21 J=1,M18
       DO 22 K=1,3
        I1=I1+1
        J2REL(I1)=J2(J,K)
   22   J3REL(I1)=J3(J,K)
   21  CONTINUE
      J23REL=.TRUE.
   26 IF(IBUG1.LT.2.OR.IBUG3.EQ.1) RETURN
      WRITE(IWRITE,304) ((J2(J,K),K=1,3),(J3(J,K),K=1,3),J=1,M18)
      RETURN
C
C --- RESET THE ARRAYS J2 AND J3 FROM STORE
C
   10 I1=0
      DO 24 J=1,M18
       DO 25 K=1,3
        I1=I1+1
        J2(J,K)=J2REL(I1)
   25   J3(J,K)=J3REL(I1)
   24  CONTINUE
      GO TO 26
      END
C***********************************************************************
      SUBROUTINE KLIM(L1,L2,L3,L4,L,KA,KB)
      K1=ABS(L1-L2)-L
      K2=L1+L2-L
      K3=ABS(L3-L4)
      K4=L3+L4
      KA=MAX(K1,K3)
      KB=MIN(K2,K4)
      RETURN
      END
C***********************************************************************
      SUBROUTINE READ2B
C     -- COMPUTES CUMULATIVE SOO(,) AND SS(,) FOR SPINxx
C{-R8}IMPLICIT REAL*8 (A-H,O-Z)
      character*6 type,mtyp
      LOGICAL IBUG0
      PARAMETER (LL75=  21+2, EPS=1.E-5, NL1=LL75,NL2=NL1+2,NL3=NL1+5)
      COMMON /ANSWER/ SOO(  60,  60),SS(  60,  60)
      COMMON/DEBUG/IBUG1,IBUG2,IBUG3,IBUG4,IBUG5,IBUG6,IBUG7,IBUG8,IBUG9
      COMMON/DEMULT/DMULT0(NL2,3),DMULT1(NL2,3),EMULT0(NL2,3),
     *       EMULT1(NL2,3),MULDSO,MULESO, KD1,KD2,KE1,KE2
      COMMON/DESSNK/D00SNK(NL2),D11SNK(NL2),E01SNK(NL2),E10SNK(NL2),
     :     MULDSS,MULDSP,MULESS,MULESP
      COMMON/DENKVK/D00N2(NL2),D00NK(NL2),D00VK(NL2),D11N2(NL2),
     1   D11NK(NL2),D11VK(NL2),E01N2(NL2),E01NK(NL2),E01VK(NL2),
     2   E10N2(NL2),E10NK(NL2),E10VK(NL2)
      COMMON/MVALUE/M1,M2,M3,M4,M5,M6,M7,M8,M9,M10,M11,M12,M13,M14,M15,
     1     M16,M17,M18,M19,M20
     *    /MKMATX/KSAME,ILIMIT,JLIMIT
     */MKMATS/RNKMAT(  60,  60,4),RNKMA2(  60,  60,4),VKMAT(  60,  60,4)
      COMMON/SSKLIM/KDS1,KDS2,KDS3,KDS4,KES1,KES2,KES3,KES4
C
C      DEFINITION OF KSAME
C      KSAME = 0  C-C ELEMENTS OF MKMAT WITH L.NE.LP
C      KSAME = 1  BOUND-BOUND MK INTEGRALS
C      KSAME = 2  BOUND-CONTINUUM; CONTINUUM ORBITAL IN ANY POSITION!
C      KSAME = 3  CONTINUUM-CONTINUUM DIAGONAL ELEMENTS OF MKMAT
C      KSAME = 4  CONTINUUM-CONTINUUM OFF-DIAGONAL ELEMENTS OF MKMAT
C
      IF(IBUG2.GT.2) WRITE(6,"(/' READ2B:',
     * 12X,'EVALUATION OF MAGNETIC INTEGRALS',9X,'KSAME =',I2)") KSAME
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
*     IF(JRELOP(4).EQ.0) GO TO 100
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      IF(IBUG2.GT.3) print "(20X,'testing READ2: MUL-,K- =',5I3,7I2)",
     * MULDSO,MULESO,MULDSS,MULDSP,KD1,KD2,KE1,KE2,KDS1,KDS2,KDS3,KDS4
C
      mtyp='      '
      IF(MULDSO.EQ.0) GO TO 2
      type='MULDSO'
      L=0
      DO 9 JJ=KD1,KD2,2
       L=L+1
       K=JJ-1
       DO IG=1,2
        DO J=1,JLIMIT
         DO I=1,ILIMIT
          RNKMAT(I,J,IG)=0.
          RNKMA2(I,J,IG)=0.
          VKMAT(I,J,IG)=0.
         ENDDO
        ENDDO
       ENDDO
C
C      GO TO (10,11,12,800),IK:
       CALL FINCGX(1,K,0)
C cor  A=D00N2(L) ! '03Nov10: why not INSIDE the nested loops!! etc!
       DO J=1,JLIMIT
        DO 8 I=1,ILIMIT
c  OR    M=I ! since off-diagonal elements can be 0
         IF ((KSAME.EQ.3 .AND. I.NE.J) .OR.
     *       (KSAME.EQ.4 .AND. I.EQ.J)) GO TO 8
         IBUG0 = I+J.eq.2 .and. IBUG2.gt.3
         IF(K) 9,17,18
   18    A=D00N2(L)
         IF(ABS(A).GE.EPS) THEN
          SOO(I,J) = SOO(I,J)+A*RNKMA2(I,J,2)
          CALL FAILMK(IBUG0,1,RNKMA2(1,1,2),TYPE,K,'RNKMA2',A,0.)
c                  OR     IF (RNKMA2(I,M,2).EQ.0.) GO TO 999
         ENDIF
         A=D00VK(L)
         IF(ABS(A).GE.EPS) THEN
          SOO(I,J) = SOO(I,J)+A*VKMAT(I,J,1)
          CALL FAILMK(IBUG0,1,VKMAT(1,1,1),TYPE,K,'-VKMAT',A,0.)
         ENDIF
         IF(M1.EQ.0.OR.M2.EQ.0) GO TO 17
          A=D11N2(L)
          IF(ABS(A).GE.EPS) THEN
           SOO(I,J) = SOO(I,J)+A*RNKMA2(I,J,1)
c          IF(IBUG0) print 99, 1,k, rnkma2(1,1,1), A;  mtyp='RNKMA2'
c          IF (RNKMA2(1,1,1).EQ.0.) GO TO 999
           CALL FAILMK(IBUG0,1,RNKMA2(1,1,1),TYPE,K,'RNKMA2',A,0.)
          ENDIF
          A=D11VK(L)
          IF(ABS(A).LT.EPS) GO TO 17
           SOO(I,J) = SOO(I,J)+A*VKMAT(I,J,2)
           CALL FAILMK(IBUG0,1,VKMAT(1,1,2),TYPE,K,' VKMAT',A,0.)
   17    A=D00NK(L)
         IF(ABS(A).GE.EPS) THEN
          SOO(I,J) = SOO(I,J)+A*RNKMAT(I,J,1)
          CALL FAILMK(IBUG0,1,RNKMAT(1,1,1),TYPE,K,'RNKMAT',A,0.)
         ENDIF
c cor    IF(M1.EQ.0.OR.M2.EQ.0) GO TO 9
         IF(M1*M2.NE.0) THEN
          A=D11NK(L)
c cor     IF(ABS(A).LT.EPS) GO TO 9
          IF(ABS(A).GE.EPS) THEN
           SOO(I,J) = SOO(I,J)+A*RNKMAT(I,J,2)
           CALL FAILMK(IBUG0,1,RNKMAT(1,1,2),TYPE,K,'RNKmat',A,0.)
          ENDIF
         ENDIF
    8    CONTINUE
        ENDDO
    9  CONTINUE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
    2 IF(MULESO.EQ.0) GO TO 100
      type='MULES0'
      L=0
      DO 43 JJ=KE1,KE2,2
       L=L+1
       K=JJ-1
       DO 47 IG=3,4
        DO 47 J=1,JLIMIT
         DO 47 I=1,ILIMIT
          RNKMAT(I,J,IG)=0.
          RNKMA2(I,J,IG)=0.
   47     VKMAT(I,J,IG)=0.
C      GO TO (57,58,50,801),IK:
       CALL FINCGX(3,K,0)
       A01K=E01NK(L)
       A01N=E01N2(L)
       A01V=E01VK(L)
       A10K=E10NK(L)
       A10N=E10N2(L)
       A10V=E10VK(L)
       DO J=1,JLIMIT
        DO 50 I=1,ILIMIT
         IF ((KSAME.EQ.3 .AND. I.NE.J) .OR.
     *       (KSAME.EQ.4 .AND. I.EQ.J)) GO TO 50
         IBUG0 = I+J.eq.2 .and. IBUG2.gt.3
         A = 0.
         IF(K) 43,44,45
   45     IF(M2.EQ.0) GO TO 46
           IF(ABS(A01N).GE.EPS) THEN
            A = A01N*RNKMA2(I,J,3)
c           IF(IBUG0) print 99, 2,k, rnkma2(1,1,3), A01N,A
c           IF (RNKMA2(1,1,3).EQ.0.) GO TO 999
            CALL FAILMK(IBUG0,2,RNKMA2(1,1,3),TYPE,K,'RNKma2',A01N,A)
           ENDIF
           IF(ABS(A01V).LT.EPS) GO TO 46
            A = A01V*VKMAT(I,J,4) + A
            CALL FAILMK(IBUG0,3,VKMAT(1,1,4),TYPE,K,'vkmat',A01V,A)
   46     IF(M1.EQ.0) GO TO 44
           IF(ABS(A10N).GE.EPS) THEN
            A = A10N*RNKMA2(I,J,4) + A
            CALL FAILMK(IBUG0,4,RNKMA2(1,1,4),TYPE,K,'rnkma2',A10N,A)
           ENDIF
           IF(ABS(A10V).LT.EPS) GO TO 44
            A = A10V*VKMAT(I,J,3) + A
            CALL FAILMK(IBUG0,5,VKMAT(1,1,3),TYPE,K,'vkmat',A10V,A)
   44     IF(M2.EQ.0) GO TO 49
           IF(ABS(A01K).LE.EPS) GO TO 49
            A = A01K*RNKMAT(I,J,4) + A
            CALL FAILMK(IBUG0,6,RNKMAT(1,1,4),TYPE,K,'rnkmat',A01K,A)
   49     IF(M1.EQ.0) GO TO 48 ! cor'03Nov12: 43 --> 48
           IF(ABS(A10K).GE.EPS) THEN
            A = A10K*RNKMAT(I,J,3) + A
            CALL FAILMK(IBUG0,7,RNKMAT(1,1,3),TYPE,K,'RNKmat',A10K,A)
           ENDIF
   48    SOO(I,J) = A + SOO(I,J)
   50    CONTINUE
        ENDDO
   43  CONTINUE
C
C 100 DO J=1,JLIMIT; DO I=1,ILIMIT
C       SOO(I,J)=SOO(I,J) ! *SIGNCH(=+1 Fano-Racah)  *CONSOO(=RAC2)
C      ENDDO; ENDDO
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*     IF(JRELOP(5).EQ.0) GO TO 200
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  100 IF(MULDSS.EQ.0) GO TO 80
      type='MULDSS'
      L=0
      DO 81 JJ=KDS1,KDS2,2
      L=L+1
      K=JJ-1
      DO J=1,JLIMIT
       DO I=1,ILIMIT
         RNKMAT(I,J,1)=0.
       ENDDO
      ENDDO
C     GO TO (83,84,86,802),IK:
      CALL FINCGX(1,K,1)
      A=D00SNK(L)
      IF(ABS(A).LT.EPS) GO TO 81
      CALL FAILMK(.FALSE.,0,RNKMAT(1,1,1),TYPE,K,'RnkMAT',A,0.)
      DO J=1,JLIMIT
       DO 79 I=1,ILIMIT
        IF ((KSAME.EQ.3 .AND. I.NE.J) .OR.
     *      (KSAME.EQ.4 .AND. I.EQ.J)) GO TO 79
        SS(I,J)=SS(I,J)+A*RNKMAT(I,J,1)
   79   CONTINUE
      ENDDO
      IF(IBUG2.gt.3) print 98, 1,k, rnkmat(1,1,1), A
   81 CONTINUE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
   80 IF(MULDSP.EQ.0) GO TO 88
      type='MULDSP'
      L=0
      DO 95 JJ=KDS3,KDS4,2
      L=L+1
      K=JJ-1
      DO J=1,JLIMIT
       DO I=1,ILIMIT
         RNKMAT(I,J,2)=0.
       ENDDO
      ENDDO
*     GO TO (90,91,89,803),IK:
      CALL FINCGX(1,K,2)
      A=D11SNK(L)
      IF(ABS(A).LT.EPS) GO TO 95
      CALL FAILMK(.FALSE.,0,RNKMAT(1,1,2),TYPE,K,'RnkMAT',A,0.)
      DO J=1,JLIMIT
       DO 94 I=1,ILIMIT
        IF ((KSAME.EQ.3 .AND. I.NE.J) .OR.
     *      (KSAME.EQ.4 .AND. I.EQ.J)) GO TO 94
        SS(I,J)=SS(I,J)+A*RNKMAT(I,J,2)
   94   CONTINUE
      ENDDO
      IF(IBUG2.gt.3) print 98, 2,k, rnkmat(1,1,2), A
   95 CONTINUE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
   88 IF(MULESS.EQ.0) GO TO 102
      type='MULESS'
      L=0
      DO 103 JJ=KES1,KES2,2
      L=L+1
      K=JJ-1
      DO J=1,JLIMIT
       DO I=1,ILIMIT
        RNKMAT(I,J,3)=0. ! marked '03Nov<_19, to ,4) Jun3, resto'04Jul11
       ENDDO
      ENDDO
*     GO TO (97,98,96,804),IK:
      CALL FINCGX(3,K,1)
      A = E01SNK(L)
      IF(ABS(A).LT.EPS) GO TO 103
      mtyp='rnkmat'
C cor if (RNKMAT(1,1,4).eq.0) go to 999 ! missing magn.int. till 04Jul11
c     if (RNKMAT(1,1,3).eq.0) go to 999:
      CALL FAILMK(.FALSE.,2,RNKMAT(1,1,3),TYPE,K,'rnkmat',A,0.)
      DO J=1,JLIMIT
       DO 89 I=1,ILIMIT
        IF ((KSAME.EQ.3 .AND. I.NE.J) .OR.
     *      (KSAME.EQ.4 .AND. I.EQ.J)) GO TO 89
C cor   SS(I,J)=SS(I,J)+A*RNKMAT(I,J,4)
        SS(I,J)=SS(I,J)+A*RNKMAT(I,J,3)
   89   CONTINUE
      ENDDO
      IF(IBUG2.gt.3) print 98, 4,k, rnkmat(1,1,3), A
  103 CONTINUE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  102 IF(MULESP.EQ.0) GO TO 200
      type='MULESP'
      L=0
      DO 111 JJ=KES3,KES4,2
      L=L+1
      K=JJ-1
      DO J=1,JLIMIT
       DO I=1,ILIMIT
         RNKMAT(I,J,4)=0.
       ENDDO
      ENDDO
*     GO TO (107,108,106,805),IK:
      CALL FINCGX(3,K,2)
      A=E10SNK(L)
      IF(ABS(A).LT.EPS) GO TO 111
* hm! if (RNKMAT(1,1,3).eq.0) go to 999 ! '04Jul11:
c     if (RNKMAT(1,1,4).eq.0) go to 999:
      CALL FAILMK(.FALSE.,2,RNKMAT(1,1,4),TYPE,K,'rnkMAT',A,0.)
      DO J=1,JLIMIT
       DO 109 I=1,ILIMIT
        IF ((KSAME.EQ.3 .AND. I.NE.J) .OR.
     *      (KSAME.EQ.4 .AND. I.EQ.J)) GO TO 109
* hm!   SS(I,J)=SS(I,J)+A*RNKMAT(I,J,3)
        SS(I,J)=SS(I,J)+A*RNKMAT(I,J,4)
  109   CONTINUE
      ENDDO
      IF(IBUG2.gt.3) print 98, 3,k, rnkmat(1,1,4), A
  111 CONTINUE
C 200 DO J=1,JLIMIT; DO I=1,ILIMIT
C       SS(I,J)=SS(I,J) ! *SIGNCH(=1 for Fano-Racah) *CONSS(=RAC2)
C      ENDDO; ENDDO
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  200 IF(IBUG2.LE.2) RETURN
       IF(ILIMIT+JLIMIT.LT.3) THEN
        print "(/20X,'SOO-,SS(1,1) =',1P,2E12.4E1)", SOO(1,1),SS(1,1)
       ELSE
        print"(9X,'SOO-,SS(1:2,1:2) =',1P,(T29,2E12.4E1,',',2E12.4E1))",
     *    SOO(1,1),SOO(1,2),SS(1,1),SS(1,2),
     *    SOO(2,1),SOO(2,2),SS(2,1),SS(2,2)
      ENDIF
      RETURN
   98 format(15x,'test: IX,K, SS,  A  =',2I3,1X,1P,2(E11.4E1,E12.4))
      END
C                        ************************
      SUBROUTINE FAILMK(IBUG0,IX,A,TYPE,K,MTYP,C,S)
      CHARACTER*6 TYPE,MTYP
      LOGICAL IBUG0
      M=1
      IF(S.EQ.0.) M=0
      IF(IBUG0) PRINT "(9x,'READ2 test: I,K, SOO,C,rt =',2I3,1X,
     *  1P,3E12.4)", IX, K, A,C,(S,L=1,M)
      IF(A.NE.0.) RETURN !!!! disabled for BB KSAME=3 tst '06Jun-Jul:
C tst IF(A.NE.9999999.) RETURN
      print"(/' READ2 stop: missing magnetic integral  ',A6,' K =',I2,
     * 2X,A6,', C =',1P,E12.5)", type,K,mtyp,C
      IF(ABS(C).GT.0.0002) STOP
      RETURN
      END
C***********************************************************************
      SUBROUTINE FINCGX(ID,K,ISS)
C
c     --- deduced from cgx modifications of recupd routines FIN2xx
c     in STG2bp and repetitive code in READ2B; 998-999 code '04Jun19/Jul
C
      CHARACTER LIT(0:7)*1
      COMMON/DEBUG/IBUG1,IBUG2,IBUG3,IBUG4,IBUG5,IBUG6,IBUG7,IBUG8,IBUG9
      COMMON  /MKMATX/KSAME,ILIMIT,JLIMIT
     */MKMATS/RNKMAT(  60,  60,4),RNKMA2(  60,  60,4),VKMAT(  60,  60,4)
      COMMON /NJLJ/NRHO,LRHO,NSIG,LSIG,NRHOP,LRHOP,NSIGP,LSIGP
      DIMENSION NA(4),LA(4)
      SAVE NA,LA, KAM, VSAVE
      DATA NA,LA,KAM/9*-1/, LIT/'s','p','d','f','g','h','i','j'/
C
      NSUM = NRHO+NRHOP + NSIG+NSIGP
      LAM = K
      IG0 = ISS/2+ID
   15  IG = IG0
   16   N1=NRHO
        L1=LRHO
        N2=NSIG
        L2=LSIG
        N3=NRHOP
        L3=LRHOP
        N4=NSIGP
        L4=LSIGP
        if(IG.eq.2.or.IG.eq.3 .or. (IG.eq.4.and.ID.eq.2))
     *                         call INTECH(N1,L1,N2,L2,0)
        if(IG.eq.2.or.IG.eq.4) call INTECH(N3,L3,N4,L4,0)
        KK = 1-ABS(MOD(L1+L3+LAM,2))*2
        IF (KK.LT.0 .AND. LAM.LT.0) GO TO 13
        IF(IBUG2.GT.2)
     *  print "(11X,'FINCGX: ID,K,ISS,IG,LAM =',5I3,'  at',4(I5,A1))",
     *     ID,K,ISS,IG,LAM, N1,LIT(L1),N2,LIT(L2), N3,LIT(L3),N4,LIT(L4)
        KGO = 1
        IF (NSUM.GT.999) THEN
         KGO = 2
         IF(NSUM.GT.1996) KGO=3
         IF(N4.ge.998) N4=0
         IF(N3.ge.998) N3=0
         IF(N2.ge.998) N2=0
         IF(N1.ge.998) N1=0
        ENDIF
C
C tst   GO TO (10,14,14),KGO !!!! tempmodd for BB KSAME=3 tst '06Jun-Jul
        GO TO (10,11,12),KGO ! IK; 800 CALL FIN2CB(ID,K,ISS);
C  10   CALL FIN2BB(ID,K,ISS):
   10   IF ((NA(1).ne.N1.or.NA(2).ne.N2.or.NA(3).ne.N3.or.NA(4).ne.N4)
     *  .or.(LA(1).ne.L1.or.LA(2).ne.L2.or.LA(3).ne.L3.or.LA(4).ne.L4)
     *  .or. KAM.ne.LAM) THEN
         CALL FIN2BB(N1,L1,N2,L2, N3,L3,N4,L4, LAM, VAL)
         VSAVE = VAL
         NA(1) = N1
         LA(1) = L1
         NA(2) = N2
         LA(2) = L2
         NA(3) = N3
         LA(3) = L3
         NA(4) = N4
         LA(4) = L4
         KAM = LAM
        ELSE
         VAL = VSAVE
        ENDIF
        if(KK.gt.0) then
          if(K.eq.LAM) then
C          RNKMAT(1,1,IG) = VAL:
           DO 5 I=1,ILIMIT
            DO 4 J=1,JLIMIT
    4        RNKMAT(J,I,IG) = 0.
    5       RNKMAT(I,I,IG) = VAL
          else if(ISS.eq.0) then
           DO 7 I=1,ILIMIT
            DO 6 J=1,JLIMIT
    6        RNKMA2(J,I,IG) = 0.
    7       RNKMA2(I,I,IG) = VAL
          endif
        else if(ISS.eq.0) then
          DO 9 I=1,ILIMIT
           DO 8 J=1,JLIMIT
    8       VKMAT(J,I,IG) = 0.
    9      VKMAT(I,I,IG) = VAL
        endif
        GO TO 13
C
C  11   CALL FIN2BC(ID,K,ISS):
   11   if(KK.gt.0) then
          if(LAM.eq.K) then
           CALL FIN2BC(N1,L1,N2,L2, N3,L3,N4,L4, LAM, RNKMAT(1,1,IG))
          else if(ISS.eq.0) then
           CALL FIN2BC(N1,L1,N2,L2, N3,L3,N4,L4, LAM, RNKMA2(1,1,IG))
          endif
        else if(ISS.eq.0) then
          CALL FIN2BC(N1,L1,N2,L2, N3,L3,N4,L4, LAM,  VKMAT(1,1,IG))
        endif
        GO TO 13
*  12   CALL FIN2CC(ID,K,ISS):
   12   if(KK.gt.0) then
          if(LAM.eq.K) then
           CALL FIN2CC(N1,L1,N2,L2, N3,L3,N4,L4, LAM, RNKMAT(1,1,IG))
          else if(ISS.eq.0) then
           CALL FIN2CC(N1,L1,N2,L2, N3,L3,N4,L4, LAM, RNKMA2(1,1,IG))
          endif
        else if(ISS.eq.0) then
          CALL FIN2CC(N1,L1,N2,L2, N3,L3,N4,L4, LAM,  VKMAT(1,1,IG))
        endif
C
   13   if(ISS.gt.0) go to 14 ! ???
        IG = IG+1
        IF(IG.LT.IG0+2) GO TO 16
       LAM = LAM-1
C
C      LAM downward selection rules (see also before label 10)
       IF (LAM.GE.-1 .AND. K-LAM.LE.2) GO TO 15
C
   14 RETURN
      END

