[数值问题] 运行时出现 Program Exception - access violation，不知道什么原因？

module gauss
  !高斯消元模块：主要为牛顿法中解线性方程所调用
  contains
  subroutine lineq(A,b,x,N)
    implicit none
    integer(kind=2) :: id_max,i,k,N!主元素标号
    real(kind=8) :: A(N,N),b(N),x(N)
    real(kind=8) :: Aup(N,N),bup(N)
    real(kind=8) :: Ab(N,N+1)  !增广矩阵
    real(kind=8) :: vtemp1(N+1),vtemp2(N+1)
    real(kind=8) :: elmax,temp
    Ab(1:N,1:N) = A
    Ab(:,N+1) = b
    do k=1,N-1
      elmax = dabs(Ab(k,k))
      id_max = k
      do i=k+1,N
        if(abs(Ab(i,k))>elmax) then
          elmax = Ab(i,k)
          id_max = i
        endif
      enddo
      vtemp1 = Ab(k,:)
      vtemp2 = Ab(id_max,:)
      Ab(k,:) = vtemp2
      Ab(id_max,:) = vtemp1
      do i=k+1,N
        temp = Ab(i,k)/Ab(k,k)
        Ab(i,:) = Ab(i,:) - temp*Ab(k,:)
      enddo
    enddo
    Aup(:,:) = Ab(1:N,1:N)
    bup(:) = Ab(:,N+1)
    call uptri(Aup,bup,x,N)
    return
  end subroutine lineq
  
  !上三角线性方程组
  subroutine uptri(A,b,x,N)
    implicit none
    integer(kind=2) :: i,j,N
    real(kind=8) :: A(N,N),b(N),x(N)
    x(N) = b(N)/A(N,N)
    do i=n-1,1,-1
      x(i) = b(i)
      do j=i+1,N
        x(i) = x(i) - A(i,j)*x(j)
      enddo
      x(i) = x(i)/A(i,i)
    enddo
    return 
   end subroutine uptri
end module gauss
  
module newton
  contains
  subroutine solve
    use gauss
    implicit  none
    integer(kind=2) :: i,itmax = 500!最大允许迭代次数
    integer(kind=2) :: p,q,N=2
    real(kind=8) :: x(2),f(2),dx(2)
    real(kind=8) :: x0(2),l
    real(kind=8),allocatable :: s(:),t(:),s0(:),t0(:)
    real(kind=8) :: df(2,2) !偏导数矩阵
    real(kind=8) :: tol,dx2
    q=(p+1)*(p+2)/2
    allocate(s(q))
    allocate(t(q))
    allocate(s0(q))
    allocate(t0(q))
    !随机给定初值
    call RANDOM_SEED()
    call RANDOM_NUMBER(x)
    tol = 1d-4  !计算精度控制
    open(unit=11,file = 'parm.txt')
    read(11,*) p,x0,l
    write(*,*) p,x0,l
    do i=1,q
      read(11,100) s0(i),t0(i) 
      write(*,100) s0(i),t0(i) 
      pause '=============================='
    enddo
    !x = x0 
    s = s0/100
    t = t0/100
100 format(E15.8E2,1x,E15.8E2)
    open(unit=11,file = 'result.txt')
    write(11,101)
101 format(/,t6,'牛顿法计算非线性方程组迭代序列',/) 
    do i=1,itmax 
      call func(f,p,x,s,t,x0,l)
      call jac(df,p,x,s,t,x0,l)
      call lineq(df,-f,dx,N)
      x = x + dx
      write(*,102) i , x*l+x0
