ClenshawCurtisNodesAndWeights Subroutine

SUBROUTINE ClenshawCurtisNodesAndWeights(N_in,xGP,wGP)
IMPLICIT NONE
!----------------------------------------------------------------------------------------------------------------------------------
! INPUT/OUTPUT VARIABLES
INTEGER,INTENT(IN)        :: N_in         !! polynomial degree, (N_in+1) CLpoints
REAL(wp),INTENT(OUT)          :: xGP(0:N_in)  !! Gauss point positions for the reference interval [-1,1]
REAL(wp),INTENT(OUT),OPTIONAL :: wGP(0:N_in)  !! Gauss point weights
!----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
INTEGER                   :: iGP,j
REAL(wp)                      :: b,theta
!==================================================================================================================================
IF(N_in.EQ.0)THEN
  xGP(0) = 0.0_wp
  wGP(0) = 2.0_wp
ELSE
  DO iGP=0,N_in
    xGP(iGP) = -COS(REAL(iGP,wp)*PP_Pi/REAL(N_in,wp))
  END DO
  xGP(0)   =-1.0_wp
  IF(MOD(N_in+1,2).EQ.1)THEN
    xGP(N_in/2)=0.0_wp
  END IF
  xGP(N_in)= 1.0_wp
  IF(PRESENT(wGP))THEN
    wGP=1.0_wp
    DO iGP=0,N_in
      theta = REAL(iGP,wp)*PP_Pi/REAL(N_in,wp)
      DO j=1,N_in/2
        b=MERGE(1.0_wp,2.0_wp,2*j.EQ.N_in)
        wGP(iGP) = wGP(iGP) - b * COS(REAL(2*j,wp)*theta) / REAL(4*j*j-1,wp)
      END DO
    END DO
    wGP(1:N_in-1)=2.0_wp*wGP(1:N_in-1)
    wGP=wGP/REAL(N_in,wp)
  END IF
END IF
END SUBROUTINE ClenshawCurtisNodesAndWeights