SUBROUTINE CDGCFD (
*
*  inputs
*
     :                  A, X,
*
*  outputs
*
     :                  GAMMA, STATUS)
*
*  Module number:
*
*  Module name:
*
*  Keyphrase:
*  ----------
*  calculate the incomplete gamma function by its continued fraction
*  representation
*
*  Description:
*  ------------
*  adapted from Numerical Recipes by William Press et. al., p. 162,
*  subroutine GCF
*  Double precision version
*
*  FORTRAN name: CDGCFD.FOR
*
*  Keywords of accessed files and tables:
*  --------------------------------------
*  Name                         I/O     Description / Comments
*
*  Subroutines Called:
*  -------------------
*  CDBS:
*       CDGLN
*  SDAS:
*       UMSPUT
*  Others:
*       None
*
*  History:
*  --------
*  Version      Date        Author          Description
*     1       11-11-88     J.-C. HSU        coding
*-------------------------------------------------------------------------------
*
*== input:
*                                       --arguments of the incomplete gamma
*                                       --function
        DOUBLE PRECISION        A, X
*
*== output:
*                                       --functional value
        DOUBLE PRECISION        GAMMA
*                                       --error status
        INTEGER                 STATUS
*
*== local:
*                                       --convergence criterion
        DOUBLE PRECISION        EPS,
*                                       --natural log of gamma(A)
     :                          GAMLN,
     :                          G, GOLD, A0, A1, B0, B1, 
     :                          FAC, AN, ANA, ANF
*                                       --loop index
        INTEGER                 ITMAX, I,
     :                          STATOK
*                                       --error message
        CHARACTER*130           CONTXT, MESS
        PARAMETER               (ITMAX = 100, EPS = 1.D-12)
*=========================begin hsp.inc=========================================
*                                       --status return code
        INTEGER         OK, ERRNUM(20)
        INTEGER         DEST, PRIO
        DATA OK /0/
        DATA ERRNUM /701, 702, 703, 704, 705, 706, 707, 708, 709, 710,
     :               711, 712, 713, 714, 715, 716, 717, 718, 719, 720/
*                                       --message destination and priority
        DATA DEST, PRIO /1, 0/
*=========================end hsp.inc===========================================
*------------------------------------------------------------------------------
*
        STATUS = OK
        CALL CDGLN (A, GAMLN)
*
*  reset "old" value, used to test for convergence
*
        GOLD = 0.D0
*
*  setting up A's and B's for evaluating the continued fraction
*
        A0 = 1.D0
        A1 = X
        B0 = 0.D0
        B1 = 1.D0
*
*  renormalization factor to prevent overflow of the partial numerators
*  and denominators
*
        FAC = 1.D0
*
        DO 10 I = 1, ITMAX
            AN = DBLE (I)
            ANA = AN - A
            A0 = (A1 + A0 * ANA) * FAC
            B0 = (B1 + B0 * ANA) * FAC
            ANF = AN * FAC
            A1 = X * A0 + ANF * A1
            B1 = X * B0 + ANF * B1
*
            IF (A1 .NE. 0.D0) THEN
*
*  renormalize
*
                FAC = 1.D0 / A1
                G = B1 * FAC
*
*  converge?
*
                IF (ABS((G - GOLD)/G) .LT. EPS) GO TO 20
                GOLD = G
            END IF
   10   CONTINUE
*
        CONTXT = 'continued fraction evalution of the incomplete '
     :                  // 'gamma function does not converge'
        STATUS = ERRNUM(2)
        GO TO 999
*
   20   GAMMA = EXP (-X + A * LOG(X) - GAMLN) * G
        GO TO 1000
*
  999   MESS = 'CDGCFD: ' // CONTXT
        CALL UMSPUT (MESS, DEST, PRIO,    STATOK)
*
 1000   RETURN
        END