[求助] 求找出问题，谢谢

program mytp12 
  implicit none 
        integer:: n
         double precision,allocatable::xk(:),wk(:)
        write(6,*)'input order of quadrature n' 
        read(5,*)n 
    allocate(xk(n),wk(n))
    call GL_xw(n,xk,wk)

 end program mytp12

        double precision function leg(n,x,dp) 
        implicit none 
            integer :: k,n 
            double precision:: dp,p(0:n),x 
         p(0) = 1.d0 
         p(1) = x 
         do k = 2, n 
            p(k)=((2.d0*k-1.d0)*x*p(k-1)-(k-1.d0)*p(k-2))/k 
         end do 
         leg=p(n) 
         dp=(n*p(n-1)-n*x*p(n))/(1-x**2) 
        end function leg 

    subroutine GL_xw(n,xk,wk)
    implicit none
    integer n,k 
         double precision xk(n),wk(n) 
        double precision, parameter :: eps = 1.d-15 
        double precision dp,dx,leg,p,pi,x 
        pi = acos(-1.d0)  
         do k = 1, n 
           x=cos((4.d0*k-1.d0)*pi/(4.d0*n+2.d0)) 
            do 
              p = leg(x,n,dp) 
              dx = p / dp 
              x = x - dx 
              if (abs (dx) < eps) exit 
             end do 
           xk(k) = x 
         wk(k)=2.d0/((1.d0-x**2)*dp*dp) 
        end do 
        do k = 1, n 
          write (6,'("xk(",i3,")= ",d15.8," wk(",i3,")= ",d15.8)') k,xk(k),k,wk(k) 
        end do 
    return
        endsubroutine
!  --------------------------------------------------
!  Silverfrost FTN95 for Microsoft Visual Studio
!  Free Format FTN95 Source File
!  --------------------------------------------------