SUBROUTINE ADD0(JT1,N1,SR1,SI1,L11,L21,JT2,N2,SR2,SI2,L12,L22,
     1                JIN,JFN,IQ,AR,AI)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION SR1(N1),SI1(N1),L11(N1),L21(N1),
     1          SR2(N2),SI2(N2),L12(N2),L22(N2)
      DATA LAMBDA/0/
      FN(J1,J2)=DBLE((J1+J1+1)*(J2+J2+1))
C
      AR=0.D0
      AI=0.D0
      DO 1000 I1=1,N1
      DO 1000 I2=1,N2
      IF (L11(I1).NE.L12(I2).OR.L21(I1).NE.L22(I2)) GO TO 1000
      IF (L11(I1).EQ.L21(I1)) THEN
        DELTA=1.D0
      ELSE
        DELTA=0.D0
      ENDIF
      CR=DELTA-(SR1(I1)*SR2(I2)+SI1(I1)*SI2(I2))
      CI=SR1(I1)*SI2(I2)-SR2(I2)*SI1(I1)
      XL=SQRT(FN(L11(I1),L21(I1)))
      LL=L11(I1)+L21(I1)
      IF (LL-2*(LL/2).NE.0) XL=-XL
      XJ=X12J(JIN,JIN,L11(I1),L21(I1),JFN,L11(I1),JFN,L21(I1),
     1        IQ,JT2,JT1,LAMBDA)
      AR=AR+XL*XJ*CR
      AI=AI+XL*XJ*CI
 1000 CONTINUE
      FACT=FN(JT1,JT2)
      AR=AR*FACT
      AI=AI*FACT
      RETURN
      END
      SUBROUTINE ADD1(JT1,N1,SR1,SI1,L11,L21,JT2,N2,SR2,SI2,L12,L22,
     1                JIN,JFN,IQ,AR,AI)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION SR1(N1),SI1(N1),L11(N1),L21(N1),
     1          SR2(N2),SI2(N2),L12(N2),L22(N2)
      DATA LAMBDA/1/
      FN(J1,J2)=DBLE((J1+J1+1)*(J2+J2+1))
C
      AR=0.D0
      AI=0.D0
      DO 1000 I1=1,N1
      DO 1000 I2=1,N2
C     SKIP ON 3J TRIANGLE RELATIONSHIPS, ELSE EVALUATE 3J'S
      IF (ABS(L11(I1)-L12(I2)).GT.LAMBDA.OR.
     1    ABS(L21(I1)-L22(I2)).GT.LAMBDA) GO TO 1000
      TJ=THREEJ(L11(I1),L12(I2),LAMBDA)*THREEJ(L21(I1),L22(I2),LAMBDA)
      IF (L11(I1).EQ.L21(I1).AND.L12(I2).EQ.L22(I2)) THEN
        DELTA=1.D0
      ELSE
        DELTA=0.D0
      ENDIF
      CR=DELTA-(SR1(I1)*SR2(I2)+SI1(I1)*SI2(I2))
      CI=SR1(I1)*SI2(I2)-SR2(I2)*SI1(I1)
C     ADJUST FOR (-I)**(L-L'-L-BAR+L-BAR')
      LL=L11(I1)-L21(I1)-L12(I2)+L22(I2)
 2001 IF (LL.GT.0) GO TO 2002
      LL=LL+400
      GO TO 2001
 2002 LL=LL-4*(LL/4)
      IF (LL.EQ.1) THEN
        TT=CR
        CR=CI
        CI=-TT
      ELSEIF (LL.EQ.2) THEN
        CR=-CR
        CI=-CI
      ELSEIF (LL.EQ.3) THEN
        TT=CR
        CR=-CI
        CI=TT
      ENDIF
      XL=SQRT(FN(L11(I1),L21(I1))*FN(L12(I2),L22(I2)))
      XJ=X12J(JIN,JIN,L12(I2),L22(I2),JFN,L11(I1),JFN,L21(I1),
     1        IQ,JT2,JT1,LAMBDA)
      AR=AR+(TJ*XL*XJ)*CR
      AI=AI+(TJ*XL*XJ)*CI
 1000 CONTINUE
      FACT=FN(JT1,JT2)
      AR=AR*FACT
      AI=AI*FACT
      RETURN
      END
      SUBROUTINE SWRITE(IU,N,S)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION S(N,N)
