#include "options.inc"
c
c--------------------------------------------------------------
c      Routines for communication between processors
c      based on MPI 1.0
c
c             C.Eden Feb. 2002 (mailto ceden@phys.ocean.dal.ca)
c--------------------------------------------------------------
c
#ifdef no_mpp
c
c--------------------------------------------------------------
c     some dummy routines, if MPI service is disabled
c--------------------------------------------------------------
c
      subroutine mpi_comm_rank(comm,my_pe,ierr)
      implicit none
      integer comm,my_pe,ierr
      ierr=0
      my_pe=0
      end subroutine mpi_comm_rank

      subroutine mpi_comm_size(comm,n_pes,ierr)
      implicit none
      integer comm,n_pes,ierr
      ierr=0
      n_pes=1
      end subroutine mpi_comm_size

      subroutine mpi_init(ierr)
      integer ierr
      end subroutine mpi_init

      subroutine mpi_finalize(ierr)
      integer ierr
      end subroutine mpi_finalize
#endif



      subroutine mpi_test()
c
c--------------------------------------------------------------
c     test some basic mpi routines, i.e. the type declaration
c     in sub_mpif.h
c--------------------------------------------------------------
c
#ifndef no_mpp
      implicit none
      integer :: ierr,my_pe=-1,all_pes
      integer :: xint,xint2
      real    :: xreal,xreal2
#include "sub_mpif.h"
c
c   get some mpi infos first
c
      print*,' performing MPI test ...'
      call mpi_comm_rank(MPI_COMM_WORLD,my_pe,ierr)
      if (ierr/=0) goto 20
      call mpi_comm_size(MPI_COMM_WORLD,all_pes,ierr)
      if (ierr/=0) goto 20
c
c   try first global barrier
c
      print*,' tyring first global barrier for PE #',my_pe
      call mpi_barrier(mpi_comm_world, ierr)
      if (ierr/=0) goto 20
c
c   try broadcasting
c
      print*,' tyring global broadcast for PE #',my_pe
      xreal = 1.0
      call mpi_bcast(xreal,1,impireal,0,mpi_comm_world,ierr)
      if (ierr/=0) goto 20
      xint = 1
      call mpi_bcast(xint,1,impiinteger,0,mpi_comm_world,ierr)
      if (ierr/=0) goto 20
c
c   check results of broadcasting
c
      if (xreal /= 1.0 ) then
       print*,' fatal: MPI test failed on broadcasting reals for ',
     &             'PE #',my_pe
       call halt_stop(' in mpi_test')
      endif
      if (xint /= 1 ) then
       print*,' fatal: MPI test failed on broadcasting integer for ',
     &             'PE #',my_pe
       call halt_stop(' in mpi_test')
      endif
      call mpi_barrier(mpi_comm_world, ierr)
      if (ierr/=0) goto 20

c
c   try global sum
c
      xreal = 2.0
      call mpi_allreduce(xreal,xreal2,1,impireal,MPI_SUM,
     &            mpi_comm_world,ierr)
      if (ierr/=0) goto 20
      xint = 2.0
      call mpi_allreduce(xint,xint2,1,impiinteger,MPI_SUM,
     &            mpi_comm_world,ierr)
      if (ierr/=0) goto 20
c
c   check results 
c
      xreal = xreal2/all_pes
      if (xreal /= 2.0 ) then
       print*,' fatal: MPI test failed on global sum (real) for ',
     &             'PE #',my_pe
       call halt_stop(' in mpi_test')
      endif
      xint = xint2/all_pes
      if (xint /= 2.0 ) then
       print*,' fatal: MPI test failed on global sum (int) for ',
     &             'PE #',my_pe
       call halt_stop(' in mpi_test')
      endif
      call mpi_barrier(mpi_comm_world, ierr)
      if (ierr/=0) goto 20

c
c   try point to point protocol
c
c   to be done
      call mpi_barrier(mpi_comm_world, ierr)
      if (ierr/=0) goto 20

      print*,' MPI test okay for PE #',my_pe,' .'

      return
      stop

 20   continue
      print*,' fatal: MPI test failed on'
      print*,' PE #',my_pe,' error message:'
      call print_mpi_error(ierr)
      call halt_stop(' in mpi_test')
#else
      print*,' performing MPI test ...'
      print*,' no MPI service present ',
     &         '(disable cpp option no_mpp to enable MPP)'
      return
#endif
      end subroutine mpi_test


      subroutine split_pes(domains,n_pes_domains,
     &   my_comm,sub_domain,sub_domain_pe0)
