!   The distribution of heat over time is described by the so called heat equation:
!   df/dt = alpha*(d2f/dx2 + d2f/dy2)   for a function f(x,y,t).
!   This formula may be discretized in a regular grid G(:,:) by computing the new
!      value G1(x,y) in a point (x,y) at each time step as:
!      G1(x,y) = G(x,y) + Ax * ( G(x+1,y) + G(x-1,y) - 2.0 * G(x,y))
! 		+ Ay * ( G(x,y+1) + G(x,y-1) - 2.0 * G(x,y))
!   For each point in the grid the next value depends on the values of the four up
!      and down, left and right adjacent points.
!
!   This solution is based on a source code written by 
!      Blaise Barney, D. Turner, George L. Gusciora
! ****************************************************************************
! Explanation of constants and variables
!   Nrows                          =  x dimension of problem grid 
!   Ncols                          =  y dimension of problem grid
!   Tsteps                         =  number of time Tsteps
!   Ax, Ay                         =  used in heat equation
!   Grid                           =  array for grids
! ---------------------------------------------------------------------------
MODULE heat2D_data
      IMPLICIT NONE
      integer, PARAMETER :: Nrows=2000,Ncols=2000
      integer, parameter :: Tsteps = 1000

      real, dimension(:,:,:), allocatable :: Grid

      real(4), parameter ::  Ax = 0.1, Ay = 0.1
      
CONTAINS
	
! ***************************************************************************
      subroutine compute (Rows, Cols, g0, g1)
! ***************************************************************************
      integer, intent(in) :: Rows, Cols
      real(4), dimension(Rows,Cols), intent(in) :: g0
      real(4), dimension(Rows,Cols), intent(out) :: g1

      integer :: i, j

      do j=2, Cols-1
        do i=2, Rows-1
            g1(i,j) = g0(i,j) &
     &       + Ax * ( g0(i+1,j) + g0(i-1,j) - 2.0 * g0(i,j)) &
     &       + Ay * ( g0(i,j+1) + g0(i,j-1) - 2.0 * g0(i,j))
        end do
      end do
   end subroutine compute


! ****************************************************************************
      subroutine inidat
! ****************************************************************************
      integer :: i,j

!     Initialize data
      Grid = 0.0

      do i=1,Nrows-1
        do j=Ncols/2-1,Ncols/2+1
          Grid(i+1,j+1,1) = float(i*(Nrows-i-1) + j*(Ncols-j-1))
        end do
      end do

      do i=Nrows/2-5,Nrows/2+5
        do j=Ncols/2-5,Ncols/2+5
          Grid(i+1,j+1,1) = float(i*(Nrows-i-1) + j*(Ncols-j-1))
        end do
      end do
      
      return
 end subroutine inidat

subroutine RealData2pgm(s1,s2,rdata,name)
!  Simple subroutine to dump real data in a PGM format
   implicit NONE
   INTEGER, INTENT(IN) :: s1, s2
   REAL, DIMENSION(s1,s2), intent(in) :: rdata
   character(*), intent(in) :: name
!
   integer :: i, j, ouni
   integer :: vp, vs
   real :: rmin, rmax
   character(80) :: fname
!
!  Write on unit 700 with  PGM format
   ouni = 700
   fname = TRIM(name)//".pgm"
   open(ouni,file=trim(fname),status="replace",iostat=vs)
   if ( vs /= 0 ) then
      print*,"!!!! Error write access to file ",trim(fname)
   end if
!  Magic code
   write(ouni,"(a)") "P2"
!  Dimensions
   write(ouni,*) s1, s2
!  Maximum value
   write(ouni,"(a)") "255"
!  Values from 0 to 255
   rmin = MINVAL(rdata); rmax = MAXVAL(rdata)
!
   vs = 0
   do i = 1, s1
      do j = 1, s2
         vp = INT ( (rdata(i,j) - rmin) * 255.0 / (rmax - rmin) )
         vs = vs + 1
         write(ouni,"(i4)",advance="no") vp
         if (vs >= 10 ) then
            write(ouni,"(a)") " "
            vs = 0
         end if
      end do
      write(ouni,"(a)") " "
      vs = 0
   end do
   close(ouni)
   return
end subroutine RealData2pgm

END MODULE heat2D_data

program heat2D
      USE heat2D_data
      IMPLICIT NONE

      CHARACTER(80) :: fname
      integer :: st, sp, it

      allocate(Grid(Nrows,0:Ncols+1,2), STAT=st)
           if ( st /= 0 ) then 
              WRITE(*,*) " Error allocating Grid(",Nrows,",",Ncols,",2)"
              WRITE(*,*) " Program terminates "
              stop
           end if


          WRITE(*,*) " Start computing ... "

          call inidat
!         Dump initial configuration
          it = 0 
          WRITE(fname,"(A,I4.4)") "tmap",it 
          call RealData2pgm(Nrows,Ncols,Grid(1,1,1),fname)

!        Compute iterations: 
          sp=1
          do it = 1, Tsteps
!          Calculate  grid points
           call compute(Nrows,Ncols,Grid(1,1,sp),Grid(1,1,3-sp))
           sp=3-sp
         end do

!           Print and show results
              
         it = Tsteps 
         WRITE(fname,"(A,I4.4)") "tmap",it 
         call RealData2pgm(Nrows,Ncols,Grid(1,1,sp),fname)

          WRITE(*,*) "  ... computed!"
   
         stop
 
end program heat2D