C
      WRITE(IU) ((S(I,J),J=1,I),I=1,N)
      RETURN
C
      ENTRY SREAD(IU,N,S,IEND)
      IEND=0
      READ(IU,END=9999) ((S(I,J),J=1,I),I=1,N)
      DO 1000 I=1,N
      DO 1000 J=1,I-1
 1000 S(J,I)=S(I,J)
      RETURN
 9999 IEND=1
      RETURN
      END
      FUNCTION X12J(J1,J2,J3,J4,L1,L2,L3,L4,K1,K2,K3,K4)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C     CALCULATION OF 12J SYMBOL OF 2ND KIND AS DEFINED BY
C       YUTSIS ET AL., EQ. (19.3), P. 63
C     CALLS J6J (MOLSCAT SHULTEN-GORDON ALGORITHM) FOR 6J'S
C     VERSION 2. CHECKS FOR REPEATED 6J INDICES (SG, 22 NOV 94)
C
      PARAMETER (MXDIM=100)
      DIMENSION X6J1(MXDIM),X6J2(MXDIM),X6J3(MXDIM),X6J4(MXDIM)
      DATA TENTH/0.1D0/
C
C     CHECK FOR K4=0.  THIS WILL ARISE FROM LOU MONCHICK'S GHM EQS.
C     ONE MIGHT WANT TO CHECK FOR OTHER ZEROS AND REARRANGE
      IF (K4.EQ.0) THEN
C       USE YUTSIS EQ. (19.10)
        IF (J4.NE.L4.OR.J3.NE.L2) THEN
           X12J=0.D0
           RETURN
        ENDIF
        X12J=SIXJ(L1,L2,K1,K3,K2,J1)*SIXJ(K2,L3,K3,J2,L4,K1)/
     1       SQRT(DBLE((L2+L2+1)*(L4+L4+1)))
        LP=J1+L1-J2-L3
        IF (LP-2*(LP/2).NE.0) X12J=-X12J
        RETURN
      ENDIF
C
C     GENERAL CASE:  GET MIN, MAX VALUES OF SUMMATION INDEX
      IXMIN=MAX(ABS(K1-K2),ABS(K3-K4),ABS(J1-J3),ABS(J2-J4))
      IXMAX=MIN(K1+K2,K3+K4,J1+J3,J2+J4)
      IF (IXMIN.GT.IXMAX) THEN
C       12J MUST VANISH IN THIS CASE
        X12J=0.D0
        RETURN
      ENDIF
C
C     CONVERT EVERYONE TO DOUBLE PRECISION FOR INPUT TO J6J
      XJ1=J1
      XJ2=J2
      XJ3=J3
      XJ4=J4
      XL1=L1
      XL2=L2
      XL3=L3
      XL4=L4
      XK1=K1
      XK2=K2
      XK3=K3
      XK4=K4
C
C     GET 4 SETS OF 6J'S  -- CHECK TO SEE IF SOME ARE SAME
      IVAL1=MXDIM
      CALL J6J(XK1,XK2,XL1,XJ3,XJ1,IVAL1,XMIN1,X6J1)
      IF (K3.EQ.K1.AND.K4.EQ.K2.AND.L2.EQ.L1) THEN
        IVAL2=IVAL1
        XMIN2=XMIN1
        DO 2001 II=1,IVAL2
 2001   X6J2(II)=X6J1(II)
      ELSE
        IVAL2=MXDIM
        CALL J6J(XK3,XK4,XL2,XJ3,XJ1,IVAL2,XMIN2,X6J2)
      ENDIF
      IF (L3.EQ.L1.AND.J4.EQ.J3.AND.J2.EQ.J1) THEN
        IVAL3=IVAL1
        XMIN3=XMIN1
        DO 2002 II=1,IVAL3
 2002   X6J3(II)=X6J1(II)
      ELSEIF (J4.EQ.J3.AND.J2.EQ.J1.AND.L3.EQ.L2.AND.
     1        K1.EQ.K3.AND.K2.EQ.K4) THEN
        IVAL3=IVAL2
        XMIN3=XMIN2
        DO 2003 II=1,IVAL3
 2003   X6J3(II)=X6J2(II)
      ELSE
        IVAL3=MXDIM
        CALL J6J(XK1,XK2,XL3,XJ4,XJ2,IVAL3,XMIN3,X6J3)
      ENDIF
      IF (L4.EQ.L2.AND.J4.EQ.J3.AND.J2.EQ.J1) THEN
        IVAL4=IVAL2
        XMIN4=XMIN2
        DO 2004 II=1,IVAL4
 2004   X6J4(II)=X6J2(II)
      ELSEIF (K3.EQ.K1.AND.K4.EQ.K2.AND.L4.EQ.L3) THEN
        IVAL4=IVAL3
        XMIN4=XMIN3
        DO 2005 II=1,IVAL4
 2005   X6J4(II)=X6J3(II)
      ELSEIF (K3.EQ.K1.AND.K4.EQ.K2.AND.L4.EQ.L1.AND.
     1        J4.EQ.J3.AND.J2.EQ.J1) THEN
        IVAL4=IVAL1
        XMIN4=XMIN1
        DO 2006 II=1,IVAL4
 2006   X6J4(II)=X6J1(II)
      ELSE
        IVAL4=MXDIM
        CALL J6J(XK3,XK4,XL4,XJ4,XJ2,IVAL4,XMIN4,X6J4)
      ENDIF
