Merge pull request #1557 from ExtremeFLOW/fix/legendre_poly

Fix legendre_poly

src/sem/speclib.f90

@@ -973,18 +973,21 @@ REAL(KIND=XP) FUNCTION PNLEG (Z,N)
     if (n.eq.0) pnleg = 1.
     RETURN
   end function PNLEG
-
+
+  !> Evaluate Legendre polynomials of degrees 0-N at point x
+  !! and store in array L
   subroutine legendre_poly(L, x, N)
-    ! Evaluate Legendre polynomials of degrees 0-N at point x
-    real(kind=rp), intent(inout):: L(1:N+1)
+    real(kind=rp), intent(inout):: L(0:N)
     real(kind=rp) :: x
     integer :: N, j
 
-    L(1) = 1.0_rp
-    L(2) = x
+    L(0) = 1.0_xp
+    if (N .eq. 0) return
+    L(1) = x
 
-    do j=3, N+1
-       L(j) = ( (2*j-1) * x * L(j-1) - (j-1) * L(j-2) ) / j 
+    do j=1, N-1
+       L(j+1) = ( (2.0_xp * real(j, xp) + 1.0_xp) * x * L(j) &
+            - real(j, xp) * L(j-1) ) / (real(j, xp) + 1.0_xp) 
     end do
   end subroutine legendre_poly