MODULE Global_data
  !Symbolic names for kind types of single- and double-precision reals: 
  INTEGER, PARAMETER :: SP = KIND(1.0) 
  INTEGER, PARAMETER :: DP = KIND(1.0D0) 
  !Frequently used mathematical constants (with precision to spare): 
  REAL(DP), PARAMETER :: PI=3.141592653589793238462643383279502884197_dp 

  ! parameters for the problem
  INTEGER :: n,m,x_print,y_print
  REAL(dp) :: dy,ymin,ymax,dx,xmin,xmax,omega,tol

  CHARACTER(LEN=50) :: results_file

CONTAINS

  SUBROUTINE input_parameters(filename)

      IMPLICIT NONE

      INTEGER,PARAMETER :: no_of_ints=4,no_of_reals=6,no_of_strings=1

      INTEGER :: the_ints(no_of_ints),i
      real(dp) :: the_reals(no_of_reals)
      CHARACTER(LEN=50) :: the_strings(no_of_strings)

      CHARACTER(LEN=*) :: filename

      OPEN(UNIT=9,file=filename)

      DO I = 1,no_of_ints
        READ(9,*)
         READ(9,*)the_ints(I)
      END DO
      DO i = 1,no_of_reals
         READ(9,*)
         READ(9,*)the_reals(i)
      ENDDO
      DO i = 1,no_of_strings
         READ(9,*)
         READ(9,*)the_strings(i)
      END DO

      n = the_ints(1) ; m = the_ints(2)
      x_print = the_ints(3) ; y_print = the_ints(4)
      xmin = the_reals(1) ; xmax = the_reals(2) ; ymin = the_reals(3)
      ymax = the_reals(4) ; omega = the_reals(5) ; tol = the_reals(6)
      results_file = the_strings(1)

      dX = (xmax - xmin)/dble(n)
      dy = (ymax - ymin)/dble(m)

      x_print = max(1,n/x_print)
      y_print = max(1,m/y_print)

      CLOSE(unit=9)

  END SUBROUTINE input_parameters

END MODULE Global_data

! This assigns the function " f(x,y) "

MODULE Poissin

  USE Global_data, wp=>dp

CONTAINS

  FUNCTION func(x,y)

     REAL(wp),INTENT(IN) :: x,y
     REAL(wp) :: func

     func = sin(2.*pi*x)*sin(2.*pi*y)/y/x

  END FUNCTION

END MODULE Poissin

! The problem we attempt to solve is:
!    d^2 U / dx^2 + d^2 U / dy^2 = f(x,y)
! We discretise the problem as
!
!   ( u_{i-1,j} - 2u_{i,j} + u_{i+1,j} )/(\Delta x)^2 + 
!         ( u_{i,j-1} - 2u_{i,j} + u_{i,j+1} )/(\Delta y)^2 = f_{i,j}
!
! so that we may write
!
! u_{i-1,j} + u_{i+1,j} - (2 + 2\beta^2)u_{i,j}
!     + \beta^2 u_{i,j-1} + \beta^2 u_{i,j+1}  = dx*dx*func(i,j)
!
!  Rearrange this equation to solve for u_{i,j}
!  Then solve using gauss-seidel point relaxation,
!
!     w = 1/(2 + 2\beta^2)*( ... )
!     u^{q+1}_{i,j} = (1-omega)u^{q}_{i,j} + omega w 

PROGRAM Poissins_equation
  
  USE global_data, wp=>dp
  IMPLICIT NONE
  
  INTEGER :: i,j

  REAL(wp) , ALLOCATABLE :: x(:),y(:),U(:,:)

  ! Read parameters from file

  CALL input_parameters('input.in')

  ! Allocate storage
  ALLOCATE(x(0:n),y(0:m),U(0:n,0:m))

  ! Assign initial values to x, y and U
  DO i = 0,n
        x(i) = xmin + dble(i)*dx
  ENDDO
  DO i = 0,m
        y(i) = ymin + dble(i)*dy
  ENDDO
  U = 0._wp

  WRITE(6,*)"Solve the Poissin equation over the domain"
  WRITE(6,'(e12.5,a,e12.5)')xmin,"< x < ",xmax
  WRITE(6,'(e12.5,a,e12.5)')ymin,"< y < ",ymax
  WRITE(6,*)"With ",n+1," nodes in the x direction and "
  WRITE(6,*)m+1," in the y direction."
  Write(6,'(a,e12.5,a,e12.5)')"Relaxation is at ",omega," and tolerance ",tol

  ! Solve using SOR scheme

  CALL Solve_with_SOR(x,y,U)

  ! Print results to file      

  OPEN(unit=10,file=TRIM(results_file)//'.dat')

  DO i = 0,n,x_print
     DO j=0,m,y_print
        WRITE(10,'(3(E12.5,1X))')X(i),Y(j),U(i,j)
     ENDDO
     WRITE(10,'(1X)')
  ENDDO

  CLOSE(UNIT=10)

END PROGRAM Poissins_equation

SUBROUTINE Solve_with_SOR(x,y,U)

  USE global_data, wp => dp
  USE poissin
  
  IMPLICIT NONE
  INTEGER :: i,j,counter,iter_max=100000

  REAL(wp),INTENT(IN) :: X(0:n),Y(0:m)
  REAL(wp),INTENT(INOUT) :: U(0:n,0:m)

  REAL(wp) :: w,error,beta,beta_sq,rho

  ! #############################################
  ! Set boundary conditions
  ! #############################################

  DO j = 0,m
     u(0,j) = 0._wp 
     u(n,j) = 0._wp
  ENDDO
  DO i = 1,n-1
     u(i,0) = 0._wp
     u(i,m) = 0._wp
  ENDDO

  ! ################################################
  !  Solve the scheme
  ! ################################################

  beta = dx/dy
  beta_sq = beta*beta

  DO counter = 1,iter_max


     ! set up error to contain the l2 norm
     !   error = sqrt(SUM_{i,j} (u^{q+1}-u^q)^2)
     error = 0._wp

     DO  i = 1,n-1

         ! interior points 

         DO j = 1,m-1

             ! calculate w from the Guass-Seidel formula, using u(), beta_sq, dx and func()
             w = ! your code here
             ! implement SOR using omega
             w = ! your code here
             error = error + (u(i,j)-w) * (u(i,j)-w) 
             u(i,j) = w

         ENDDO

     ENDDO

     ! Check for convergence on the L2 norm
     IF(sqrt(error) < tol)THEN
         WRITE(6,*)"Converged after ",counter,"iterations."
         exit
     ENDIF
  
  ENDDO

  IF(counter>iter_max)THEN
        WRITE(6,*)"Not converged with error ",error
        WRITE(6,*)"after ",iter_max," iterations."
        stop
  ENDIF

END SUBROUTINE Solve_with_SOR