program laplace
      integer imax,jmax,im1,im2,jm1,jm2,it,itmax
      ! parameter (imax=2001,jmax=2001)
      integer, parameter :: imax=501
      integer, parameter :: jmax=501
      integer, parameter :: im1=imax-1
      integer, parameter :: im2=imax-2
      integer, parameter :: jm1=jmax-1
      integer, parameter :: jm2=jmax-2
      integer, parameter :: itmax=100
      real*8 u(imax,jmax),du(imax,jmax),umax,dumax,tol
      ! parameter (umax=10.0,tol=1.0e-3)
      parameter (umax=10.0, tol=1.0e-2)
      real t(2), e
      integer clock_count1, clock_count2, clock_max, clock_rate

      ! Initialize
      do j=1,jmax
         do i=1,imax-1
            u(i,j)=0.0
            du(i,j)=0.0
         enddo
         u(imax,j)=umax
      enddo

      dumax = 1000.0
      it=1

      ! e = etime(t)
      call system_clock ( clock_count1, clock_rate, clock_max )
      ! http://www.math.iastate.edu/burkardt/f_src/timer/timer.html 
      ! Main computation loop
      do while (dumax.ge.tol)
         dumax=0.0
         do j=2,jm1
            do i=2,im1
               du(i,j)=0.25*(u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1))-u(i,j)
               dumax=max(dumax,abs(du(i,j)))
               u(i,j)=u(i,j)+du(i,j)
            enddo
         enddo
      !   write (1,*) it,dumax
         it=it+1
      enddo
      
      !  e = etime(t)
      call system_clock ( clock_count2, clock_rate, clock_max )
      write (1,*) 'Time=',(clock_count2-clock_count1)/ real(clock_rate)
      write (1,*) it, dumax
      print *, 'Time = ',(clock_count2-clock_count1)/ real (clock_rate) 
      print *, it, dumax 
      stop
      end