SUBROUTINE CDGAMD (
*
*  inputs
*
     :                  A, X,
*
*  outputs
*
     :                  GAMMA, STATUS)
*
*  Module number:
*
*  Module name:
*
*  Keyphrase:
*  ----------
*  calculate the incomplete gamma function
*
*  Description:
*  ------------
*  adapted from Numerical Recipes by William Press et. al., p. 161,
*  subroutine GAMMP
*  Double precision version
*
*  FORTRAN name: CDGAMD.FOR
*
*  Keywords of accessed files and tables:
*  --------------------------------------
*  Name                         I/O     Description / Comments
*
*  Subroutines Called:
*  -------------------
*  CDBS:
*       CDGSRD, CDGCFD
*  SDAS:
*       UMSPUT
*  Others:
*       None
*
*  History:
*  --------
*  Version      Date        Author          Description
*     1       11-11-88     J.-C. HSU        coding
*-------------------------------------------------------------------------------
*
*== input:
*                                       --arguments of the gamma function
        DOUBLE PRECISION        A, X
*
*== output:
*                                       --functional value
        DOUBLE PRECISION        GAMMA
*                                       --error status
        INTEGER                 STATUS
*
*== local:
*                                       --error message
        CHARACTER*130           CONTXT, MESS
        INTEGER                 STATOK
*=========================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
*
        IF (X .LT. 0.D0 .OR. A .LE. 0.D0) THEN
            CONTXT = 'negative second argument or non-positive first '
     :                  // 'argument'
            STATUS = ERRNUM(1)
            GO TO 999
        END IF
*
*  use the series representation
*
        IF (X .LT. A+1.D0) THEN
            CALL CDGSRD (A, X,    GAMMA, STATUS)
        ELSE
*
*  use the continued fraction representation
*
            CALL CDGCFD (A, X,    GAMMA, STATUS)
            GAMMA = 1.D0 - GAMMA
        END IF
*
        IF (STATUS .NE. OK) THEN
            CONTXT = 'cannot evaluate gamma function'
            GAMMA = 0.D0
            GO TO 999
        END IF
*
        GO TO 1000
*
  999   MESS = 'CDGAMD: ' // CONTXT
        CALL UMSPUT (MESS, DEST, PRIO,    STATOK)
*
 1000   RETURN
        END