C
C     SUM OVER IX=IXMIN,IXMAX
      SUM=0.D0
      DO 1000 IX=IXMIN,IXMAX
      IND1=1+IX-INT(XMIN1+TENTH)
      IF (IND1.LT.1.OR.IND1.GT.IVAL1) GO TO 1000
      IND2=1+IX-INT(XMIN2+TENTH)
      IF (IND2.LT.1.OR.IND2.GT.IVAL2) GO TO 1000
      IND3=1+IX-INT(XMIN3+TENTH)
      IF (IND3.LT.1.OR.IND3.GT.IVAL3) GO TO 1000
      IND4=1+IX-INT(XMIN4+TENTH)
      IF (IND4.LT.1.OR.IND4.GT.IVAL4) GO TO 1000
      SUM=SUM+X6J1(IND1)*X6J2(IND2)*X6J3(IND3)*X6J4(IND4)*DBLE(IX+IX+1)
 1000 CONTINUE
C     GET PARITY(L1-L2-L3+L4)
      LP=L1-L2-L3+L4
      IF (LP-2*(LP/2).NE.0) SUM=-SUM
      X12J=SUM
      RETURN
      END
      FUNCTION SIXJ(J1,J2,J5,J4,J3,J6)
C
C     CALCULATES 6-J SYMBOL:   _(J1 J2 J3 )_
C                               (J4 J5 J6 )
C     INTERFACE TO J6J ROUTINE.
C     MODIFIED BY S. GREEN 20 AUG 93; PASS DIMENSION OF XJ6J FOR CHECKING
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(MXDIM=200)
      DIMENSION XJ6J(MXDIM)
      IVAL=MXDIM
      CALL J6J(         DBLE(J2),DBLE(J3),
     1         DBLE(J4),DBLE(J5),DBLE(J6),
     3         IVAL,XJ1MIN,XJ6J)
      IND=1+J1-INT(XJ1MIN+0.1D0)
      SIXJ=0.D0
      IF(IND.GE.1 .AND. IND.LE.IVAL) SIXJ=XJ6J(IND)
      RETURN
      END
      SUBROUTINE J6J(J2,J3,L1,L2,L3,IVAL,J1MIN,D6J)
      IMPLICIT DOUBLE PRECISION (A-H,J-Z)
      DIMENSION D6J(2)
      DATA ZERO/0.D0/,TENTH/0.1D0/,HALF/0.5D0/,ONE/1.D0/,TWO/2.D0/,
     $ CONST/1.0D-12/
      E(J1S)=SQRT((J1S-MJ23S)*(J23P1S-J1S)*(J1S-ML23S)*(L23P1S-J1S))
      F(J1,JJP1)=(TWO*J1+ONE)*(JJP1*(FACT-JJP1-TWO*LLP1)+FACT2)
