bernob Subroutine

      subroutine bernob(n,Bn)

      integer,intent(in) :: n !! Serial number
      real(wp),intent(out) :: Bn(0:n) !! `Bn`

      real(wp) :: r1 , r2 , s
      integer :: k , m

      real(wp),parameter :: tol = 1.0e-15_wp
      integer,parameter :: maxiter = 10000

      Bn(0) = 1.0_wp
      Bn(1) = -0.5_wp
      Bn(2) = 1.0_wp/6.0_wp
      r1 = (2.0_wp/twopi)**2
      do m = 4 , n , 2
         r1 = -r1*(m-1)*m/(twopi*twopi)
         r2 = 1.0_wp
         do k = 2 , maxiter
            s = (1.0_wp/k)**m
            r2 = r2 + s
            if ( s<tol ) exit
         enddo
         Bn(m) = r1*r2
      enddo

   end subroutine bernob