module bmimodule
   use iso_c_binding
   use iso_c_utils
   use mpi
   implicit none
   save
   integer isize,irank
   integer, target :: comm, commwork
   real*8, allocatable ::  x(:)
   integer, target :: nx
   logical :: mpi_ready

contains
   subroutine decomp(n, numprocs, myid, s, e)

      !  From the mpich-distribution: MPE_Decomp1d
      !
      !  n        : integer, input (for example: number of tasks or
      !                               dimension of an array)
      !  numprocs : integer, input (number of processes)
      !  myid     : integer, input (MPI id of this process)
      !  s        : integer, output
      !  e        : integer, output
      !
      !  The purpose of this subroutine is to get an as equal as possible
      !  division of the integers 1..n among numprocs intervals.
      !  A common use is to divide work among MPI processes.
      !
      !  example:
      !
      !  integer numprocs, ier, myid, n, s, e, i
      !  call MPI_Comm_size(MPI_COMM_WORLD, numprocs, ier)
      !  call MPI_Comm_rank(MPI_COMM_WORLD, myid, ier)
      !  n = 1000
      !  call decomp(n, nprocs, myid, s, e)
      !  do i = s,e
      !    ...
      !  enddo
      !
      integer, intent(in) :: n, numprocs, myid
      integer, intent(out):: s, e
      integer nlocal
      integer deficit

      nlocal  = n / numprocs
      s       = myid * nlocal + 1
      deficit = mod(n,numprocs)
      s       = s + min(myid,deficit)
      if (myid .lt. deficit) then
         nlocal = nlocal + 1
      endif
      e = s + nlocal - 1
      if (e .gt. n .or. myid .eq. numprocs-1) then
         e = n
      endif

   end subroutine decomp


   integer(c_int) function set_var(c_var_name, xptr) result(ier) bind(C, name="set_var")
      character(kind=c_char), intent(in) :: c_var_name(*)
      type(c_ptr), value, intent(in)     :: xptr
      character(len=strlen(c_var_name))  :: var_name

      integer, pointer :: ipointer
      real*8,  pointer :: xpointer(:)
      real*8, target   :: dummyx(0)

      integer          :: i, ierror, rank,sz,s,e,d
      integer, allocatable :: sendcounts(:),displs(:)

      var_name = char_array_to_string(c_var_name)
      print *,irank,'xbmi set_var var_name = ',var_name

      select case(var_name)
       case("mpicomm")
         call c_f_pointer(xptr, ipointer)
         comm = ipointer
         call mpi_comm_rank(comm, rank, ierror)
         call mpi_comm_size(comm, sz, ierror)
         i = 0
         if (rank .eq. 0) i = sz
         call mpi_bcast(i, 1, mpi_integer, 0,comm,ierror)
         print *,'xbmi after setting comm and bcast, i = ',i
       case("nx")
         print *,'xbmi about to set nx'
         call c_f_pointer(xptr, ipointer)
         nx = ipointer
         print *,'xbmi set nx to ',nx
       case("x")
         ! process 0 is doing nothing, will not come here
         ! process 1 will receive total array
         ! this array will be distributed among p 1..isize-1
         if (irank .eq. 1) then
            call c_f_pointer(xptr, xpointer, [nx])
            ! xpointer now contains the x-values
            print *,irank,'xbmi get_var first and last elements',xpointer(1),xpointer(nx)
         else
            xpointer => dummyx
         endif

         allocate(displs(isize-1))
         allocate(sendcounts(isize-1))

         d = 0
         do i=1,isize-1
            call decomp(nx, isize-1, i-1, s, e)
            print *,irank,'xbmi nx, isize-1, i,s,e:',nx,isize-1,1,s,e
            sendcounts(i) = e - s + 1
            displs(i) = d
            d = d + sendcounts(i)
         enddo

         print *,'xbmi set_var displs:',displs
         print *,'xbmi set_var sendcounts',sendcounts

         ! compute start and end indices (s,e) in xpointer to
         ! copy to x

         call decomp(nx, isize-1, irank-1, s, e)
         if (allocated(x)) deallocate(x)
         allocate(x(e-s+1))
         print *,irank,'xbmi-1 nx, isize-1, i,s,e:',nx,isize-1,1,s,e,size(x)

         call mpi_scatterv(xpointer,sendcounts,displs,mpi_double_precision,x,e-s+1, &
         &                 mpi_double_precision,0,commwork,ierror)

         print *,irank,'xbmi get_var-1 first and last elements',x(1),x(size(x))

      end select
      ier = 0
   end function set_var

   subroutine get_var_type(c_var_name, c_type_name)  bind(C, name="get_var_type")
      character(kind=c_char), intent(in) :: c_var_name(*)
      character(kind=c_char), intent(out):: c_type_name(MAXSTRINGLEN)

      character(len=strlen(c_var_name)) :: var_name
      character(len=MAXSTRINGLEN)       :: type_name

      var_name = char_array_to_string(c_var_name)

      print *,irank,'xbmi get_var_type var_name = ',var_name
      select case(var_name)
       case("x")
         type_name = "double"
       case("nx")
         type_name = "int"
       case("comm")
         type_name = "int"
      end select

      c_type_name = string_to_char_array(trim(type_name))

   end subroutine get_var_type

   subroutine get_var_rank(c_var_name, rank) bind(C, name="get_var_rank")
      character(kind=c_char), intent(in) :: c_var_name(*)
      integer(c_int), intent(out)        :: rank

      ! The fortran name of the attribute name
      character(len=strlen(c_var_name)) :: var_name

      ! Store the name
      var_name = char_array_to_string(c_var_name)
      print *,irank,'xbmi get_var_rank var_name = ',var_name
      select case (var_name)
       case("x")
         rank = 1
       case("comm")
         rank = 0
       case("nx")
         rank  = 0
      end select

   end subroutine get_var_rank

   subroutine get_var_shape(c_var_name, shape) bind(C, name="get_var_shape")

      character(kind=c_char), intent(in) :: c_var_name(*)
      integer(c_int), intent(inout)      :: shape(6)

      character(len=strlen(c_var_name)) :: var_name

      var_name = char_array_to_string(c_var_name)
      print *,irank,'xbmi get_var_shape var_name = ',var_name
      shape = (/0, 0, 0, 0, 0, 0/)

      select case(var_name)
       case("x")
         if (irank.eq.1) then
            shape(1) = nx
         else
            shape(1) = 1
         endif
       case("comm")
         shape(1)=0
       case("nx")
         shape(1)=0
      end select

   end subroutine get_var_shape

   subroutine get_var(c_var_name, xptr) bind(C, name="get_var")

      ! Return a pointer to the variable

      character(kind=c_char), intent(in) :: c_var_name(*)
      type(c_ptr), intent(inout) :: xptr

      real*8, allocatable, target, save :: data(:)
      real*8, target, save :: dummydata(1)
      integer, allocatable :: displs(:), recvcounts(:)
      integer i,s,e,d,ierror
      character(len=strlen(c_var_name)) :: var_name

      var_name = char_array_to_string(c_var_name)
      print *,irank,'xbmi get_var var_name = ',var_name

      select case(var_name)
       case("comm")
         xptr = c_loc(comm)
       case("nx")
         xptr = c_loc(nx)
       case("x")
         if (allocated(data)) then
            deallocate(data)
         endif
         if (irank .eq. 1) then
            allocate(data(nx))
         else
            allocate(data(0))
         endif
         allocate(displs(isize-1))
         allocate(recvcounts(isize-1))
         d = 0
         do i=1,isize-1
            call decomp(nx, isize-1, i-1, s, e)
            recvcounts(i) = e - s + 1
            displs(i) = d
            d = d + recvcounts(i)
         enddo
         print *,'xbmi get_var displs:',displs
         print *,'xbmi get_var recvcounts',recvcounts
         call mpi_gatherv(x, size(x), mpi_double_precision, data, recvcounts, &
         &                displs,mpi_double_precision, 0, commwork, ierror)

         if(irank .eq. 1) then
            xptr = c_loc(data)
         else
            xptr = c_loc(dummydata)
         endif
      end select

   end subroutine get_var


   integer function finale()
      integer buf, ierror
      if (irank .eq. 1) then
         print *,irank,'xbmi sending message ...'
         call mpi_send(buf, 1, mpi_integer, 0, 20, comm, ierror)
         print *,irank,'xbmi sent message ...'
      endif
      if (.not. mpi_ready) then
         call mpi_finalize(ierror)
      endif
      finale = 0
   end function finale

   integer(c_int) function finalize() result(ierr) bind(C, name="finalize")
      print *,irank,'xbmi in finalize'
      ierr = finale()
   end function finalize

   integer(c_int) function initialize(c_configfile) result(ierr) bind(C, name="initialize")

      implicit none

      ! Variables
      character(kind=c_char), intent(in) :: c_configfile(*)

      if (len(c_configfile) .gt. -1) then
         ierr = init()
      endif

   end function initialize

   integer function init()
      integer ierror, color
      call mpi_initialized(mpi_ready,ierror)
      if (.not. mpi_ready) then
         call mpi_init(ierror)
      endif
      call mpi_comm_size(comm, isize, ierror)
      call mpi_comm_rank(comm, irank, ierror)

      print *,'bmi.f90 init: size:',isize,'rank:',irank

      ! create communicator commwork containing all processes minus process 0
      if (irank .eq. 0) then
         color = 0
      else
         color = 1
      endif

      call mpi_comm_split(comm, color, 0, commwork, ierror)

      if (irank .eq. 0) then
         commwork = mpi_comm_null
      endif
      ! process 0 will do nothing, except hanging in recv from process 1
      ! upon receipt of this message, process 0 will finalize
      if (irank .eq. 0) then
         call waitformessage
      endif

      init = 0

   end function init

   subroutine waitformessage
      integer buf, ierror
      print *,irank,'xbmi waiting for message ...'
      call mpi_recv(buf, 1, mpi_integer, 1, 20, comm, mpi_status_ignore, ierror)
      print *,irank,'xbmi got message'
      if (.not. mpi_ready) then
         call mpi_finalize(ierror)
      endif
   end subroutine waitformessage

end module bmimodule