!<CONTACT EMAIL="petri@pik-potsdam.de"> Stefan Petri
!</CONTACT>
!<REVIEWER EMAIL="">
!
!</REVIEWER>
!-----------------------------------------------------------------------
!<OVERVIEW>
!  Test program for development of the Fortran-to-python wrapper
!  for the atmospheric model Aeolus2.
!  This program tests MPI process management and
!  passing data, especially arrays, forth and back
!  betweeen Fortran and python.
!</OVERVIEW>

!<DESCRIPTION>
! mpiifort -g -warn -warn noerrors -o test_MPI_Spawn test_MPI_Spawn.F90
!
!  Fortran-to-python wrapper for the atmospheric model Aeolus2
!  Since it was not feasible to get the python interpreter running
!  compiled and linked into the FMS/Fortran executable, we now
!  start it as separate executable via MPI_Comm_spawn,
!  and communicate with that via quasi-remote-procedure calls,
!  i.e. all Fortran variables/arrays are packed into MPI
!  messages and sent to the Python side, where they are unpacked into
!  Python variables resp. numpy arrays. And vice vers. Sigh.
!</DESCRIPTION>

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!                                                                   !!
!!                   GNU General Public License                      !!
!!                                                                   !!
!! This file is part of the Flexible Modeling System (FMS).          !!
!!                                                                   !!
!! FMS is free software; you can redistribute it and/or modify       !!
!! it and are expected to follow the terms of the GNU General Public !!
!! License as published by the Free Software Foundation.             !!
!!                                                                   !!
!! FMS is distributed in the hope that it will be useful,            !!
!! but WITHOUT ANY WARRANTY; without even the implied warranty of    !!
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the     !!
!! GNU General Public License for more details.                      !!
!!                                                                   !!
!! You should have received a copy of the GNU General Public License !!
!! along with FMS; if not, write to:                                 !!
!!          Free Software Foundation, Inc.                           !!
!!          59 Temple Place, Suite 330                               !!
!!          Boston, MA  02111-1307  USA                              !!
!! or see:                                                           !!
!!          http://www.gnu.org/licenses/gpl.txt                      !!
!!                                                                   !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

program test_MPI_Spawn

  use mpi ! , only: ....
  use ifport, only: getpid

  implicit none

  !include "mpif.h"

  integer :: atmos_comm, pe, root_pe, npes
  integer :: inter_comm
  integer :: error

  !!
  !! main()
  !!

  call MPI_INIT(error)
  call MPI_COMM_RANK( MPI_COMM_WORLD, pe,   error )
  call MPI_COMM_SIZE( MPI_COMM_WORLD, npes, error )

  write(*,*) 'test_MPI_Spawn.F90 on PE',pe,' of ',npes,'PEs'

  call interface_check(MPI_COMM_WORLD)

  call finish(1, 2)

contains

  subroutine interface_check(atmos_comm_in)
    integer, intent(in) :: atmos_comm_in

    character(len=MPI_MAX_LIBRARY_VERSION_STRING) :: mpi_version_string

    ! some magical values for testing the interface between Fortran90 and python code.
    ! we pass them from Fortran90 to python functions and test in there if they were
    ! received correctly.
    integer, parameter :: deadbeef = X'eadbeef' ! X'deadbeef' yields overflow error on 32bit machines
    integer, parameter :: fse = 4711
    real,    parameter :: onetwothree = 1234567809.0987654321D0, &
         pitest = 3.1415926535897932384626433832795029D0 ! from linux <math.h>
    logical :: true = .true., false = .false., flag
    integer :: kind_i = KIND(deadbeef), kind_r = KIND(pitest), kind_l = KIND(true)
    real, dimension(2,3) :: testarray2d
    real, dimension(2,3) :: response2d ! to store the data that was sent back from the python side
    !real, dimension(2,3,4) :: testarray3d

    integer :: i, j, k, ierror
    integer :: inter_comm_size, universe_size = -1, remote_pe, remote_size = -1
    integer :: mypid
    character(len=25) :: command, argv(5)
    integer, dimension(:), allocatable :: remotepids
    integer :: status(MPI_STATUS_SIZE)
    real :: zef ! space to receive 08.15

#ifdef I_MPI_VERSION
#ifdef MPICH_VERSION
    write (*,*) 'compiled with I_MPI_VERSION', I_MPI_VERSION, 'and MPICH_VERSION', MPICH_VERSION
#else
    write (*,*) 'compiled with I_MPI_VERSION', I_MPI_VERSION
#endif
#else
#ifdef MPICH_VERSION
    write (*,*) 'compiled with MPICH_VERSION', MPICH_VERSION
