!     Last change:  PD   16 Apr 2001    0:35 am
PROGRAM elliptic 
 USE local
 USE ellip_solvers
 IMPLICIT NONE
 INTEGER(i4) :: i,j,method
 REAL(r8) :: ymin,ymax,xmin,xmax,init,omega
 REAL(r8) :: len_rat,h_rat

  epsilon = 1.e-6
  nx = 40
  ny = 40
  len_rat = 1.8
  xmin = -1.0
  xmax = 1.0
  ymin = -1.0
  ymax = 1.0
  dx = (xmax-xmin)/nx
  dy = (ymax-ymin)/ny
  h_rat = dy/dx
  scale = (len_rat*h_rat)**2
  ALLOCATE(w(0:nx,0:ny),r(0:nx,0:ny),x(0:nx),y(0:ny))
  DO i=0,nx
   x(i) = xmin + i*dx
  END DO
  DO i=0,ny
   y(i) = ymin + i*dy
  END DO        
!-----------------------------------------------------------------------
!  Initialize the solution vector and the matrix which represents the RHS of
!  the dimensionless Poisson's equation.
!-----------------------------------------------------------------------
  WRITE(*,*) 'Enter the initial value for the Dirichlet conditions:'
  READ(*,*) init
  w = init
  r = -1.
!-----------------------------------------------------------------------
!  Initialize the Dirichlet boundary conditions.
!-----------------------------------------------------------------------
   w(0,:) = 0.0
   w(nx,:) = 0.0
   w(:,0) = 0.0
   w(:,ny) = 0.0
!-----------------------------------------------------------------------
!  Call the desired subroutine. 
!-----------------------------------------------------------------------
  WRITE(*,*) 'Enter the desired iterative method: '
  WRITE(*,*) '(1) Jacobi,  (2) Gauss-Seidel,  (3) SOR'
  READ(*,*) method

  IF (method.eq.1) THEN
   CALL jacobi()
  ELSE 
   IF(method.eq.2) THEN
    gauss = .true.
   ELSE
    WRITE(*,*) 'Enter a value for omega'
    READ(*,*) omega
    gauss = .false.
   END IF
   CALL gauss_sor(gauss,omega)
  END IF
!-----------------------------------------------------------------------
!  Upon completion of respective method write out the contours of w and
!  the values of w at (i=5,j=5), (i=20,j=5), and (i=20,j=20).
!-----------------------------------------------------------------------
  OPEN(unit=88,file='contour.dat',status='unknown')
  OPEN(unit=99,file='error.dat',status='unknown')
  DO j=0,ny
   DO i=0,nx
    WRITE(88,2) x(i),1.8*y(j),w(i,j)
   END DO
  END DO
  WRITE(*,*) 'The number of iterations was:  ', niter
  DO i=1,niter
   WRITE(99,1) i, error(i)
  END DO
  CLOSE(88)
  CLOSE(99)

  WRITE(*,*) '       w(5,5)        w(20,5)         w(20,20)'
  WRITE(*,2) w(5,5), w(20,5), w(20,20)
1 FORMAT(' ',i6,f16.8)
2 FORMAT(' ',3(f16.8))
  DEALLOCATE(w,r)
 STOP
END PROGRAM elliptic