102   format(i5,2f16.10)      
      dx2 = dsqrt(dx(1)**2+dx(2)**2)
      if (dx2 < tol) exit
    enddo
    return
  end subroutine solve
  
  subroutine func(f,p,x,s,t,x0,l)
    implicit none
    real(kind=8) :: x(2),f(2)
    real(kind=8) :: x0(2),l
    integer(kind=2) :: p
    integer :: i,j,k=0
    real(kind=8) :: s((p+1)*(p+2)/2),t((p+1)*(p+2)/2)
    f = 0
    !f(1) = s(1)+x0(1)-10+(s(2)+l)*x(1)+s(3)*x(2)+s(4)*x(1)**2+s(5)*x(1)*x(2)+s(6)*x(2)**2+&
    !&s(7)*x(1)**3+s(8)*x(1)**2*x(2)+s(9)*x(1)*x(2)**2+s(10)*x(2)**3!+&
    !&s(11)*x(1)**4+s(12)*x(1)**3*x(2)+s(13)*x(1)**2*x(2)**2+s(14)*x(1)*x(2)**3+s(15)*x(2)**4
    !f(2) = t(1)+x0(2)-4.8+t(2)*x(1)+(t(3)+l)*x(2)+t(4)*x(1)**2+t(5)*x(1)*x(2)+t(6)*x(2)**2+&
    !&t(7)*x(1)**3+t(8)*x(1)**2*x(2)+t(9)*x(1)*x(2)**2+t(10)*x(2)**3!+&
    !&t(11)*x(1)**4+t(12)*x(1)**3*x(2)+t(13)*x(1)**2*x(2)**2+t(14)*x(1)*x(2)**3+t(15)*x(2)**4
    !s(2)=s(2)+l
    !t(3)=t(3)+l
    !s(1)=s(1)+x0(1)-10
    !t(1)=t(1)+x0(1)-4.8
    do i=0,p
      do j=0,i
        k = k+1
        f(1)=f(1)+s(k)*x(1)**(i-j)*x(2)**j
      enddo
    enddo
    k=0    
    do i=0,p
      do j=0,i
        k = k+1
        f(2)=f(2)+t(k)*x(1)**(i-j)*x(2)**j
      enddo
    enddo
  end subroutine func
  subroutine jac(df,p,x,s,t,x0,l)
    implicit none
    real(kind=8) :: x(2),df(2,2)
    integer(kind=2) :: p,i,j,k,m
    real(kind=8) :: s((p+1)*(p+2)/2),t((p+1)*(p+2)/2)
    real(kind=8) :: x0(2),l
    !df(1,1) = s(2)+l+2*s(4)*x(1)+s(5)*x(2)+3*s(7)*x(1)**2+2*s(8)*x(1)*x(2)+s(9)*x(2)**2!+&
              !&4*s(11)*x(1)**3+3*s(12)*x(1)**2*x(2)+2*s(13)*x(1)*x(2)**2+s(14)*x(2)**3
    !df(2,1) = s(3)+s(5)*x(1)+2*s(6)*x(2)+s(8)*x(1)**2+2*s(9)*x(1)*x(2)+3*s(10)*x(2)**3!+&
              !&s(12)*x(1)**3+2*s(13)*x(1)**2*x(2)+3*s(14)*x(1)*x(2)**2+4*s(15)*x(2)**3
    !df(1,2) = t(2)+2*t(4)*x(1)+t(5)*x(2)+3*t(7)*x(1)**2+2*t(8)*x(1)*x(2)+t(9)*x(2)**2!+&
              !&4*t(11)*x(1)**3+3*t(12)*x(1)**2*x(2)+2*t(13)*x(1)*x(2)**2+t(14)*x(2)**3
    !df(2,2) = t(3)+l+t(5)*x(1)+2*t(6)*x(2)+t(8)*x(1)**2+2*t(9)*x(1)*x(2)+3*t(10)*x(2)**3!+&
              !&s(12)*x(1)**3+2*s(13)*x(1)**2*x(2)+3*s(14)*x(1)*x(2)**2+4*s(15)*x(2)**3
    !s(2)=s(2)+l
    !t(3)=t(3)+l           
    df=0         
    !df(1,1)
    m=0
    do i=0,p-1
      k = i+1
      m = m+1
      do j=0,i
        m=m+1
        df(1,1)=df(1,1)+k*s(m)*x(1)**(i-j)*x(2)*j
        k=k-1 
      enddo
    enddo
    !df(1,2)
    m=0
    do i=0,p-1
      k = i+1
      m = m+1
      do j=0,i
        m=m+1
        df(1,2)=df(1,2)+k*t(m)*x(1)**(i-j)*x(2)*j
        k=k-1
      enddo
    enddo
    !df(2,1)
    m=1
    do i=0,p-1
      k = 1
      m = m+1
      do j=0,i
        m=m+1        
        df(2,1)=df(2,1)+k*s(m)*x(1)**(i-j)*x(2)*j
        k=k+1
      enddo
    enddo
    !df(2,2)
    m=1
    do i=0,p-1
      k = 1
      m = m+1
      do j=0,i
        m=m+1        
        df(2,2)=df(2,2)+k*t(m)*x(1)**(i-j)*x(2)*j
        k=k+1
      enddo
    enddo
  end subroutine jac
end module newton
  
  
program main
  use newton
  call solve
end program main