program main

c*********************************************************************72
c
cc MAIN is the main program for HEAT.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    08 September 2013
c
c  Author:
c
c    John Burkardt
c
      include 'mpif.h'

      double precision a
      parameter ( a = 0.0D+00 )
      double precision b
      parameter ( b = 1.0D+00 )
      integer id
      integer ierr
      integer n
      integer p
      double precision x_max
      double precision x_min
c
c  Startup:
c
      call MPI_Init ( ierr )

      call MPI_Comm_rank ( MPI_COMM_WORLD, id, ierr )

      call MPI_Comm_size ( MPI_COMM_WORLD, p, ierr )
c
c  Determine the portion of [A,B] to be assigned to processor ID.
c
      n = 12
      i = 0

      x_min = ( dble ( p * n + 1 - id * n - i ) * a   
     &        + dble (             id * n + i ) * b ) 
     &        / dble ( p * n + 1              )

      i = n + 1

      x_max = ( dble ( p * n + 1 - id * n - i ) * a   
     &        + dble (             id * n + i ) * b ) 
     &        / dble ( p * n + 1              )

      call heat_part ( n, p, id, x_min, x_max )
c
c  Shut down.
c
      call MPI_Finalize ( ierr )

      stop
      end
      subroutine heat_part ( n, p, id, x_min, x_max )

c*********************************************************************72
c
cc HEAT_PART solves the heat equation over "part" of the domain.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    08 September 2013
c
c  Author:
c
c    John Burkardt
c
      include 'mpif.h'

      integer n

      double precision cfl
      double precision h(0:n+1)
      double precision h_new(0:n+1)
      integer i
      integer id
      integer ierr
      integer j
      integer j_max
      integer j_min
      double precision k
      integer p
      integer status(MPI_STATUS_SIZE)
      double precision t
      double precision t_del
      double precision t_max
      double precision t_min
      double precision t_new
      integer tag
      double precision wtime
      double precision x(0:n+1)
      double precision x_del
      double precision x_max
      double precision x_min

      k = 0.002 / dble ( p )
c
c  Time parameters:
c
      j_min = 0
      j_max = 100
      t_min = 0.0
      t_max = 10.0
      t_del = ( t_max - t_min ) / dble ( j_max - j_min )
c
c  Space parameters.
c
      x_del = ( x_max - x_min ) / dble ( n + 1 )
      do i = 0, n + 1
        x(i) = ( dble (         i ) * x_max   
     &         + dble ( n + 1 - i ) * x_min ) 
     &         / dble ( n + 1     )
      end do
c
c  Set the initial value of H.
c
      h(0:n+1) = 95.0
c
c  Check the CFL condition.
c
      cfl = k * t_del / x_del / x_del

      if ( 0.5 .le. cfl ) then
        write ( *, '(a)' ) '  CFL condition failed.'
        stop
      end if
c
c  Each execution of this loop computes the solution at the next time.
c
      wtime = MPI_Wtime ( )

      do j = 1, j_max
c
c  Determine new time.
c
        t_new = ( dble (         j - j_min ) * t_max   
     &          + dble ( j_max - j         ) * t_min ) 
     &          / dble ( j_max     - j_min )
c
c  To set H_NEW(1:N), update the temperature based on the four point stencil.
c
        do i = 1, n
          h_new(i) = h(i) + t_del * ( 
     &      k * ( h(i-1) - 2.0 * h(i) + h(i+1) ) / x_del / x_del 
     &      + 2.0 * sin ( x(i) * t ) )
        end do
c
c  Send MESSAGE #1, containing H_NEW(N) to neighbor ID+1.
c
        tag = 1

        if ( id .lt. p - 1 ) then
          call MPI_Send ( h_new(n), 1, MPI_DOUBLE_PRECISION, id+1, tag, 
     &      MPI_COMM_WORLD, ierr )
        end if
c
c  Receive MESSAGE #1, containing H_NEW(0) from neighbor ID-1.
c
        if ( 0 .lt. id ) then
          call MPI_Recv ( h_new(0), 1, MPI_DOUBLE_PRECISION, id-1, tag, 
     &      MPI_COMM_WORLD, status, ierr )
        else
          h_new(0) = 100.0 + 10.0 * sin ( t_new )
        end if
c
c  Send MESSAGE #2, containing H_NEW(1) to neighbor ID-1.
c
        tag = 2

        if ( 0 .lt. id ) then
          call MPI_Send ( ?, ?, ?, ?, ?, 
     &      ?, ? )
        end if
c
c  Receive MESSAGE #2, containing H_NEW(N+1) from neighbor ID+1.
c
        if ( id .lt. p - 1 ) then
          call MPI_Recv ( ?, ?, ?, ?, ?, 
     &      ?, ?, ? )
        else
          h_new[n+1] = ?
        end if
c
c  Update time and temperature.
c
        t = t_new

        do i = 0, n + 1
          h(i) = h_new(i)
        end do

      end do

      wtime = MPI_Wtime ( ) - wtime

      if ( id .eq. 0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a,g14.6)' ) 
     &    '  Wall clock elapsed seconds = ', wtime
      end if
c
c  Print solution for last value of T at end of computation.
c  
      write ( *, '(i2, 2x, a,   f7.2)' ) id, 'T=', t
      write ( *, '(i2, 2x, a, 14f7.2)' ) id, 'X=', x(0:n+1)
      write ( *, '(i2, 2x, a, 14f7.2)' ) id, 'H=', h(0:n+1)

      return
      end