PROGRAM test_chi_squared
IMPLICIT NONE
INTEGER, PARAMETER :: doub_prec = SELECTED_REAL_KIND(15, 60)
INTEGER            :: ndf
REAL (doub_prec)   :: chi2

INTERFACE
  FUNCTION chi_squared(ndf, chi2) RESULT(prob)
    IMPLICIT NONE
    INTEGER, PARAMETER           :: doub_prec = SELECTED_REAL_KIND(15, 60)
    INTEGER, INTENT(IN)          :: ndf
    REAL (doub_prec), INTENT(IN) :: chi2
    REAL (doub_prec)             :: prob
  END FUNCTION chi_squared
END INTERFACE

DO
  WRITE(*, *) 'Enter no. of degrees of freedom: '
  READ(*, *) ndf
  WRITE(*, *) 'Enter chi-squared value: '
  READ(*, *) chi2
  WRITE(*, *) 'Chi-squared probability = ', chi_squared(ndf, chi2)
END DO
STOP

END PROGRAM test_chi_squared


FUNCTION chi_squared(ndf, chi2) RESULT(prob)

IMPLICIT NONE
INTEGER, PARAMETER           :: doub_prec = SELECTED_REAL_KIND(15, 60)
INTEGER, INTENT(IN)          :: ndf
REAL (doub_prec), INTENT(IN) :: chi2
REAL (doub_prec)             :: prob

INTERFACE
  FUNCTION gammad(x, p) RESULT(gamma_prob)
    IMPLICIT NONE
    INTEGER, PARAMETER           :: doub_prec = SELECTED_REAL_KIND(15, 60)
    REAL (doub_prec), INTENT(IN) :: x, p
    REAL (doub_prec)             :: gamma_prob
  END FUNCTION gammad
END INTERFACE

REAL (doub_prec) :: half = 0.5d0, x, p

x = half * chi2
p = half * REAL(ndf)
prob = gammad(x, p)
RETURN

END FUNCTION chi_squared