src/math.f90: use the intrinsic gamma function

- remove the `gamm` function
- rewrite `fact` in terms of `gamma`

src/dispersion.f90

@@ -904,7 +904,6 @@ end subroutine antihermitian
 !
 !
 subroutine ssbi(zz,n,l,fsbi)
-  use math, only : gamm
   implicit none
 ! local constants
   integer, parameter :: lmx=20,nmx=lmx+2
@@ -918,7 +917,7 @@ subroutine ssbi(zz,n,l,fsbi)
   real(wp_) :: c0,c1,sbi
 !
   do m=n,l+2
-    c0=one/gamm(dble(m)+1.5_wp_)
+    c0=one/gamma(dble(m)+1.5_wp_)
     sbi=c0
     do k=1,50
       c1=c0*0.25_wp_*zz**2/(dble(m+k)+0.5_wp_)/dble(k)

src/math.f90

@@ -84,42 +84,17 @@ contains
     end if
   end function catand
 
+
   function fact(k)
+    ! Factorial function k!
+    ! Note: the result is real
+
     implicit none
     integer, intent(in) :: k
-    real(wp_) :: fact
-    integer :: i
-!   Factorial function
-    fact=zero
-    if(k<0) return
-    fact=one
-    if(k==0) return
-    do i=1,k
-      fact=fact*i
-    end do
-  end function fact
+    real(wp_)           :: fact
 
-  function gamm(xx)
-    implicit none
-    real(wp_) :: gamm
-    real(wp_), intent(in) :: xx
-!   Returns the value Gamma(xx) for xx > 0.
-    INTEGER :: j
-    real(wp_) :: ser,tmp,x,y
-    real(wp_), parameter :: stp=2.5066282746310005_wp_
-    real(wp_), dimension(6), parameter :: cof=(/76.18009172947146_wp_, &
-    -86.50532032941677_wp_,24.01409824083091_wp_,-1.231739572450155_wp_, &
-    .1208650973866179e-2_wp_,-.5395239384953e-5_wp_/)
-    x=xx
-    y=x
-    tmp=x+5.5_wp_
-    tmp=(x+0.5_wp_)*log(tmp)-tmp
-    ser=1.000000000190015_wp_
-    do  j=1,6
-      y=y+1._wp_
-      ser=ser+cof(j)/y
-    end do
-    gamm=exp(tmp)*(stp*ser/x)
-  end function gamm
+    fact = gamma(real(1 + k, kind=wp_))
 
-end module math
+   end function fact
+
+end module math