C
C  THIS ROUTINE CALCULATES THE 6-J COEFFICIENTS FOR ALL PERMISSIBLE
C  VALUES OF J1 FOR FIXED VALUES OF J2, J3, L1, L2, AND L3 USING THE
C  RECURSIVE ALGORITHM OF K. SCHULTEN AND R. G. GORDON, J. MATH. PHYS.
C  VOL. 16, P. 1961, (1975).
C  PROGRAMMED BY D. E. FITZ, 10/22/79
C  MODIFIED BY S. GREEN 20 AUG 93 TO TEST DIMENSION ON D6J
C
      MXDIM=IVAL
      JJP2=J2*(J2+ONE)
      JJP3=J3*(J3+ONE)
      LLP1=L1*(L1+ONE)
      LLP2=L2*(L2+ONE)
      LLP3=L3*(L3+ONE)
      MJ23S=(J2-J3)**2
      ML23S=(L2-L3)**2
      J23P1S=(J2+J3+ONE)**2
      L23P1S=(L2+L3+ONE)**2
      FACT2=(LLP2-LLP3)*(JJP2-JJP3)
      FACT=JJP2+JJP3+LLP2+LLP3
      J1MIN=MAX(ABS(J2-J3),ABS(L2-L3))
      J1MAX=MIN(J2+J3,L2+L3)
      IVAL=INT(J1MAX-J1MIN+ONE+TENTH)
      IF (IVAL.GT.MXDIM) THEN
        WRITE(6,*) 'J6J: ARRAY D6J TOO SMALL. NEEDS ',IVAL,' BUT ONLY ',
     1   MXDIM,' SUPPLIED'
        STOP
      ENDIF
C
C  TEST FOR OTHER TRIANGULAR INEQUALITES.
C
      IL1=INT(TWO*L1+TENTH)
      IL2=INT(TWO*L2+TENTH)
      IL3=INT(TWO*L3+TENTH)
      IJ2=INT(TWO*J2+TENTH)
      IJ3=INT(TWO*J3+TENTH)
      IF((IJ2.LE.IL1+IL3.AND.IJ2.GE.IABS(IL1-IL3)).AND.
     $ (IJ3.LE.IL1+IL2.AND.IJ3.GE.IABS(IL1-IL2))) GO TO 11
      DO 12 I=1,IVAL
   12 D6J(I)=ZERO
      RETURN
C
   11 INMID=(IVAL+3)/2
      SGNV=J2+J3+L2+L3
      SGN=ONE
      ISIGN=INT(SGNV+TENTH)
      IF(MOD(ISIGN,2).NE.0) SGN=-ONE
      D6J(1)=HALF
C
C  UPWARD RECURSION.
C
      IF(IVAL.EQ.1) GO TO 40
      JJP1=J1MIN*(J1MIN+ONE)
      F1=F(J1MIN,JJP1)
      J1=J1MIN+ONE
      J1S=J1*J1
      E2=E(J1S)
      IF(J1MIN.LT.TENTH) GO TO 15
      D6J(2)=-F1*D6J(1)/(E2*J1MIN)
      GO TO 16
   15 D6J(2)=-HALF*(LLP2+JJP2-LLP1)*D6J(1)/SQRT(JJP2*LLP2)
   16 SCALE=D6J(2)
      IF(IVAL.EQ.2) GO TO 40
      DO 21 IJ2=3,INMID
      JJP1=J1*(J1+ONE)
      F1=F(J1,JJP1)
      J1=J1+ONE
      E1=E2
      J1S=J1*J1
      E2=E(J1S)
   21 D6J(IJ2)=-(F1*D6J(IJ2-1)+J1*E1*D6J(IJ2-2))/(E2*(J1-ONE))
      SCALE=D6J(INMID)
      IEXC=5
      IF(ABS(SCALE).GT.CONST) GO TO 18
      INMID=INMID-1
      SCALE=D6J(INMID)
      IEXC=3
      GO TO 30
   18 IF(IVAL.EQ.3) GO TO 40
C
C  DOWNWARD RECURSION.
C
   30 D6J(IVAL)=HALF
      J1=J1MAX
      J1S=J1*J1
      JJP1=J1*(J1+ONE)
      F1=F(J1,JJP1)
      E1=E(J1S)
      D6J(IVAL-1)=-F1*D6J(IVAL)/(E1*(J1+ONE))
      IEND=IVAL-INMID
      IF(IVAL.LE.IEXC) GO TO 31
      DO 32 IJ2=2,IEND
      J1=J1-ONE
      E2=E1
      J1S=J1*J1
      JJP1=J1*(J1+ONE)
      E1=E(J1S)
      F1=F(J1,JJP1)
   32 D6J(IVAL-IJ2)=-(J1*E2*D6J(IVAL-IJ2+2)+F1*D6J(IVAL-IJ2+1))/
     $ (E1*(J1+ONE))
