      FUNCTION RD(X,Y,Z)
C
C     * NUMERICAL RECIPES ROUTINE USED TO EVALUATE INCOMPLETE
C     * ELLIPTIC INTEGRALS. CALLED BY GWLOOKU THROUGH FUNCTIONS
C     * EF KF 
C     * "PAUSE" STATEMENT IS REPLACED WITH WITH A CALL TO THE 
C     * "XIT" SUBROUTINE.
C     (C) COPR. 1986-92 NUMERICAL RECIPES SOFTWARE #=!E=#,)]UBCJ
C
C     * MAR 24/99 J SCINOCCA
C
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)
      REAL RD,X,Y,Z,ERRTOL,TINY,BIG,C1,C2,C3,C4,C5,C6

      PARAMETER (ERRTOL=0.05E0,TINY=1.E-50,BIG=4.5E42,
     1           C1=3.E0/14.E0,C2=1.E0/6.E0,C3=9.E0/22.E0,C4=3.E0/26.E0,
     2           C5=0.25E0*C3,C6=1.5E0*C4)

      REAL ALAMB,AVE,DELX,DELY,DELZ,EA,EB,EC,ED,EE,FAC,
     1     SQRTX,SQRTY,SQRTZ,SUM,XT,YT,ZT
C
      IF(MIN(X,Y).LT.0.E0.OR.MIN(X+Y,Z).LT.TINY
     1                 .OR.MAX(X,Y,Z).GT.BIG) THEN
        WRITE(6,'(A)') ' INVALID ARGUMENTS IN RD'
        CALL                                       XIT ('RD',-1)
      ENDIF

      XT=X
      YT=Y
      ZT=Z
      SUM=0.E0
      FAC=1.E0

    1 CONTINUE
        SQRTX=SQRT(XT)
        SQRTY=SQRT(YT)
        SQRTZ=SQRT(ZT)
        ALAMB=SQRTX*(SQRTY+SQRTZ)+SQRTY*SQRTZ
        SUM=SUM+FAC/(SQRTZ*(ZT+ALAMB))
        FAC=0.25E0*FAC
        XT=0.25E0*(XT+ALAMB)
        YT=0.25E0*(YT+ALAMB)
        ZT=0.25E0*(ZT+ALAMB)
        AVE=0.2E0*(XT+YT+3.E0*ZT)
        DELX=(AVE-XT)/AVE
        DELY=(AVE-YT)/AVE
        DELZ=(AVE-ZT)/AVE
      IF(MAX(ABS(DELX),ABS(DELY),ABS(DELZ)).GT.ERRTOL)GOTO 1

      EA=DELX*DELY
      EB=DELZ*DELZ
      EC=EA-EB
      ED=EA-6.E0*EB
      EE=ED+EC+EC
      RD=3.E0*SUM+FAC*(1.E0+ED*(-C1+C5*ED-C6*DELZ*EE)
     1                 +DELZ*(C2*EE+
     2                        DELZ*(-C3*EC+DELZ*C4*EA)))
     3               /(AVE*SQRT(AVE))

      RETURN
      END
