C     Elementary multi-grid solution of the Poisson equation in a
C     rectangle.  Adapted from a program by W. O'Farrell.

      program findiff

      real  x(128,128),  xold(128,128),  charge(128,128)
      real x2(256,256), x2old(256,256), charge2(256,256)
      real epsilon
      integer tm, time1, time2, i, j

      x = 0.0
      x(:,1) = 100.0
      x(1,:) = 100.0
      charge = 0.1
      i = 0
      time1 = 0
      epsilon = 1
      call mpTimerStart()
      do while (epsilon .GT. .02)
         xold = x
         x(2:127,2:127) = ((x(1:126,2:127)+x(3:128,2:127)
     c      +x(2:127,1:126)+x(2:127,3:128))*.25)+charge(2:127,2:127)
         i = i + 1
         epsilon = maxval(abs(x-xold))
         if (mod(i,1000) .eq. 0) then
            tm = mpTimerElapsed()
            time1 = time1 + tm
            print *,i,epsilon, tm
            call mpTimerStart
         end if         
      end do
      print *,i, epsilon, time1

      print *, ' Refining grid'
      charge2 = 0.025           ! Reflects finer spacing
      x2(::2,::2) = x           ! Interpolate interior
      x2(2::2,:) = x2(::2,:)    ! by copying
      x2(:,2::2) = x2(:,::2)
      x2(:,1) = 1.0            ! Install boundary values
      x2(1,:) = 1.0
      x2(:,256) = 0.0
      x2(256,:) = 0.0
      j = 0
      time2 = 0
      epsilon = 1
      call mpTimerStart()
      do while (epsilon .GT. .002)
         x2old = x2
         x2(2:255,2:255) = ((x2(1:254,2:255)+x2(3:256,2:255)
     c      +x2(2:255,1:254)+x2(2:255,3:256))*.25)+charge2(2:255,2:255)
         j = j + 1
         epsilon = maxval(abs(x2-x2old))
         if (mod(j,1000) .eq. 0) then
            tm = mpTimerElapsed()
            time2 = time2 + tm
            print *,j,epsilon,tm
            call mpTimerStart()
         end if         
      end do
      print *,j, epsilon, time2

      print *
      print *, ' Coarse solution (63:66, 63:66):'
      print *, x(63:66,63:66)
      print *
      print *, ' Fine solution, (127:130, 127:130):'
      print *, x2(127:130,127:130)
      print *
      print *,'iterations = ', i, j
      stop
      end

