PROGRAM TestMachEps

        IMPLICIT NONE
        INTEGER,PARAMETER :: dp = 8
        REAL (kind=dp)    :: eps,MachEps

        eps = MachEps()

        Write(*,'(//,1PD20.12,//)') eps

      END PROGRAM TestMachEps

      FUNCTION MachEps()

        IMPLICIT NONE
        INTEGER,PARAMETER :: dp = 8
        REAL (kind=dp)    :: MachEps
        REAL (kind=dp)    :: eps,EinsPlusEps
     
        eps=0.5 
        EinsPlusEps=1.5
        DO WHILE (EinsPlusEps .NE. 1.0D0)
          eps=0.5D0*eps
          EinsPlusEps=1.0D0+eps
        ENDDO
        MachEps=2.0D0*eps
        RETURN
      END FUNCTION MachEps