      program grid_test

c simple grid relaxation program to demo
c FORGE90 distributed memory parallelizer
      parameter (nmax = 500)
      real  w(2 * (nmax + 2) * (nmax + 2))
      real*8    second,timein
      open (unit = 5,file = 'grid.inp')
      open (unit = 6,file = 'grid.out')
10      print *,'Enter grid size (max=',nmax,'):'

c     read *,n
      n = 300
      if (n .lt. 1 .or. n .gt. nmax) then 
          print *,'n=',n,' is out of range.'
         goto 10
       endif
      print *,'Enter ncycles:'

c     read *,ncycles
      ncycles = 20
      call iterate(n,ncycles,w(1))
      end

      subroutine iterate (n, ncycles, grid)
      real*8    second,timein
      parameter (nmax = 500)
      real  grid(n + 2,n + 2,2)
      real  t(nmax + 2)
      real  dummy(2)
      iold = 1
      inew = 2

c initialize grid to 1.0
      timein = second ()
      do  j = 1,n + 2
         do  i = 1,n + 2
            grid(i,j,iold) = 1.0
         enddo
      enddo

c put a spike at the center
      grid(n / 2 + 2,n / 2 + 2,iold) = 50.0

c put a smaller spike at lower right corner
      grid(n + 1,n + 1,iold) = 40.0
      grid(1,n + 1,iold) = 40.0
      grid(n + 1,1,iold) = 40.0
      grid(1,1,iold) = 40.0
      do  it = 1,ncycles

CAPR$ IGNORE ALL INHIBITORS
CAPR$ IGNORE SYNC COM
CAPR$ DOPAR
         do  j = 2,n + 1
            do  i = 2,n + 1

c use exponential of averarage of logs to increase cpu time
               grid(i,j,inew) = exp ((alog (grid(i - 1,j - 1,iold)) + al
     +og (grid(i - 1,j,iold)) + alog (grid(i - 1,j + 1,iold)) + alog (gr
     +id(i,j - 1,iold)) + alog (grid(i,j,iold)) + alog (grid(i,j + 1,iol
     +d)) + alog (grid(i + 1,j - 1,iold)) + alog (grid(i + 1,j,iold)) + 
     +alog (grid(i + 1,j + 1,iold))) / 9.0)
            enddo
         enddo

c copy boundaries so they wrap as a torus
         do  j = 2,n + 1
            grid(1,j,inew) = grid(n + 1,j,inew)
         enddo
         do  j = 2,n + 1
            grid(n + 2,j,inew) = grid(2,j,inew)
         enddo

CAPR$ DOPAR
         do  i = 1,n + 2
            grid(i,1,inew) = grid(i,n + 1,inew)
         enddo

CAPR$ DOPAR
         do  i = 1,n + 2
            grid(i,n + 2,inew) = grid(i,2,inew)
         enddo

c switch iold and inew
         iold = 3 - iold
         inew = 3 - inew
      enddo
      gmax =  (-(1.0e30))
      gmin = 1.0e30
      ilo = 2
      ihi = n + 1
      if (ihi - ilo + 1 .gt. 40) then 
          is = (ihi - ilo + 1 - 40) / 2
         ilo = ilo + is
         ihi = ihi - is
       endif
      do  j = ilo,ihi
         do  i = ilo,ihi
            if (gmax .lt. grid(i,j,iold))             gmax = grid(i,j,io
     +ld)
 
            if (gmin .gt. grid(i,j,iold))             gmin = grid(i,j,io
     +ld)
 
         enddo
      enddo
      if (gmax .eq. gmin)       gmax = gmin + 1.0
 

c print histogram
      do  j = ilo,ihi
         do  i = ilo,ihi
            t(i) = grid(i,j,iold)
         enddo
         print 900,(int ((t(i) - gmin) * 9.0 / (gmax - gmin)), i = ilo,i
     +hi)
900             format(1000i1)
      enddo

c compute and print checksum
      checksum = 0.0
      do  j = 2,n + 1
         do  i = 2,n + 1
            checksum = checksum + grid(i,j,iold)
         enddo
      enddo
      print *,'Check sum = ',checksum

c print CPU time (node 0 for parallel version)
C     print *,"Time=",float(mclock())/100.0
      print *,'Time=',second () - timein
      end

