genere_coeff_lagrange Subroutine

subroutine genere_coeff_lagrange()
!! subroutine generating coefficients for kth-order interpolation
implicit none

   integer(kind=I4) :: i, j, k, n, c
   real(kind=R8)    :: coeff

   do

      write(*,*) 'n' ; read(*,*) n ; if (n==0) exit
      write(*,*) 'k' ; read(*,*) k
      write(*,*) 'c' ; read(*,*) c

      do i = 0, n

         coeff = 1.0d0
         do j = 0, n
            if (j==i) cycle
            coeff = coeff * ( 0.5_R8 * ( 2 * k - 2 * j + 1 ) )/( i - j )
         enddo
         write(*,*) coeff*c

      enddo

   enddo

return
endsubroutine genere_coeff_lagrange