c
c--------------------------------------------------------------
c     split all PEs into several groups, which
c     are the sub domains
c--------------------------------------------------------------
c
      implicit none
      integer, intent(in)  :: domains                ! number of domains
      integer, intent(in)  :: n_pes_domains(domains) ! number of pes in domains
      integer, intent(out) :: sub_domain             ! number of sub domain
      integer, intent(out) :: my_comm                ! local communicator
      integer, intent(out) :: sub_domain_pe0(domains)! rank of leading PE of 
                                                     ! sub domain in global 
                                                     ! address space
#ifdef no_mpp
      sub_domain=0
      my_comm=1
      sub_domain_pe0=0

#else
# include "sub_mpif.h"
      integer ierr,all_pes,my_pe
      integer color,key,n,k
      call mpi_barrier(mpi_comm_world, ierr)
      call mpi_comm_size(MPI_COMM_WORLD,all_pes,ierr)
      call mpi_comm_rank(MPI_COMM_WORLD,my_pe,ierr)


      k=0
      do n=1,domains
        k=k+n_pes_domains(n)
      enddo
      if (k /= all_pes) then
        if (my_pe == 0) then
         print*,' Number of all PEs (=',all_pes,')'
         print*,' does not match number of PEs in sub domains'
         do n=1,domains
          print*,' nr of PEs in domain ',n,' = ',n_pes_domains(n)
         enddo
        endif
        call halt_stop(' in split_pes ')
      endif


      k=0
      color=-1
      do n=1,domains
        sub_domain_pe0(n)=k
        k=k+n_pes_domains(n)
        if (my_pe < k.and. color == -1 ) then
          color=n
          key=my_pe-(k-n_pes_domains(n))
c          key=0
          sub_domain=n-1
        endif
      enddo
      call mpi_comm_split(MPI_COMM_WORLD,color,key,my_comm,ierr)
      call print_mpi_error(ierr)
c      print*,'color=',color,'key=',key,'my_pe=',my_pe,'comm',my_comm
      call mpi_barrier(mpi_comm_world, ierr)


#endif
      end subroutine split_pes



      subroutine print_mpi_error(ierr)
c
c--------------------------------------------------------------
c     print out MPI error message
c--------------------------------------------------------------
c
#ifndef no_mpp
      implicit none
      integer :: ierr,ierr2,len
      character*80 :: string
      if (ierr/=0) then
       call MPI_ERROR_STRING(ierr,string, LEN, ierr2)
       print*,' MPI ERROR:',ierr
       print*,'           ',string(1:len)
      endif
#else
      implicit none
      integer :: ierr
#endif
      end subroutine print_mpi_error



      subroutine halt_stop(string)
c--------------------------------------------------------------
c     controlled stop
c--------------------------------------------------------------
      implicit none
      character*(*) :: string
      integer :: ierr,code,my_pe
#ifndef no_mpp
#include "sub_mpif.h"
#endif

#ifdef no_mpp
      print*,string
#else
      call mpi_comm_rank(MPI_COMM_WORLD,my_pe,ierr)
      print*,' global pe #',my_pe,' : ',string
      print*,' global pe #',my_pe,' aborting '
#endif

#ifdef no_mpp
      call abort()
#else
      code=99
      call MPI_ABORT(mpi_comm_world, code, IERR)
#endif

      end subroutine halt_stop



      subroutine barrier
c
c--------------------------------------------------------------
c     A barrier for the local sub domain
c--------------------------------------------------------------
c
#ifndef no_mpp
      use spflame_module
      implicit none
      integer :: ierr

      call mpi_barrier(my_comm, ierr)
      call print_mpi_error(ierr)
#endif
      end subroutine barrier




      subroutine bcast_real(x,len,pe)
c
c--------------------------------------------------------------
c     Broadcast a real vector to the local sub domain
c--------------------------------------------------------------
c
#ifndef no_mpp
      use spflame_module
      implicit none
      integer len,ierr,pe
      real x(len)
#include "sub_mpif.h"
      call mpi_barrier(my_comm, ierr)
      call mpi_bcast(x,len,impireal,pe,my_comm,ierr)
#else
      implicit none
      integer :: len,pe
      real :: x
#endif
      end subroutine bcast_real

      subroutine bcast_integer(x,len,pe)
c
c--------------------------------------------------------------
c     Broadcast a integer vector to the local sub domain
c--------------------------------------------------------------
c
#ifndef no_mpp
      use spflame_module
      implicit none
      integer len,ierr,pe
      integer x(len)
#include "sub_mpif.h"
      call mpi_barrier(my_comm, ierr)
      call mpi_bcast(x,len,impiinteger,pe,my_comm,ierr)
#else
      implicit none
      integer :: len,x,pe
#endif
      end subroutine bcast_integer

      subroutine bcast_char(x,len,pe)
