--- MITgcm_contrib/cg2d_bench/ini_mp.F	2006/05/12 21:58:05	1.1
+++ MITgcm_contrib/cg2d_bench/ini_mp.F	2006/05/12 22:24:08	1.2
@@ -1,3 +1,4 @@
+C	$Id: ini_mp.F,v 1.2 2006/05/12 22:24:08 ce107 Exp $	
       SUBROUTINE INI_MP
 
 C     Initialise multi-processing
@@ -8,17 +9,25 @@
 #include "EEPARAMS.h"
 
 #include "JAM_INFO.h"
-#include "MPI_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
-      INTEGER myTwoProcRank
+#ifdef DECOMP2D
+      integer dimens(2), coords(2), comm_cart
+      logical periods(2)
+      periods(1) = .true.
+      periods(2) = .true.
+#endif
 
       myXGlobalLo   =1
       myYGlobalLo   =1
@@ -36,14 +45,48 @@
 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
-
-      numberOfProcs = mpi_np
-      myProcId      = mpi_pid
+#endif
 
 #endif
 
@@ -64,7 +107,7 @@
       CALL GET_JAM_self_address(   jam_pid )
       CALL GET_JAM_partition_size( jam_np  )
       jam_exchKey = 100
-      dummyVal    = 0.
+      dummyVal    = 0. _d 0
       CALL JAM_barrier_start( dummyVal )
       CALL JAM_barrier_done(  dummyVal )
 
@@ -82,9 +125,15 @@
 
 #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)