LegendreGaussLobatto Subroutine

  subroutine LegendreGaussLobatto(N,nodes,weights)
    implicit none
    integer    :: N
    real(real64) :: nodes(0:N)
    real(real64) :: weights(0:N)
    ! Local
    real(real64) :: nodes_local(0:N)
    real(real64) :: weights_local(0:N)
    real(real64) :: delta,q,qprime,lN
    integer  :: j,kIt

    if(N == 1) then

      nodes_local(0) = -1.0_real64
      weights_local(0) = 1.0_real64
      nodes_local(1) = 1.0_real64
      weights_local(1) = 1.0_real64

    else

      nodes_local(0) = -1.0_real64
      weights_local(0) = 2.0_real64/(real(N,real64)*(real(N,real64)+1.0_real64))
      nodes_local(N) = 1.0_real64
      weights_local(N) = weights_local(0)

      do j = 1,((N+1)/2-1)

        nodes_local(j) = -cos((real(j,real64)+0.25_real64)*pi/real(N,real64)- &
                              3.0_real64/(8.0_real64*real(N,real64)*pi*(real(j,real64)+0.25_real64)))

        do kIt = 1,newtonMax

          call LegendreQandL(N,nodes_local(j),q,qprime,lN)

          delta = -q/qprime
          nodes_local(j) = nodes_local(j)+delta
          if(abs(delta) <= TOL*nodes_local(j)) exit

        enddo

        call LegendreQandL(N,nodes_local(j),q,qprime,lN)

        weights_local(j) = 2.0_real64/(real(N,real64)*(real(N,real64)+1.0_real64)*lN*lN)
        weights_local(N-j) = weights_local(j)
        nodes_local(N-j) = -nodes_local(j)

      enddo

    endif

    if(mod(real(N,real64),2.0_real64) == 0.0_real64) then

      call LegendreQandL(N,0.0_real64,q,qprime,lN)

      nodes_local(N/2) = 0.0_real64
      weights_local(N/2) = 2.0_real64/(real(N,real64)*(real(N,real64)+1.0_real64)*lN*lN)

    endif

    do j = 0,N
      nodes(j) = real(nodes_local(j),real64)
      weights(j) = real(weights_local(j),real64)
    enddo

  endsubroutine LegendreGaussLobatto