SUBROUTINE CDGSRD (
*
*  inputs
*
     :                  A, X,
*
*  outputs
*
     :                  GAMMA, STATUS)
*
*  Module number:
*
*  Module name:
*
*  Keyphrase:
*  ----------
*  calculate the incomplete gamma function by its series representation
*
*  Description:
*  ------------
*  adapted from Numerical Recipes by William Press et. al., p. 162,
*  subroutine GSER
*  Double precision version
*
*  FORTRAN name: CDGSRD.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,
     :                          AP, SUM, DEL
*                                       --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)
*
        IF (X .LE. 0) THEN
            IF (X .LT. 0) THEN
                CONTXT = 'second argument of the incomplete gamma '
     :                   // 'function is negative'
                STATUS = ERRNUM(1)
                GO TO 999
            ELSE
                GAMMA = 0.D0
                GO TO 1000
            END IF
        END IF
*
        AP = A
        SUM = 1.D0 / A
        DEL = SUM
*
        DO 10 I = 1, ITMAX
            AP = AP + 1.D0
            DEL = DEL * X / AP
            SUM = SUM + DEL
            IF (ABS(DEL) .LT. ABS(SUM)*EPS) GO TO 20
   10   CONTINUE
*
        CONTXT = 'series evalution of the incomplete gamma function '
     :                  // 'does not converge'
        STATUS = ERRNUM(2)
        GO TO 999
*
   20   GAMMA = SUM * EXP (-X + A * LOG(X) - GAMLN)
        GO TO 1000
*
  999   MESS = 'CDGSRD: ' // CONTXT
        CALL UMSPUT (MESS, DEST, PRIO,    STATOK)
*
 1000   RETURN
        END