ChebyshevGaussLobatto Subroutine

  subroutine ChebyshevGaussLobatto(N,nodes,weights)
    implicit none
    integer       :: N
    real(real64)    :: nodes(0:N)
    real(real64)    :: weights(0:N)
    ! LOCAL
    integer    :: j

    do j = 0,N

      weights(j) = pi/real(N,real64)
      nodes(j) = -cos(pi*real(j,real64)/real(N,real64))

    enddo

    weights(0) = weights(0)*0.5_real64
    weights(N) = weights(N)*0.5_real64

  endsubroutine ChebyshevGaussLobatto