#endif
#endif

    call MPI_Get_library_version(mpi_version_string, i, ierror)
    write (*,*) 'running on library version ', trim(mpi_version_string)

    atmos_comm = atmos_comm_in
    call MPI_Comm_rank(atmos_comm, pe, ierror)
    call MPI_Comm_size(atmos_comm, npes, ierror)
    mypid = getpid()

    call MPI_Comm_get_attr(atmos_comm, MPI_UNIVERSE_SIZE, universe_size, flag);
    if (.not. flag) write(*,*) 'This MPI does not support UNIVERSE_SIZE.'
    write(*,*) 'pid', mypid, 'task', pe, 'of', npes, 'in a universe of size', universe_size

    ! Now spawn the python-implemented atmos model
    command = 'python3'
    argv(1) = '-u'
    argv(2) = '-m'
    argv(3) = 'mpi4py'
    argv(4) = 'test_MPI_Spawn.py'
    argv(5) = ' '
    call MPI_Comm_spawn(command, argv, npes, MPI_INFO_NULL, 0, atmos_comm, inter_comm, MPI_ERRCODES_IGNORE, ierror)
    call MPI_Comm_size(inter_comm, inter_comm_size, ierror)
    call MPI_Comm_rank(inter_comm, remote_pe, ierror)
    call MPI_Comm_remote_size(inter_comm, remote_size, ierror)
    write (*,*) 'pid',mypid,'task',pe,'is task',remote_pe,'of remotesize',remote_size,'in inter_comm'
    !call MPI_Barrier(inter_comm, ierror)
    ! let python code print its hello message
    !call MPI_Barrier(inter_comm, ierror)

    allocate(remotepids(remote_size))
    do i = 1,remote_size
       call MPI_Recv(remotepids(i), 1, MPI_INT, MPI_ANY_SOURCE, 42, inter_comm, status, ierror)
       write (*,*) 'F90 task',pe,'received pid',remotepids(i),'from task',status(MPI_SOURCE)
    end do
    !call MPI_Barrier(inter_comm, ierror)
    do i = 1,remote_size
       call MPI_Send(mypid, 1, MPI_INT, i-1, 43, inter_comm, ierror)
    end do

    write(*,*) 'aeolus2fms_interface_check entering first barrier' ; call flush(6)
    call MPI_Barrier(inter_comm, ierror)

    ! now call out into the python - implemented heart of Aeolus2
    ! first assemble some parameters with well-known values to verify
    ! that the F90<->python interface is OK.
    do j=1,3
       do i=1,2
          testarray2d(i,j) = i+j/10.0
          write (*,*) i,j,testarray2d(i,j)
       end do
    end do
    !do k=1,4
    !   do j=1,3
    !      do i=1,2
    !         testarray3d(i,j,k) = i*100.0+j+k/10.0
    !         write (*,*) i,j,k,testarray3d(i,j,k)
    !      end do
    !   end do
    !end do

    write(*,*) 'MAXEXPONENT ', MAXEXPONENT(pitest) ; call flush(6)

    write(*,*) 'interface_check sending deadbeef' ; call flush(6)
    call MPI_Send(deadbeef, 1, MPI_INT, pe, 44, inter_comm, ierror)
    write(*,*) 'interface_check sending pitest' ; call flush(6)
    call MPI_Send(pitest, 1, MPI_DOUBLE, pe, 45, inter_comm, ierror)
    !ierror = pyargs%setitem(3, true); errcheck
    !ierror = pyargs%setitem(4, false); errcheck
    write(*,*) 'interface_check sending testarray2d' ; call flush(6)
    call MPI_Send(testarray2d, size(testarray2d), MPI_DOUBLE, pe, 46, inter_comm, ierror)
    !call MPI_Send(testarray3d, size(testarray3d), MPI_DOUBLE, pe, 46, inter_comm, ierror)
    !ierror = pyargs%setitem(6, kind_i); errcheck
    !ierror = pyargs%setitem(7, kind_r); errcheck
    !ierror = pyargs%setitem(8, kind_l); errcheck
    !
    !ierror = call_py_noret(aeolus2py, "interface_check", pyargs)
    !errcheck

    write(*,*) 'interface_check receiving one double' ; call flush(6)
    call MPI_Recv(zef, 1, MPI_DOUBLE, MPI_ANY_SOURCE, 47, inter_comm, status, ierror)
    if (8.15 /= zef) then
       write(*,*) 'interface_check expected 8.15, received', zef; call flush(6)
    end if
    call MPI_Recv(response2d, size(response2d), MPI_DOUBLE, MPI_ANY_SOURCE, 48, inter_comm, status, ierror)
    if (any(testarray2d /= response2d)) then
       write(*,*) 'interface_check response2d does not match testarray2d' ; call flush(6)
       do j=1,3
          do i=1,2
             write (*,*) i,j,testarray2d(i,j),response2d(i,j)
          end do
       end do
    end if

    deallocate(remotepids)

    write(*,*) 'interface_check entering final barrier' ; call flush(6)
    call MPI_Barrier(inter_comm, ierror)
    write(*,*) 'interface_check finished' ; call flush(6)
  end subroutine interface_check


  subroutine finish(Time_seconds, Time_days)
    integer, intent(in) :: Time_seconds, Time_days

    integer :: ierror
    integer, dimension(2) :: aeolus2params
    aeolus2params(1) = Time_seconds
    aeolus2params(2) = Time_days
    call MPI_Send(aeolus2params, size(aeolus2params), MPI_INT, pe, 72, inter_comm, ierror)
    write(*,*) 'aeolus2fms_finish sent params' ; call flush(6)

  end subroutine finish

end program test_MPI_Spawn
