C	$Id: ini_mp.F,v 1.2 2006/05/12 22:24:08 ce107 Exp $	
      SUBROUTINE INI_MP

C     Initialise multi-processing
      IMPLICIT NONE

C     == Global data ==
#include "SIZE.h"
#include "EEPARAMS.h"

#include "JAM_INFO.h"

#ifdef ALLOW_MPI
#include "mpif.h"
#include "MPI_INFO.h"
#endif

C     == Local variables ==
#ifdef USE_JAM_INIT
      Real*8  dummyVal
      INTEGER myTwoProcRank
#endif
      INTEGER rc
      CHARACTER*(MAX_LEN_FNAM) fnam
#ifdef DECOMP2D
      integer dimens(2), coords(2), comm_cart
      logical periods(2)
      periods(1) = .true.
      periods(2) = .true.
#endif

      myXGlobalLo   =1
      myYGlobalLo   =1
      myProcId      =0
      numberOfProcs =1
      Nx            = sNx
      Ny            = sNy
      standardMessageUnit=6

#ifdef USE_MPI_INIT
C     MPI Initialisation
      CALL MPI_INIT( rc )
C     Get my proc. number
      CALL MPI_COMM_RANK( MPI_COMM_WORLD, mpi_pid, rc )
C     Get total count of procs.
      CALL MPI_COMM_SIZE( MPI_COMM_WORLD, mpi_np , rc )

      numberOfProcs = mpi_np
      myProcId      = mpi_pid
      comm_use      = MPI_COMM_WORLD
#ifdef DECOMP2D
c
       OPEN(UNIT=21,FILE="decomp.touse")
       read(21,*) dimens(1)
       read(21,*) dimens(2)
       if ((dimens(1)*dimens(2)) .ne. mpi_np) then
          if (mpi_pid .eq. 0) then
             write(0,*) 'Letting MPI choose a decomposition'
          endif
          dimens(1) = 0
          dimens(2) = 0
          call mpi_dims_create(mpi_np, 2, dimens, rc)
       endif

      call mpi_cart_create(MPI_COMM_WORLD, 2, dimens, periods, .true.,
     $     comm_cart, rc)
      if (mpi_pid .eq. 0) then
         write(0,*) 'using a ', dimens(1), ' x ', dimens(2),
     $        ' decomposition'
      endif
      call mpi_comm_rank(comm_cart, mpi_pid, rc)
      call mpi_cart_coords(comm_cart, mpi_pid, 2, coords, rc)
C     Set up connectivity
      call mpi_cart_shift(comm_cart, 0, 1, mpi_southId, mpi_northId, rc)
      call mpi_cart_shift(comm_cart, 1, 1, mpi_westId, mpi_eastId, rc)
c     create the derived datatype
      call mpi_type_vector(sNy,OLx,sNx+OLx*2,_MPI_TYPE_REAL,ewslice, rc)
      call mpi_type_commit(ewslice, rc)
      call mpi_type_vector(OLy,sNx,sNx+OLx*2,_MPI_TYPE_REAL,nsslice, rc) 
      call mpi_type_commit(nsslice, rc)
c     replace the communicator used
      comm_use = comm_cart
#else
C     Set up connectivity
      mpi_northId = mpi_pid+1
      IF ( mpi_northId .EQ. mpi_np ) mpi_northId = 0
      mpi_southId = mpi_pid-1
      IF ( mpi_southId .LT. 0      ) mpi_southId = mpi_np-1
#endif

#endif


#ifdef USE_JAM_INIT
C     JAM initialisation. This should work with or without
C     MPI. If we don't use MPI we have to start procs. by hand
C     though!

      myTwoProcRank=MOD(myProcId,2)
      IF ( myTwoProcRank .EQ. 0 ) THEN
       myTwoProcRank = 1
      ELSE
       myTwoProcRank = 0
      ENDIF
      CALL JAM_collective_init( myTwoProcRank )

      CALL GET_JAM_self_address(   jam_pid )
      CALL GET_JAM_partition_size( jam_np  )
      jam_exchKey = 100
      dummyVal    = 0. _d 0
      CALL JAM_barrier_start( dummyVal )
      CALL JAM_barrier_done(  dummyVal )

C     Set up connectivity
      jam_northId = jam_pid+1
      IF ( jam_northId .EQ. jam_np ) jam_northId = 0
      jam_southId = jam_pid-1
      IF ( jam_southId .LT. 0      ) jam_southId = jam_np-1

      numberOfProcs = jam_np
      myProcId      = jam_pid

      mpi_northId   = jam_northId
      mpi_southId   = jam_southId

#endif

#if defined(ALLOW_MPI) && defined(DECOMP2D)
      myXGlobalLo   = coords(1)*sNx+1
      myYGlobalLo   = coords(2)*sNy+1
      Nx            = dimens(1)*sNx
      Ny            = dimens(2)*sNy
#else
      myYGlobalLo   = myProcId*sNy+1
      Ny            = numberOfProcs*sNy
#endif
      IF ( numberOfProcs .GT. 1 ) THEN
       WRITE(fnam,'(A7,I6.6)') 'STDOUT.',myProcId
       OPEN(UNIT=standardMessageUnit,FILE=fnam)
      ENDIF


      RETURN
      END