!
!   Block_array
!
!   An example of MPI derived data type
!
!
module mpi_stuff
    implicit none
    include 'mpif.h'
    !
    type  ::  msgtype
        sequence
        real(8) :: o
    end type msgtype
    !
    integer ::  mpi_type      !  New MPI data type
    !
    !   Fortran structure to be translated in MPI derived data type 
    !
    type  ::  vect_msgtype
        sequence
        type(msgtype), dimension(12) :: v
    end type vect_msgtype
    !
    integer ::  mpi_ncvect_type  !  New MPI data type
    !
    integer, dimension(MPI_STATUS_SIZE) :: v_status
    integer ::  error_code, proc_rank, mpi_procs = 0, mpi_world
    !
contains
    subroutine start_mpi()
        implicit none
        if (mpi_procs.gt.0) return
        call mpi_init(error_code)
        mpi_world = MPI_COMM_WORLD
        call mpi_comm_rank(mpi_world, proc_rank, error_code)
        call mpi_comm_size(mpi_world, mpi_procs, error_code)
        return
    end subroutine start_mpi
    !
    !   MPI derived type definition: 
    !
    subroutine build_ncvect_type()
        implicit none
        integer, parameter :: nblocks = 4
        integer, parameter :: clblock = 2
        integer, parameter :: lblock = 3
        integer :: elldatatype, rc
        !
        elldatatype = MPI_DOUBLE_PRECISION
        !
        call mpi_type_vector(nblocks, &
                             clblock,     &
                             lblock,     &
                             elldatatype, &
                             mpi_ncvect_type, rc )
        if(rc.ne.mpi_success) write(*, *)"Error mpi_type_vector ", rc
        !
        !   Before using mpi_ncvect_type it has to be MPI registered
        !
        call mpi_type_commit(mpi_ncvect_type, rc)
        if (rc.ne.mpi_success) write(*, *)"Error mpi_type_commit ", rc
        return
    end subroutine build_ncvect_type
    !
    subroutine release_ncvect_type()
        implicit none
        integer  ::  rc
        !
        ! Just release derived data type 
        !
        call mpi_type_free(mpi_ncvect_type, rc)
        if(rc.ne.mpi_success) write(*, *)"Error mpi_type_free ", rc
        return
    end subroutine release_ncvect_type
    !
    !
    subroutine initialize_ncvect_data(given)
        implicit none
        type(vect_msgtype), intent(out) :: given
        integer ::  k
        do k = 1, 12
           given%v(k)%o = k*1.d3
        end do
        return
    end subroutine initialize_ncvect_data
    !
    subroutine zero_ncvect_data(taken)
        implicit none
        type(vect_msgtype), intent(out) :: taken
        integer ::  k
        do k = 1, 12
           taken%v(k)%o = 0.0d0
        end do
        return
    end subroutine zero_ncvect_data
    !
    subroutine update_vect_data(given, taken)
        implicit none
        type(vect_msgtype), intent(in)  :: taken
        type(vect_msgtype), intent(out) :: given
        !
        given%v(1:12)%o = taken%v(1:12)%o + 1.0d0
        return
    end subroutine update_vect_data
    !
end module mpi_stuff
!
!   Example program: round robin
!
program block_array
    use mpi_stuff
    implicit none
    
    integer :: next_proc, msg_tag = 111, prec_proc
    integer :: i
    type(vect_msgtype) :: given, taken
    
    call start_mpi()
    
    call build_ncvect_type()
    
    next_proc = mod(proc_rank+1, mpi_procs)
    prec_proc = mod(mpi_procs+proc_rank-1, mpi_procs)
    
    call initialize_ncvect_data(given)
    
    if(proc_rank.eq.0) &
       call mpi_send(given, 1, mpi_ncvect_type, next_proc, msg_tag, &
                  mpi_world, error_code)
                  
    call mpi_recv(taken, 1, mpi_ncvect_type, prec_proc, msg_tag, &
                  mpi_world, v_status, error_code)
                  
    if(proc_rank.ne.0) then
        call update_vect_data(given, taken)
        call mpi_send(given, 1, mpi_ncvect_type, next_proc, msg_tag, &
                  mpi_world, error_code)
    else
        write(*, *)"Sent data ="
        do i = 1, 12
           write(*, "(1x, g8.2)", advance="NO") given%v(i)%o
           if ( mod(i, 6) == 0 ) write(*, "(1x)")
        end do
        
        write(*, *)"Received data ="
        do i = 1, 12
           write(*, "(1x, g8.2)", advance="NO") taken%v(i)%o
           if ( mod(i, 6) == 0 ) write(*, "(1x)")
        end do
    end if
    call mpi_finalize(error_code)
    stop
end program block_array