C **********************************************
SUBROUTINE TCGCOF(C,JA,JB,JC,MA,MB,IND,IPRINT,J1,J2,J3)
C CLEBSCH-GORDON COEFFICIENT
C SEE QUANTUM MECHANICS BY A.MESSIAH EQ.(C21) AND (C12)
DOUBLE PRECISION C,C1,C2,C3,WA,WB,W1,W2,W3
C CHARACTER A,MSR(64)
CHARACTER A*8,MSR(64)*8
DIMENSION IST(64),MSQ(64),MSP(64),NSS(64,2),LC(2),LM(2)
DIMENSION A(6),J(3),M(3),IB(3,3),IA(5)
DATA (IST(I),I=1,10)/2,3,5,7,11,13,17,19,23,29/
DATA KETA,KETB/64,64/
C=0.0
IND=0
MT=10
DO 12 IP=1,MT
12 NSS(IP,1)=0
LC(1)=0
LM(1)=0
JQ=JA+JB+JC+2
JP=JQ/2
IF(JP*2.NE.JQ) GO TO 3
J(1)=JA
J(2)=JB
J(3)=JC
M(1)=MA
M(2)=MB
M(3)=MA+MB
DO 1 K=1,3,1
IQ=J(K)+M(K)
IB(1,K)=IQ/2
IF(IB(1,K)*2.NE.IQ) GO TO 3
IF(IB(1,K).LT.0) GO TO 2
IB(2,K)=(J(K)-M(K))/2
IF(IB(2,K).LT.0) GO TO 2
KA=MOD(K,3)+1
KB=MOD(K+1,3)+1
IB(3,K)=(J(K)+J(KA)-J(KB))/2
IF(IB(3,K).LT.0) GO TO 2
1 CONTINUE
IA(1)=IB(3,1)
IA(2)=IB(2,1)
IA(3)=IB(1,2)
IA(4)=(JB-JC-MA)/2
IA(5)=(JA-JC+MB)/2
KA=MIN0(IA(1),IA(2),IA(3))+1
KB=MAX0(IA(4),IA(5),0)+1
IF(JP.LE.IST(MT)) GO TO 91
NN=JP
MTA=MT
NT=IST(MT)+1
DO 92 NNI=NT,NN
DO 93 I=1,MTA
IF(MOD(NNI,IST(I)).EQ.0) GO TO 92
93 CONTINUE
MTA=MTA+1
IF(MTA.GT.KETA) GO TO 94
IST(MTA)=NNI
92 CONTINUE
NT=MT+1
DO 95 I=NT,MTA
95 NSS(I,1)=0
MT=MTA
91 CONTINUE
DO 5 NX=KB,KA,1
LM(2)=1-MOD(NX+1,2)*2
DO 6 IP=1,MT
6 NSS(IP,2)=0
LC(2)=0
III=NX-1
IF(III.LE.1) GO TO 37
DO 34 II=2,III
N=II
NT=1
L=0
DO 32 I=NT,MT
36 IF(N.EQ.1) GO TO 35
IF(MOD(N,IST(I)).NE.0) GO TO 33
NSS(I,2)=NSS(I,2)-1
N=N/IST(I)
IF(NSS(I,2).NE.0) L=I
GO TO 36
33 CONTINUE
IF(NSS(I,2).NE.0) L=I
32 CONTINUE
35 IF(I.GE.LC(2)) LC(2)=L
34 CONTINUE
37 CONTINUE
DO 7 KC=1,5
III=IABS(IA(KC)-NX+1)
IF(III.LE.1) GO TO 47
DO 44 II=2,III
N=II
NT=1
L=0
DO 42 I=NT,MT
46 IF(N.EQ.1) GO TO 45
IF(MOD(N,IST(I)).NE.0) GO TO 43
NSS(I,2)=NSS(I,2)-1
N=N/IST(I)
IF(NSS(I,2).NE.0) L=I
GO TO 46
43 CONTINUE
IF(NSS(I,2).NE.0) L=I
42 CONTINUE
45 IF(I.GE.LC(2)) LC(2)=L
44 CONTINUE
47 CONTINUE
7 CONTINUE
IF(LM(1).NE.0) GO TO 9
DO 8 IP=1,MT
8 NSS(IP,1)=NSS(IP,2)
LM(1)=LM(2)
LC(1)=LC(2)
DO 811 JLP=2,KETB
811 MSP(JLP)=0
MSP(1)=1
NNP=1
GO TO 5
9 CONTINUE
DO 812 JLQ=2,KETB
812 MSQ(JLQ)=0
MSQ(1)=1
NNQ=1
LB=MAX0(LC(1),LC(2))
IF (LB.EQ.0) GO TO 888
LL=0
DO 88 IP=1,LB
LA=NSS(IP,1)-NSS(IP,2)
IF(LA) 82,81,83
82 LA=-LA
DO 821 ILQ=1,LA
DO 822 JLQ=1,NNQ
822 MSQ(JLQ)=MSQ(JLQ)*IST(IP)
DO 823 JLQ=1,NNQ
JWQ=MSQ(JLQ)/10000
IF(JWQ.EQ.0) GO TO 823
MSQ(JLQ)=MSQ(JLQ)-JWQ*10000
MSQ(JLQ+1)=MSQ(JLQ+1)+JWQ
823 CONTINUE
IF(JWQ.EQ.0) GO TO 821
NNQ=NNQ+1
IF(NNQ.GE.KETB) GO TO 891
821 CONTINUE
GO TO 81
83 NSS(IP,1)=NSS(IP,2)
DO 831 ILP=1,LA
DO 832 JLP=1,NNP
832 MSP(JLP)=MSP(JLP)*IST(IP)
DO 833 JLP=1,NNP
JWP=MSP(JLP)/10000
IF(JWP.EQ.0) GO TO 833
MSP(JLP)=MSP(JLP)-JWP*10000
MSP(JLP+1)=MSP(JLP+1)+JWP
833 CONTINUE
IF(JWP.EQ.0) GO TO 831
NNP=NNP+1
IF(NNP.GE.KETB) GO TO 892
831 CONTINUE
81 CONTINUE
IF(NSS(IP,1).NE.0) LL=IP
88 CONTINUE
LC(1)=LL
888 CONTINUE
NNR=MAX0(NNP,NNQ)+1
NRR=NNR-1
ISP=LM(1)*LM(2)
ILP=0
DO 801 JLP=1,NRR
JWP=MSP(JLP)+MSQ(JLP)*ISP
IF(JWP.LT.0) GO TO 802
IF(JWP.LT.10000) GO TO 803
JWP=JWP-10000
MSP(JLP+1)=MSP(JLP+1)+1
GO TO 803
802 JWP=JWP+10000
MSP(JLP+1)=MSP(JLP+1)-1
803 MSP(JLP)=JWP
IF(JWP.EQ.0) GO TO 801
IF(ILP.EQ.0) ILQ=JLP
ILP=JLP
801 CONTINUE
IF(ILP.NE.0) GO TO 804
IF(MSP(NNR).NE.0) GO TO 804
LM(1)=0
LC(1)=0
DO 89 IP=1,MT
89 NSS(IP,1)=0
GO TO 5
804 IF(MSP(NNR).LT.0) GO TO 805
NNP=ILP
IF(MSP(NNR).GT.0) NNP=NNR
IF(NNP.GE.KETB) GO TO 893
GO TO 840
805 LM(1)=-LM(1)
MSP(ILQ)=MSP(ILQ)-1
ILP=0
DO 806 JLP=ILQ,NRR
MSP(JLP)=9999-MSP(JLP)
IF(MSP(JLP).NE.0) ILP=JLP
806 CONTINUE
MSP(NNR)=0
NNP=ILP
840 IF((NNP.EQ.1).AND.(MSP(1).EQ.1)) GO TO 5
L=0
DO 22 I=1,MT
26 IF((NNP.EQ.1).AND.(MSP(1).EQ.1)) GO TO 25
NNQ=0
JWP=0
DO 841 JLP=1,NNP
ILP=NNP-JLP+1
JWR=MSP(ILP)+JWP*10000
JWQ=JWR/IST(I)
JWP=JWR-JWQ*IST(I)
MSQ(ILP)=JWQ
IF(JWQ.NE.0.AND.NNQ.EQ.0) NNQ=ILP
841 CONTINUE
IF(JWP.NE.0) GO TO 23
NSS(I,1)=NSS(I,1)+1
NNP=NNQ
DO 842 JLQ=1,NNQ
842 MSP(JLQ)=MSQ(JLQ)
IF(NSS(I,1).NE.0) L=I
GO TO 26
23 CONTINUE
IF(NSS(I,1).NE.0) L=I
22 CONTINUE
25 IF(I.GE.LC(1)) LC(1)=L
NNQ=NNP+1
DO 843 JLP=NNQ,KETB
843 MSP(JLP)=0
5 CONTINUE
IF(LM(1).EQ.0) GO TO 16
MMM=LC(1)
IF(MMM.EQ.0) GO TO 15
DO 10 IP=1,MMM
10 NSS(IP,1)=2*NSS(IP,1)
15 CONTINUE
N=JC+1
IF(N.EQ.1) GO TO 77
NT=1
L=0
DO 72 I=NT,MT
76 IF(N.EQ.1) GO TO 75
IF(MOD(N,IST(I)).NE.0) GO TO 73
NSS(I,1)=NSS(I,1)+1
N=N/IST(I)
IF(NSS(I,1).NE.0) L=I
GO TO 76
73 CONTINUE
IF(NSS(I,1).NE.0) L=I
72 CONTINUE
75 IF(I.GE.LC(1)) LC(1)=L
77 CONTINUE
III=JP
IF(III.LE.1) GO TO 57
DO 54 II=2,III
N=II
NT=1
L=0
DO 52 I=NT,MT
56 IF(N.EQ.1) GO TO 55
IF(MOD(N,IST(I)).NE.0) GO TO 53
NSS(I,1)=NSS(I,1)-1
N=N/IST(I)
IF(NSS(I,1).NE.0) L=I
GO TO 56
53 CONTINUE
IF(NSS(I,1).NE.0) L=I
52 CONTINUE
55 IF(I.GE.LC(1)) LC(1)=L
54 CONTINUE
57 CONTINUE
DO 11 KA=1,3,1
DO 11 KB=1,3,1
III=IB(KA,KB)
IF(III.LE.1) GO TO 67
DO 64 II=2,III
N=II
NT=1
L=0
DO 62 I=NT,MT
66 IF(N.EQ.1) GO TO 65
IF(MOD(N,IST(I)).NE.0) GO TO 63
NSS(I,1)=NSS(I,1)+1
N=N/IST(I)
IF(NSS(I,1).NE.0) L=I
GO TO 66
63 CONTINUE
IF(NSS(I,1).NE.0) L=I
62 CONTINUE
65 IF(I.GE.LC(1)) LC(1)=L
64 CONTINUE
67 CONTINUE
11 CONTINUE
16 CONTINUE
IND=1
INDA=0
IF(LM(1).EQ.0) GO TO 971
J2=1
J3=1
C2=1.0
C3=1.0
IF(NNP.LE.2) GO TO 992
INDA=1
J1=1
C1=0.0
DO 993 JLP=1,NNP,2
ILP=NNP-JLP+1
IF(ILP.EQ.1) GO TO 994
WA=1.0
IF(ILP.GT.2) WA=10000.0D0**(ILP-2)
WB=MSP(ILP)*10000+MSP(ILP-1)
C1=C1+WB*WA
GO TO 993
994 WB=MSP(ILP)
C1=C1+WB
993 CONTINUE
GO TO 995
992 J1=MSP(2)*10000+MSP(1)
C1=1.0
995 CONTINUE
LA=LC(1)
IF(LA.EQ.0) GO TO 972
DO 988 IP=1,LA
988 NSS(IP,2)=NSS(IP,1)
DO 981 IP=1,LA
JWP=MOD(IABS(NSS(IP,2)),2)
IF(JWP.EQ.0) GO TO 982
NSS(IP,2)=NSS(IP,2)+1
IF(FLOAT(J3)*FLOAT(IST(IP)).LT.9999000D1) GO TO 983
INDA=3
WB=J3
C3=WB*C3
J3=IST(IP)
GO TO 982
983 J3=J3*IST(IP)
982 LB=IABS(NSS(IP,2))/2
IF(LB.EQ.0) GO TO 981
DO 984 ILP=1,LB
IF(NSS(IP,1).LT.0) GO TO 985
IF(FLOAT(J1)*FLOAT(IST(IP)).LT.9999000D1) GO TO 986
INDA=1
WB=J1
C1=C1*WB
J1=IST(IP)
GO TO 984
986 J1=J1*IST(IP)
GO TO 984
985 IF(FLOAT(J2)*FLOAT(IST(IP)).LT.9999000D1) GO TO 987
INDA=2
WB=J2
C2=C2*WB
J2=IST(IP)
GO TO 984
987 J2=J2*IST(IP)
984 CONTINUE
981 CONTINUE
972 J1=J1*LM(1)
MSP(NNP)=MSP(NNP)*LM(1)
W1=J1
W2=J2
W3=J3
C=((C1/C2)*(W1/W2))/DSQRT(W3*C3)
IF(INDA.NE.0) IND=2
GO TO 973
3 IND=-1
2 CONTINUE
INDA=0
971 J1=0
J2=1
J3=1
C=0.0
973 CONTINUE
IF(IPRINT.LE.0) GO TO 974
MC=MA+MB
WRITE(6,101) JA,JB,JC,MA,MB,MC,C,IND
101 FORMAT(/4H CG(,3(I4,2H/2),2X,3(I4,2H/2),1H),D16.7,6H IND=,I2)
IF(IPRINT.LE.1) GO TO 974
IF(INDA.NE.0) GO TO 975
WRITE(6,102) J1,J2,J3
102 FORMAT(4H = (,I11,1H/,I10,7H)/SQRT(,I10,1H))
IF(IPRINT.LE.2) GO TO 974
975 IF(IND.LE.0) GO TO 974
MMM=LC(1)
IF(MMM.EQ.0) GO TO 974
DO 976 JLP=1,NNP
ILP=NNP-JLP+1
WRITE(MSR(JLP),106) MSP(ILP)
106 FORMAT(I5)
IF(ILP.EQ.NNP) GO TO 976
READ(MSR(JLP),107) (A(III),III=1,5)
107 FORMAT(5A1)
DO 977 III=2,5
IF(A(III).NE.' ') GO TO 978
A(III)='0'
977 CONTINUE
978 CONTINUE
WRITE(MSR(JLP),107) (A(III),III=1,5)
976 CONTINUE
WRITE(6,103) (MSR(ILP),ILP=1,NNP)
103 FORMAT(1H ,10A5)
WRITE(6,104) (NSS(IP,1),IP=1,MMM)
104 FORMAT(1H ,10I6)
974 RETURN
891 NSTOP=891
GO TO 894
892 NSTOP=892
GO TO 894
893 NSTOP=893
GO TO 894
94 NSTOP=94
894 WRITE(6,601) NSTOP
601 FORMAT(24H ERROR IN TCGCOF, NSTOP=,I5)
IND=-2
C=0.0
J1=0
J2=1
J3=1
RETURN
END