!     Last change:  PD   16 Apr 2001    0:42 am
MODULE ellip_solvers
!-----------------------------------------------------------------------
!  This module contains subroutines which solve Poissons equation using
!  2nd order central differences.  The methods employed are Jacobi,
!  Gauss-Seidel and SOR iterations.
!-----------------------------------------------------------------------                
 USE local
 IMPLICIT NONE
 INTEGER(i4), PARAMETER :: maxiter=10000
 INTEGER(i4) :: nx,ny,niter
 REAL(r8), DIMENSION(:,:), ALLOCATABLE :: w,r
 REAL(r8), DIMENSION(:), ALLOCATABLE :: x,y
 REAL(r8):: epsilon,dx,dy,error(maxiter),scale
 LOGICAL :: gauss                 

 CONTAINS

 SUBROUTINE jacobi()
!-----------------------------------------------------------------------
!  This subroutine uses the Jacobi method of iteration which uses only
!  previous grid points to update the next value of w.
!  w(i,j) = .25*(v(i+1,j) + v(i-1,j) + v(i,j-1) + v(i,j+1) - r(i,j)*(dx**2)).  
!  The array r represents the RHS of Poisson's equation.  
!-----------------------------------------------------------------------              
 USE local
 IMPLICIT NONE
 INTEGER(i4) :: i,j
 REAL(r8), DIMENSION(:,:), ALLOCATABLE :: v

  ALLOCATE(v(0:nx,0:ny))
  niter = 1
  error = 0.
  DO 
!-----------------------------------------------------------------------
!  Update each of the grid points for the next iteration.
!-----------------------------------------------------------------------
   v = w
   DO j=1,ny-1
    DO i=1,nx-1
     w(i,j) = .5*(scale*(v(i+1,j) + v(i-1,j)) + v(i,j-1) + v(i,j+1) - r(i,j)*dy**2)/(1.+scale)
     error(niter) = error(niter) + ABS(w(i,j) - v(i,j))
    END DO
   END DO
   IF((abs(error(niter)).lt.epsilon).or.(niter.ge.maxiter)) EXIT
   niter = niter+1
  END DO
  DEALLOCATE(v)
 RETURN
END SUBROUTINE jacobi

SUBROUTINE gauss_sor(gauss,omega)
!-----------------------------------------------------------------------
!  This subroutine uses either the Gauss-Seidel method of iteration which entails
!  using the most recent value of w for each grid point in the calculation.
!  w(i,j) = .25*(w(i+1,j) + w(i-1,j) + w(i,j-1) +w(i,j+1) - r(i,j)*(dx**2)), or
!  the SOR method of iteration which uses a the most recent values of w and 
!  the slope of the calculation to better estimate the value of w(k+1).  
!  w(i,j) = w(i,j) + omega*temp where temp = .25*(w(i+1,j) + w(i-1,j) + w(i,j-1) 
!  + w(i,j+1) - r(i,j)*(dx**2)) - w(i,j).  The value of omega can be adjusted to 
!  optimize the calculation.     
!-----------------------------------------------------------------------
 USE local
 IMPLICIT NONE
 INTEGER(i4) :: i,j
 REAL(r8), INTENT(IN) :: omega
 LOGICAL, INTENT(IN) :: gauss
 REAL(r8) :: temp
!-----------------------------------------------------------------------
!  Initialize the general variables needed for the convergence calculation.
!-----------------------------------------------------------------------
  niter = 1
  error = 0.0
  DO
   IF((abs(error(niter)).lt.epsilon).or.(niter.ge.maxiter)) EXIT
   niter = niter+1                     
  END DO
 RETURN
END SUBROUTINE gauss_sor
END MODULE ellip_solvers