c
c--------------------------------------------------------------
c     Broadcast a character vector to the local sub domain
c--------------------------------------------------------------
c
#ifndef no_mpp
      use spflame_module
      implicit none
      integer len,ierr,pe
      character x(*)
#include "sub_mpif.h"
      call mpi_barrier(my_comm, ierr)
      call mpi_bcast(x,len,mpi_character,pe,my_comm,ierr)
#else
      implicit none
      integer :: len,pe
      character x(*)
#endif
      end subroutine bcast_char


      subroutine bcast_logical(x,len,pe)
c
c--------------------------------------------------------------
c     Broadcast a logical vector to the local sub domain
c--------------------------------------------------------------
c
#ifndef no_mpp
      use spflame_module
      implicit none
      integer len,ierr,pe
      logical x(len)
#include "sub_mpif.h"
      call mpi_barrier(my_comm, ierr)
      call mpi_bcast(x,len,mpi_logical,pe,my_comm,ierr)
#else
      implicit none
      integer :: len,pe
      logical :: x
#endif
      end subroutine bcast_logical





      subroutine bcast_world_real(x,len,pe)
c
c--------------------------------------------------------------
c     Broadcast a real vector to all sub domains
c--------------------------------------------------------------
c
#ifndef no_mpp
      use spflame_module
      implicit none
      integer len,ierr,pe
      real x(len)
#include "sub_mpif.h"
      call mpi_barrier(mpi_comm_world, ierr)
      call mpi_bcast(x,len,impireal,pe,mpi_comm_world,ierr)
#else
      implicit none
      integer :: len,pe
      real :: x
#endif
      end subroutine bcast_world_real

      subroutine bcast_world_integer(x,len,pe)
c
c--------------------------------------------------------------
c     Broadcast an integer vector to all sub domains
c--------------------------------------------------------------
c
#ifndef no_mpp
      use spflame_module
      implicit none
      integer len,ierr,pe
      integer x(len)
#include "sub_mpif.h"
      call mpi_barrier(mpi_comm_world, ierr)
      call mpi_bcast(x,len,impiinteger,pe,mpi_comm_world,ierr)
#else
      implicit none
      integer :: len,pe
      integer x
#endif
      end subroutine bcast_world_integer

      subroutine bcast_world_logical(x,len,pe)
c
c--------------------------------------------------------------
c     Broadcast a logical vector to all sub domains
c--------------------------------------------------------------
c
#ifndef no_mpp
      use spflame_module
      implicit none
      integer len,ierr,pe
      logical x(len)
#include "sub_mpif.h"
      call mpi_barrier(mpi_comm_world, ierr)
      call mpi_bcast(x,len,mpi_logical,pe,mpi_comm_world,ierr)
#else
      implicit none
      integer :: len,pe
      logical x
#endif
      end subroutine bcast_world_logical



      subroutine recv_integer(x,len,pe,tag)
c
c--------------------------------------------------------------
c     Reveive an integer vector from PE pe in local sub domain
c--------------------------------------------------------------
c
      use spflame_module
      implicit none
      integer len,ierr,pe,tag
      integer x(len)
#ifndef no_mpp
#include "sub_mpif.h"
      integer :: status(MPI_STATUS_SIZE)
      call mpi_recv(x,len,MPI_INTEGER,pe,tag,my_comm,Status,ierr)
#endif
      end subroutine recv_integer

      subroutine send_integer(x,len,pe,tag)
c
c--------------------------------------------------------------
c     Send an integer vector to PE pe in local sub domain
c--------------------------------------------------------------
c
      use spflame_module
      implicit none
      integer len,ierr,pe,tag
      integer x(len)
#ifndef no_mpp
#include "sub_mpif.h"
      call mpi_send (x,len,MPI_INTEGER,pe,tag,my_comm,ierr)
#endif
      end subroutine send_integer


      subroutine recv_world_integer(x,len,pe,tag)
c
c--------------------------------------------------------------
c     Receive an integer vector from PE pe in global space
c--------------------------------------------------------------
c
      use spflame_module
      implicit none
      integer len,ierr,pe,tag
      integer x(len)
#ifndef no_mpp
#include "sub_mpif.h"
      integer :: status(MPI_STATUS_SIZE)
      call mpi_recv(x,len,MPI_INTEGER,pe,tag,
     &              mpi_comm_world,Status,ierr)
#endif
      end subroutine recv_world_integer

      subroutine send_world_integer(x,len,pe,tag)
c
c--------------------------------------------------------------
c     Send an integer vector to PE pe in global space
c--------------------------------------------------------------
c
      use spflame_module
      implicit none
      integer len,ierr,pe,tag
      integer x(len)
#ifndef no_mpp
#include "sub_mpif.h"
      call mpi_send (x,len,MPI_INTEGER,pe,tag,
     &               mpi_comm_world,ierr)
