!+ Function to calculate the gamma-function.
!
! Method:
!	The gamma function for the fractional part of the argument is
!	evaluated using Weierstrass's formula. The integral part is
!	handled by recurrence. A correction is applied to Weierstrass's
!	formula to reduce the influence of the truncation: this is
!	derived by imposing the condition that Gamma(1)=1. The 
!	coefficient in this truncation could be precalculated, but it
!	depends on N_LIMIT and is therefore calculated here so that
!	N_LIMIT may be changed easily and consistently. With 
!	N_LIMIT=10,000, the error should be less than one part in 10^7.
!	This routine is intended only for positive arguments and is
!	not especially efficient. 
!	
!
! Current Owner of Code: J. M. Edwards
!
! History:
!	 Version 	   Date	  		Comment
!   
!	 1.0		   12-04-95		First Version under RCS
!						(J. M. Edwards)
!
! Description of Code:
!   FORTRAN 77 with extensions listed in documentation. 
!
!- ---------------------------------------------------------------------
      FUNCTION GAMMA_FNC(IERR, Z)
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/error.pcf'
      INCLUDE 'incl/def_std_io.icf'
!
!     DUMMY ARGUMENTS
      INTEGER	!, INTENT(OUT)
     &     IERR
!             ERROR FLAG
      REAL*8	!, INTENT(IN)
     &     Z
!             POINT OF EVALUATION
      REAL*8	!, INTENT(OUT)
     &     GAMMA_FNC
!             GAMMA FUNCTION
!
!     LOCAL VARIABLES
      INTEGER
     &     N_LIMIT
!             LIMITING NUMBER OF TERMS
     &   , Z_INT
!             INTEGER PART OF Z
     &   , I
!             LOOP VARAIBLE
      REAL*8
     &     EXCLUDE_ZERO
!             REGION TO EXCLUDE NEAR 0.0
     &   , EULER
!             EULER'S CONSTANT
     &   , Z_REDUCED
!             REDUCED PART OF Z
     &   , U
!             TEMPORARY VARIABLE
     &   , V
!             TEMPORARY VARIABLE
     &   , U_1
!             TEMPORARY VARIABLE
     &   , V_1
!             TEMPORARY VARIABLE
     &   , CORRECTION
!             CORRECTION TO ALLOW FOR
!             TRUNCATION
!
!     DEFINE EULER'S CONSTANT
      DATA EULER/5.77215665D-01/
!     DEFINE NUMERICAL CONSTANTS.
      DATA EXCLUDE_ZERO/5.0D-02/
      DATA N_LIMIT/10000/
!
!
!     REDUCE THE ARGUMENT TO A VALUE IN (0,1)
      IF (Z.LE.0.0D+00) THEN
         WRITE(IU_ERR, '(/A)') 
     &      '***ERROR: THE GAMMA FUNCTION CANNOT '
     &      //'BE EVALUATED BELOW 0.0.'
         IERR=I_ERR_FATAL
         RETURN
      ENDIF
!
      Z_INT=INT(Z)
      Z_REDUCED=Z-Z_INT
!     THERE WILL BE NUMERICAL DIFFICULTIES NEAR 0, 
!     SO WE EXCLUDE SUCH CASES.
      IF (Z_REDUCED.LT.EXCLUDE_ZERO) THEN
         Z_REDUCED=Z_REDUCED+1.0D+00
         Z_INT=Z_INT-1
      ENDIF
      U=0.0D+00
      U_1=0.0D+00
      V=1.0D+00
      V_1=1.0D+00
      DO I=1, N_LIMIT
         U=U+Z_REDUCED/DBLE(I)
         V=V*(1.0D+00+Z_REDUCED/DBLE(I))
         U_1=U_1+1.0D+00/DBLE(I)
         V_1=V_1*(1.0D+00+1.0D+00/DBLE(I))
      ENDDO
      CORRECTION=1.0D+00+Z_REDUCED**2
     &   *(1.0D+00-EXP(-EULER)/(V_1*EXP(-U_1)))
      GAMMA_FNC=1.0D+00/(Z_REDUCED*V*EXP(EULER*Z_REDUCED-U))
      GAMMA_FNC=GAMMA_FNC*CORRECTION
!
!     CONVERT FROM THE REDUCED FUNCTION TO THE FULL FUNCTION.
      IF (Z_INT.GE.1) THEN
         DO I=1, Z_INT
            GAMMA_FNC=GAMMA_FNC*(Z-DBLE(I))
         ENDDO
      ELSE IF(Z_INT.EQ.-1) THEN
!        WHEN Z IS SMALL Z_REDUCED WILL HAVE BEEN RAISED TO AVOID
!        ILL-CONDITIONING, SO THE CONVERSION GOES THE OTHER WAY.
         GAMMA_FNC=GAMMA_FNC/Z
      ENDIF
!
!
      RETURN
      END
