LegendrePolynomialAndDerivative Subroutine

SUBROUTINE LegendrePolynomialAndDerivative(N_in,x,L,Lder)
IMPLICIT NONE
!----------------------------------------------------------------------------------------------------------------------------------
! INPUT/OUTPUT VARIABLES
INTEGER,INTENT(IN) :: N_in   !! (IN)  polynomial degree, (N+1) CLpoints
REAL(wp),INTENT(IN)    :: x      !! (IN)  coordinate value in the interval [-1,1]
REAL(wp),INTENT(OUT)   :: L      !! (OUT) Legedre polynomial evaluated at \f$ \xi: L_N(\xi), \partial/\partial\xi L_N(\xi) \f$
REAL(wp),INTENT(OUT)   :: Lder   !! (OUT) Legedre polynomial deriv. evaluated at \f$ \xi: L_N(\xi), \partial/\partial\xi L_N(\xi) \f$
!----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
INTEGER :: iLegendre
REAL(wp)    :: L_Nm1,L_Nm2 ! L_{N_in-2},L_{N_in-1}
REAL(wp)    :: Lder_Nm1,Lder_Nm2 ! Lder_{N_in-2},Lder_{N_in-1}
!==================================================================================================================================
IF(N_in .EQ. 0)THEN
  L=1.0_wp
  Lder=0.0_wp
ELSEIF(N_in .EQ. 1) THEN
  L=x
  Lder=1.0_wp
ELSE ! N_in > 1
  L_Nm2=1.0_wp
  L_Nm1=x
  Lder_Nm2=0.0_wp
  Lder_Nm1=1.0_wp
  DO iLegendre=2,N_in
    L=(REAL(2*iLegendre-1,wp)*x*L_Nm1 - REAL(iLegendre-1,wp)*L_Nm2)/REAL(iLegendre,wp)
    Lder=Lder_Nm2 + REAL(2*iLegendre-1,wp)*L_Nm1
    L_Nm2=L_Nm1
    L_Nm1=L
    Lder_Nm2=Lder_Nm1
    Lder_Nm1=Lder
  END DO !iLegendre=2,N_in
END IF ! N_in
!normalize
L=L*SQRT(REAL(N_in,wp)+0.5_wp)
Lder=Lder*SQRT(REAL(N_in,wp)+0.5_wp)
END SUBROUTINE LegendrePolynomialAndDerivative