program TestFactorial
 use iso_fortran_env, only : wp => real32
 implicit none
 ! integer, parameter :: wp=10 ! Try this too
 logical :: b0
 integer :: m0, m1
 write (*,'(a)') 'Testing the factorial function ...'
 call Factorial_Limits (m0, m1)
 ! m0 and m1 are the largest values such that using real (kind=wp):
 ! For 0<=k<m0 Factorial(k) can be represented exactly.
 ! For 0<=k<m1 Factorial(k) is below the overflow limit.
 write (*,'(a,2i5)') 'm0, m1:', m0, m1
 b0 = .false.
 ! b0 has the value (we have found a discrepancy between two ways of
 ! computing the factorial).
 block
 integer :: i
 real (kind=wp) :: t0, t1
 do i = 0, m0-1
 if (i==0) then
 t0 = 1
 else
 t0 = i*t0
 end if
 t1 = gamma(i+1.0_wp)
 ! t0 = i! computed using the recursion.
 ! t1 = i! computed via the gamma intrinsic function.
 if (t1/=t0) then
 b0 = .true.
 write (*,'(a,1x,i3,2(1x,1pg10.2))') 'error:', i, &
 spacing(t0), (t1-t0)/spacing(t0)
 end if
 end do
 end block
 if (.not.b0) then
 write (*,'(a)') '... PASS'
 end if
 stop
contains
 subroutine Factorial_Limits (m0, m1)
 integer, intent (out) :: m0, m1
 ! Determine the limits for computation of the Factorial function
 ! using real (kind=wp).
 ! m0: for 0<=k<m0, Factorial(k) can be represented exactly.
 ! m1: for 0<=k<m1, Factorial(k) is below the overflow limit.
 integer :: m0loc, i, j
 real (kind=wp) :: t0, t1
 logical :: b0
 i = 0 ; t0 = 1 ; t1 = 1 ; b0 = .true.
 m0loc = huge(0)
 ! Invariant: 0 <= i.
 ! t0 = Factorial(i) (subject to roundoff for large i).
 ! t1 = Factorial(i)/2^k where 2^k is the largest power of two that
 ! leaves an (odd) integer result.
 ! b0 = ((i+1)*t0 will not overflow).
 ! m0loc is the smallest nonnegative integer value i for which we
 ! have deterimined that Factorial(i) cannot be represented exactly.
 do while (b0)
 i = i+1
 t0 = i*t0
 j = i
 do while (modulo(j,2).eq.0)
 j = j/2
 end do
 t1 = j*t1
 if (1<spacing(t1)) then
 ! Factorial(i) cannot be represented exactly.
 m0loc = min(i,m0loc)
 end if
 b0 = log(t0)+log(real(i+1,kind=wp)) < log(huge(t0))
 end do
 m0 = min(i+1,m0loc)
 m1 = i+1
 return
 end subroutine Factorial_Limits
end program

bash-5.0$ gfortran TestFactorial.f90 
bash-5.0$ ./a.out
Testing the factorial function ...
m0, m1: 14 35
error: 2 2.38E-07 1.0 
error: 3 4.77E-07 1.0 
error: 5 7.63E-06 1.0 
error: 9 3.12E-02 1.0 
error: 11 4.0 1.0 
error: 12 32. 1.0 
error: 13 5.12E+02 2.0