#endif
      end subroutine send_world_integer


      subroutine global_sum(x)
c
c--------------------------------------------------------------
c     Do a sum of real x over all PEs in sub domain
c--------------------------------------------------------------
c
#ifndef no_mpp
      use spflame_module
      implicit none
#include "sub_mpif.h"
      real    :: x, x_sym,x_sym2
      integer :: ierr
c     skip this copying overhead  here
      x_sym = x
      call mpi_allreduce(x_sym,x_sym2,1,impireal,MPI_SUM,my_comm,ierr)
      call print_mpi_error(ierr)
      x = x_sym2
#else
      implicit none
      real :: x
#endif
      end subroutine global_sum

      subroutine global_sum_int(x)
c
c--------------------------------------------------------------
c     Do a sum of integer x over all PEs in sub domain
c--------------------------------------------------------------
c
#ifndef no_mpp
      use spflame_module
      implicit none
#include "sub_mpif.h"
      integer    :: x, x_sym,x_sym2
      integer :: ierr
c     skip this copying overhead  here
      x_sym = x
      call mpi_allreduce(x_sym,x_sym2,1,impiinteger,
     &                   MPI_SUM,my_comm,ierr)
      call print_mpi_error(ierr)
      x = x_sym2
#else
      implicit none
      real :: x
#endif
      end subroutine global_sum_int



      subroutine global_sum_vec(x,len)
c
c--------------------------------------------------------------
c     Do a sum of real vecter x(len) over all PEs in sub domain
c--------------------------------------------------------------
c
#ifndef no_mpp
      use spflame_module
      implicit none
#include "sub_mpif.h"
      integer :: ierr,len
      real    :: x(len), x_sym(len),x_sym2(len)
c     skip this copying overhead  here
      x_sym(1:len) = x(1:len)
      call mpi_allreduce(x_sym,x_sym2,len,impireal,MPI_SUM,my_comm,ierr)
      call print_mpi_error(ierr)
      x(1:len) = x_sym2(1:len)
#else
      implicit none
      real :: x
      integer :: len
#endif
      end subroutine global_sum_vec



      subroutine global_lor(x)
c
c--------------------------------------------------------------
c     Do a OR of logical x over all PEs in sub domain
c--------------------------------------------------------------
c
#ifndef no_mpp
      use spflame_module
      implicit none
#include "sub_mpif.h"
      logical    :: x, x_sym,x_sym2
      integer :: ierr
c     skip this copying overhead  here
      x_sym = x
      call mpi_allreduce(x_sym,x_sym2,1,mpi_logical,
     &                   MPI_LOR,my_comm,ierr)
      call print_mpi_error(ierr)
      x = x_sym2
#else
      implicit none
      logical :: x
#endif
      end subroutine global_lor


      subroutine global_lor_vec(x,len)
c
c--------------------------------------------------------------
c     Do a OR of logical x over all PEs in sub domain
c--------------------------------------------------------------
c
#ifndef no_mpp
      use spflame_module
      implicit none
#include "sub_mpif.h"
      logical    :: x(len), x_sym(len),x_sym2(len)
      integer :: ierr,len
c     skip this copying overhead  here
      x_sym = x
      call mpi_allreduce(x_sym,x_sym2,len,mpi_logical,
     &                   MPI_LOR,my_comm,ierr)
      call print_mpi_error(ierr)
      x = x_sym2
#else
      implicit none
      logical :: x
      integer :: len
#endif
      end subroutine global_lor_vec




      subroutine global_max(x)
c
c--------------------------------------------------------------
c     Get the max of real x over all PEs in sub domain
c--------------------------------------------------------------
c
#ifndef no_mpp
      use spflame_module
      implicit none
#include "sub_mpif.h"
      real    :: x, x_sym,x_sym2
      integer :: ierr
c     skip this copying overhead  here
      x_sym = x
      call mpi_allreduce(x_sym,x_sym2,1,impireal,MPI_MAX,my_comm,ierr)
      call print_mpi_error(ierr)
      x = x_sym2
#else
      implicit none
      real :: x
#endif
      end subroutine global_max

      subroutine global_min(x)
c
c--------------------------------------------------------------
c     Get the a min of real x over all PEs in sub domain
c--------------------------------------------------------------
c
#ifndef no_mpp
      use spflame_module
      implicit none
#include "sub_mpif.h"
      real    :: x, x_sym,x_sym2
      integer :: ierr
c     skip this copying overhead  here
      x_sym = x
      call mpi_allreduce(x_sym,x_sym2,1,impireal,MPI_MIN,my_comm,ierr)
      x = x_sym2
#else
      implicit none
      real :: x
#endif
      end subroutine global_min


