program poisson

   implicit none

   integer linear_dim_grid

   ! Dimension of coefficient matrix
   integer dim_A

   ! Auxilliary variables
   integer i, j, k, k_nn

    real(kind=8), dimension(:,:), allocatable :: grid
    real(kind=8), dimension(:,:), allocatable :: coefficient_matrix

    
   linear_dim_grid = 5

    dim_A = linear_dim_grid * linear_dim_grid


    allocate(grid(linear_dim_grid, linear_dim_grid))
    allocate(coefficient_matrix(dim_A, dim_A))

    ! Initialise vanishing matrix elements
    grid = 0.0d0
    coefficient_matrix = 0.0d0
     
    do i = 1, linear_dim_grid
       do j = 1, linear_dim_grid
          ! Current index in basis of LES
          k = (i-1)*linear_dim_grid + (j-1) + 1
          ! Rechter Nachbar
          if ((j+1).le.linear_dim_grid) then
             k_nn = (i-1)*linear_dim_grid + (j+1-1) + 1
             coefficient_matrix(k,k_nn) = 1.0d0
          endif
          ! Linker Nachbar
          if ((j-1).gt.0) then
             k_nn = (i-1)*linear_dim_grid + (j-1-1) + 1
             coefficient_matrix(k,k_nn) = 1.0d0
          endif
          ! Unterer Nachbar
          if ((i+1).le.linear_dim_grid) then
             k_nn = (i+1-1)*linear_dim_grid + (j-1) + 1
             coefficient_matrix(k,k_nn) = 1.0d0
          endif
          ! Oberer Nachbar
          if ((i-1).gt.0) then
             k_nn = (i-1-1)*linear_dim_grid + (j-1) + 1
             coefficient_matrix(k,k_nn) = 1.0d0
          endif
          ! Mich selbst
          coefficient_matrix(k,k) = -4.0d0
      enddo
   enddo

   do i = 1, dim_A
      write(*,'(25f5.0)') (coefficient_matrix(i,j),j=1, dim_A)
   enddo

end program