C
C  MATCH UPWARD AND DOWNWARD RECURSIVE RESULTS BY SCALING.
C
   31 SCALE=SCALE/D6J(INMID)
      DO 33 IJ2=INMID,IVAL
   33 D6J(IJ2)=SCALE*D6J(IJ2)
C
C  NORMALIZE RESULTS AND SET PHASE.
C
   40 SUM=ZERO
      DO 41 IJ2=1,IVAL
      J1=J1MIN+DBLE(IJ2-1)
   41 SUM=SUM+(TWO*J1+ONE)*D6J(IJ2)**2
      RNORM=ONE/SQRT(SUM*(TWO*L1+ONE))
      IF((SGN*D6J(IVAL)).LT.ZERO) RNORM=-RNORM
      DO 42 IJ2=1,IVAL
   42 D6J(IJ2)=D6J(IJ2)*RNORM
      RETURN
      END
      FUNCTION THREEJ (J1,J2,J3)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C     COMPUTATION OF SPECIAL WIGNER 3J COEFFICIENT WITH
C     VANISHING PROJECTIONS.  SEE EDMONDS, P. 50.
C
C     THIS VERSION EVALUATES BINOM AND PARITY IN-LINE
C       SHOULD IMPROVE EFFICIENCY, ESPECIALLY ON CRAY;
C       ALSO GIVES IMPROVEMENT ON AMDAHL  (SG: 20 DEC 92)
C
C     STATEMENT FUNCTION FOR DELTA ASSOCIATED W/ RACAH AND SIXJ SYMBOLS
C     DELTA(I,J,K)= SQRT(1.D0/ ( BINOM(I+J+K+1,I+J-K) *
C    1                 BINOM(K+K+1,I-J+K) * DBLE(K+J-I+1) )  )
C
      I1=J1+J2+J3
      IF (I1-2*(I1/2).NE.0)    GO TO 8
    1 I2=J1-J2+J3
      IF (I2) 8,2,2
    2 I3=J1+J2-J3
      IF (I3) 8,3,3
    3 I4=-J1+J2+J3
      IF (I4) 8,4,4
    4 I5=I1/2
      I6=I2/2
      SIGN=1.D0
      IF (I5-2*(I5/2).NE.0) SIGN=-SIGN
C   7 THREEJ=SIGN*DELTA(J1,J2,J3)*BINOM(I5,J1)*BINOM(J1,I6)
C     B1,B2 ARE BINOM ASSOCIATED W/ DELTA
      N=J1+J2+J3+1
      M=J1+J2-J3
      NM = N-M
      MNM = MIN(NM,M)
      IF(MNM.LE.0) THEN
        B1=1.D0
      ELSE
        FN = N+1
        F = 0.D0
        B = 1.D0
        DO 101 I = 1,MNM
        F = F+1.D0
        C = (FN-F)*B
  101   B = C/F
        B1 = B
      ENDIF
      N=J3+J3+1
      M=J1-J2+J3
      NM = N-M
      MNM = MIN(NM,M)
      IF(MNM.LE.0) THEN
        B2=1.D0
      ELSE
        FN = N+1
        F = 0.D0
        B = 1.D0
        DO 102 I = 1,MNM
        F = F+1.D0
        C = (FN-F)*B
  102   B = C/F
        B2 = B
      ENDIF
      DELTA=SQRT(1.D0/(B1*B2*(J3+J2-J1+1)))
C     B3=BINOM(I5,J1),  B4=BINOM(J1,I6)
      N=I5
      M=J1
      NM = N-M
      MNM = MIN(NM,M)
      IF(MNM.LE.0) THEN
        B3=1.D0
      ELSE
        FN = N+1
        F = 0.D0
        B = 1.D0
        DO 103 I = 1,MNM
        F = F+1.D0
        C = (FN-F)*B
  103   B = C/F
        B3 = B
      ENDIF
      N=J1
      M=I6
      NM = N-M
      MNM = MIN(NM,M)
      IF(MNM.LE.0) THEN
        B4=1.D0
      ELSE
        FN = N+1
        F = 0.D0
        B = 1.D0
        DO 104 I = 1,MNM
        F = F+1.D0
        C = (FN-F)*B
  104   B = C/F
        B4 = B
      ENDIF
      THREEJ=SIGN*DELTA*B3*B4
      RETURN
    8 THREEJ=0.D0
      RETURN
      END