implicit real*8 (a-h,o-z) 
1     print*,' enter j1,j2,j3,m1,m2,m3'
      read(*,*)a1,a2,a3,b1,b2,b3
      print*,clebsh(a1,a2,a3,b1,b2,b3)
	go to 1
      end
C=====================================================================
C
      real*8 function threej(a1,a2,a3,b1,b2,b3)
C
C     Returns WIGNER 3J-COEFFICIENTS (for Clebsh coeficients change
C     name to clebsh and delete lines where b3=-b3)
C     a1,a2,a3 are j1,j2,j3 and b1,b2,b3 are m1,m2,m3.
C     array f(n+1)=n!/a**n, where a is an arbitrary number
C     that cancels in the calculation of clebsh; its purpose
C     is to prevent under/overflows in array f
C      
      implicit real*8 (a-h,o-z) 
      dimension f(85)
      save f
      data ifirst /1/
C
      if (ifirst .eq. 1) then
	f(1)=1.0
	do 10 i=1,84
10      f(i+1)=f(i)*dfloat(i)/15.
	ifirst=0
      endif
C
C  delete this line to run this program for clebsh coefficients      
      b3=-b3 
C
      clebsh=0. 
      threej=0.
      if(abs(b1)-a1 .gt. 0.01)  goto 100
      if(abs(b2)-a2 .gt. 0.01)  goto 100
      if(abs(b3)-a3 .gt. 0.01)  goto 100
      if(abs(b1+b2-b3).gt.0.01) goto 100
      if(a3-a1-a2 .gt. 0.01)    goto 100 
      if(abs(a1-a2)-a3 .gt. 0.01)  goto 100
      l1=a1+a2-a3+1.01
      l2=a1+b1+1.01
      l3=a2+b2+1.01
      l4=a2-b2+1.01
      l5=a3+a2-a1+1.01
      l6=a1-b1+1.01
      l7=a3+b3+1.01
      l8=a3-b3+1.01
      l9=a3+a1-a2+1.01
      l10=a3+a1+a2+1.01
      m1=l1
      m2=l6
      m3=l3
      m4=(l7+l8-l3-l4+l2-l6)/2+1
      m5=(l7+l8-l2-l6-l3+l4)/2+1
      m6=1
      zw=sqrt(f(l1)*f(l2)*f(l6))
      tw=sqrt(f(l9)/f(l10)*(2.0*a3+1.)/float(l10))
      sw=tw*sqrt(f(l3)*f(l7)*f(l8)*f(l4)*f(l5))
      min=max0(-m4,-m5,-m6)+2
      max=min0(m1,m2,m3)
      fr=-(-1.)**min
      if(max.lt.min)  return
      do m=min,max
	 mw=m-1
	 clebsh=clebsh+sw*fr/f(m4+mw)/(f(m2-mw)/zw)/f(m3-mw)/
     1                   f(m1-mw)/f(m5+mw)/f(m6+mw)
	 fr=-fr
      enddo
      threej=(-1)**nint(a1-a2+b3)/sqrt(2.*a3+1.)*clebsh
100   continue
C
C  delete this line to run this program for clebsh coefficients      
      b3=-b3 
C
      return
      end
C
C=====================================================================