FUNCTION chyper(point, kk, ll, mm, nn, ifault) RESULT(fn_val)
 
! Code converted using TO_F90 by Alan Miller
! Date: 2002-09-24  Time: 22:57:18

!     Cumulative hypergeometric probabilities

!     ALGORITHM AS R77  APPL. STATIST. (1989), VOL.38, NO.1
!     Replaces AS 59 and AS 152
!     Incorporates AS R86 from vol.40(2)

!     Auxiliary routines required: ALNFAC (AS 245), ALNORM (AS 66)

IMPLICIT NONE
INTEGER, PARAMETER    :: dp = SELECTED_REAL_KIND(15, 100)

LOGICAL, INTENT(IN)   :: point
INTEGER, INTENT(IN)   :: kk
INTEGER, INTENT(IN)   :: ll
INTEGER, INTENT(IN)   :: mm
INTEGER, INTENT(IN)   :: nn
INTEGER, INTENT(OUT)  :: ifault
REAL(dp)              :: fn_val

INTEGER   :: k, l, m, n, i, j, nl, kl, mnkl,
REAL(dp)  :: p, pt,  alnfac,  mean, sig, alnorm, arg
LOGICAL   :: dir
REAL(dp), PARAMETER  :: zero = 0.0_dp, half = 0.5_dp, one = 1.0_dp
INTEGER, PARAMETER   :: mvbig = 1000, mbig = 600
REAL(dp), PARAMETER  :: elimit = -88.0_dp, sxteen = 16.0_dp, scale = 1.0D35, &
                      rootpi = 2.506628274631001_dp, hundrd = 100.0_dp

k = kk + 1
l = ll + 1
m = mm + 1
n = nn + 1
dir = .true.

!     Check arguments are within permitted limits

ifault = 1
fn_val = zero
IF (n < 1 .OR. m < n .OR. k < 1 .OR. k > m) RETURN

ifault = 2
IF (l < 1 .OR. k-l > m-n) RETURN
IF (.NOT. point) fn_val = one
IF (l > n .OR. l > k) RETURN
ifault = 0
fn_val = one
IF (k == 1 .OR. k == m .OR. n == 1 .OR. n == m) RETURN
IF (.NOT. point .AND. ll == MIN(kk, nn)) RETURN

p = REAL(nn, KIND=dp) / REAL(mm - nn, KIND=dp)
IF (REAL(MIN(kk, mm-kk)) > sxteen * MAX(p, one/p) .AND.  &
      mm > mvbig .AND. elimit > -hundrd) THEN
  
!     Use a normal approximation
  
  mean = REAL(kk) * REAL(nn) / REAL(mm)
  sig = SQRT(mean * (REAL(mm-nn) / REAL(mm)) * (REAL(mm-kk) / (REAL(mm-1))))
  IF (point) THEN
    arg = -half * (((REAL(ll) - mean) / sig)**2)
    fn_val = zero
    IF (arg >= elimit) fn_val = EXP(arg) / (sig * rootpi)
  ELSE
    fn_val = alnorm((REAL(ll) + half - mean) / sig, .false.)
  END IF
  
ELSE
  
!     Calculate exact hypergeometric probabilities.
!     Interchange K and N if this saves calculations.
  
  IF (MIN(k-1, m-k) > MIN(n-1, m-n)) THEN
    i = k
    k = n
    n = i
  END IF
  IF (m-k < k-1) THEN
    dir = .NOT. dir
    l = n - l + 1
    k = m - k + 1
  END IF
  
  IF (mm > mbig) THEN
    
!     Take logarithms of factorials.
    
    p = alnfac(nn) - alnfac(mm) + alnfac(mm-kk) + alnfac(kk) +  &
        alnfac(mm-nn) - alnfac(ll) - alnfac(nn-ll) - alnfac(kk-ll)  &
        - alnfac(mm-nn-kk+ll)
    fn_val = zero
    IF (p >= elimit) fn_val = EXP(p)
  ELSE
    
!     Use Freeman/Lund algorithm
    
    DO  i = 1, l-1
      fn_val = fn_val * REAL(k-i) * REAL(n-i) / (REAL(l-i) * REAL(:: m-i))
    END DO
    IF (l /= k) THEN
      j = m - n + l
      DO  i = l, k-1
        fn_val = fn_val * REAL(j-i) / REAL(m-i)
      END DO
    END IF
    
  END IF
  
  IF (point) RETURN
  IF (fn_val == zero) THEN
    
!     We must recompute the point probability since it has underflowed.
    
    IF (mm <= mbig) p = alnfac(nn) - alnfac(mm) + alnfac(kk) +  &
        alnfac(mm-nn) - alnfac(ll) - alnfac(nn-ll) - alnfac(kk-ll) -  &
        alnfac(mm-nn-kk+ll) + alnfac(mm-kk)
    p = p + LOG(scale)
    IF (p < elimit) THEN
      ifault = 3
      IF (ll > REAL(nn*kk + nn + kk +1)/(mm +2)) fn_val = one
      RETURN
    ELSE
      p = EXP(p)
    END IF
  ELSE
    
!     Scale up at this point.
    
    p = fn_val * scale
  END IF
  
  pt = zero
  nl = n - l
  kl = k - l
  mnkl = m - n - kl + 1
  IF (l <= kl) THEN
    DO  i = 1, l-1
      p = p * REAL(l-i) * REAL(mnkl-i) / (REAL(nl+i) * REAL(kl+i))
      pt = pt + p
    END DO
    IF (p == zero) ifault = 3
  ELSE
    dir = .NOT. dir
    DO  j = 0, kl-1
      p = p * REAL(nl-j) * REAL(kl-j) / (REAL(l+j) * REAL(mnkl+j))
      pt = pt + p
    END DO
    IF (p == zero) ifault = 3
  END IF
  
  IF (dir) THEN
    fn_val = fn_val + (pt / scale)
  ELSE
    fn_val = one - (pt / scale)
  END IF
  
END IF

END FUNCTION chyper