      PROGRAM CL

      IMPLICIT NONE
      REAL*8 chisqr,CONF_LEVEL,conflevel
      INTEGER nfree

      WRITE(*,'(a,$)')'      nfree?  '
      READ(*,*)nfree
      WRITE(*,'(a,$)')' chi-square?  '
      READ(*,*)chisqr
      WRITE(*,*)
      IF(chisqr.LT.(REAL(nfree)/2.d0))THEN ! make a guess as to whether this 
        chisqr = chisqr*REAL(nfree) ! is a total chi^2 
        WRITE(*,*)'   * total '//
     &    'chi-square = ',chisqr
      ELSE                      ! or a reduced chi^2/dof
        WRITE(*,'(a,f6.3)')'  * reduced '//
     &    'chi-square =',chisqr/REAL(nfree)
      END IF

      conflevel = CONF_LEVEL(chisqr,nfree)
      IF(conflevel.GT.0.001d0)THEN
        WRITE(*,'(a,f8.4,a)')'    * '//
     &    'Confidence level = ',
     &    conflevel*100.d0,'%'
      ELSE
        WRITE(*,'(a,1pe11.4,a)')'    * '//
     &    'Confidence level = ',
     &    conflevel*100.d0,'%'
      END IF
      STOP
      END
***************************************************************************
      REAL*8 FUNCTION CONF_LEVEL(x,n)
***************************************************************************
      IMPLICIT NONE
      INTEGER n,m,i
      DOUBLE PRECISION x,srtopi,emyo2,sum,term,fi,srty,value,const,
     &     anu,an9,xp,z,DERF,FREQ
      LOGICAL first_call
      DATA first_call/.true./
 
      FREQ = 0.5d0*(1.d0+DERF(x/SQRT(2.d0)))
      srtopi = sqrt(2.d0/acos(-1.d0))
      CONF_LEVEL = 0.d0
      IF((n.LE.0).OR.(x.LT.0.d0))RETURN
      emyo2 = EXP((-0.5d0)*x)
      sum = 1.d0
      term = 1.d0
      m = INT(REAL(n)/2.d0)
      IF(x.LE.1000)THEN         ! non-asymptotic calculation
        IF(2*m.EQ.n)THEN        ! if "n" is even
          IF(m.NE.1)THEN
            DO i=2,m
              fi = REAL(i)-1.d0
              term = 0.5d0*term*x/fi
              sum = sum+term
            END DO
          END IF
          CONF_LEVEL = REAL(emyo2*sum)
        ELSE                    ! if "n" is odd
          srty = sqrt(x)
          FREQ = 0.5d0*(1.d0+DERF(srty/sqrt(2.d0)))
          value = 2.d0*(1.d0-FREQ)
          IF(n.EQ.1)THEN
            const = 0.d0
          ELSE
            const = srtopi*srty*emyo2
            IF(n.NE.3)THEN
              DO i=1,m-1
                fi = REAL(i)
                term = term*x/(2.d0*fi+1.d0)
                sum = sum+term
              END DO
            END IF
          END IF
          CONF_LEVEL = REAL(const*sum+value)
        END IF
      ELSE                      ! use asymptotic formula
        anu = 1.d0/REAL(n)
        an9 = anu/4.5d0
        xp  = 1.d0/3.d0
        z   = ((x*anu)**xp-(1.d0-an9))/sqrt(an9)
        FREQ = 0.5d0*(1.d0+DERF(z/sqrt(2.d0)))
        CONF_LEVEL = REAL(1.d0-FREQ)
      END IF
      RETURN
      END
