recursive  function fact(n) result (fact1)

    implicit none
    real(qp) ::  fact1
    real(qp), intent(in) :: n

    if (n.eq.0.or.n.eq.1) then
       fact1 = 1
       return
    else
       fact1 = n*fact(n-1)
      return
    end if

  end function fact