program main

!*****************************************************************************80
!
!! MAIN is the main program for HEAT.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    08 September 2013
!
!  Author:
!
!    John Burkardt
!
  use mpi

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

  call MPI_Comm_rank ( MPI_COMM_WORLD, id, ierr )

  call MPI_Comm_size ( MPI_COMM_WORLD, p, ierr )
!
!  Determine the portion of [A,B] to be assigned to processor ID.
!
  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 )
!
!  Terminate MPI.
!
  call MPI_Finalize ( ierr )
!
!  Terminate.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'HEAT:'
  write ( *, '(a)' ) '  Normal end of execution.'

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

!*****************************************************************************80
!
!! HEAT_PART solves the heat equation over "part" of the domain.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    08 September 2013
!
!  Author:
!
!    John Burkardt
!
  use mpi

  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 )
!
!  Time parameters:
!
  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 )
!
!  Space parameters.
!
  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
!
!  Set the initial value of H.
!
  h(0:n+1) = 95.0
!
!  Check the CFL condition.
!
  cfl = k * t_del / x_del / x_del

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

  do j = 1, j_max
!
!  Determine new time.
!
    t_new = ( dble (         j - j_min ) * t_max   &
            + dble ( j_max - j         ) * t_min ) &
            / dble ( j_max     - j_min )
!
!  To set H_NEW(1:N), update the temperature based on the four point stencil.
!
    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
!
!  Send MESSAGE #1, containing H_NEW(N) to neighbor ID+1.
!
    tag = 1

    if ( id < p - 1 ) then
      call MPI_Send ( h_new(n), 1, MPI_DOUBLE_PRECISION, id+1, tag, &
        MPI_COMM_WORLD, ierr )
    end if
!
!  Receive MESSAGE #1, containing H_NEW(0) from neighbor ID-1.
!
    if ( 0 < 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
!
!  Send MESSAGE #2, containing H_NEW(1) to neighbor ID-1.
!
    tag = 2

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

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

  end do

  wtime = MPI_Wtime ( ) - wtime

  if ( id == 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,g14.6)' ) '  Wall clock elapsed seconds = ', wtime
  end if
!
!  Print solution for last value of T at end of computation.
!  
  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