C -*- Mode: Fortran; -*-
c-----------------------------------------------------------------
c     Ravi Samtaney
c     Copyright 2001
c     Princeton Plasma Physics Laboratory
c     All Rights Reserved
c-----------------------------------------------------------------
c     $Log: main.F,v $
c     Revision 1.3  2005/10/02 19:09:32  samtaney
c     Put in several options.
c
c     Revision 1.2  2005/09/23 19:06:43  samtaney
c     Added lines for easier port to parallel machines.
c
c     Revision 1.1.1.1  2005/09/12 17:28:10  samtaney
c     Original source.
c
c     Revision 1.2  2004/12/09 14:37:33  samtaney
c     Changed PERIODIC to ETAPERIODIC; periods(3) changed to .true.
c
c     Revision 1.1.1.1  2004/05/28 19:33:27  samtaney
c     Original source - copied from 3D Unsplit code.
c
c     Revision 1.1.1.1  2004/04/09 14:58:44  samtaney
c     Original unimesh Unsplit MHD solver.
c
c-----------------------------------------------------------------
c       Code to solve 3D ideal MHD equations
c-----------------------------------------------------------------
c
	program MHD
	use mesh_parms
	use mesh_common
	use iounits
#ifdef PARALLEL
	use mpistuff
#endif

	logical:: error
	integer:: num_nodes
	integer:: ipxnum
	integer:: ipynum
	integer:: ipznum
	integer:: ipx, ipy, ipz
	integer:: lpio

C       =======safety checks


	integer:: ndims
	logical:: periods(3), reorder
	integer:: NProcAxis(3)
c
#ifdef DYNAMIC
	namelist/meshsize/nx,ny,nz
	namelist/meshprocs/xprocs,yprocs,zprocs
c
c        open(imesh,file='mesh.inp',form='formatted')
c	read(imesh,meshsize)
c	read(imesh,meshprocs)
c        close(imesh)
	nx=16
	ny=16
	nz=16
	xprocs=1
	yprocs=1
	zprocs=1
#endif
C       ---- Initialize the communicator, call it Comm3D

	ndims = 3                        
#ifdef DYNAMIC
	NPROCS=XPROCS*YPROCS*ZPROCS
	nxlocal=nx/xprocs
	nylocal=ny/yprocs
	nzlocal=nz/zprocs
	nxlsize=nxlocal
	nylsize=nylocal
	nzlsize=nzlocal
#else
	nxlocal=nxlsize
	nylocal=nylsize
	nzlocal=nzlsize
#endif
#ifdef PARALLEL
	NProcAxis(1) = XPROCS            
	NProcAxis(2) = YPROCS            
	NProcAxis(3) = ZPROCS
	periods(1) = .false.             
#ifdef ETAPERIODIC
	periods(2) = .true.
#else
	periods(2) = .false.
#endif
	periods(3) = .true.
	reorder = .false.
	master = 0

	call MPI_init(ierr)          
	call MPI_cart_create( MPI_COMM_WORLD, ndims, NProcAxis, periods,
     1       reorder, comm3D, ierr)
	call MPI_comm_rank( comm3D, my_id, ierr)
	call MPI_comm_size( comm3D, num_nodes, ierr)
	
	write(6,*) 'Started proc #', my_id, num_nodes

        if (num_nodes .lt. NPROCS) then
           write(istderr,*) 'SORRY: PROGRAM ONLY WORKS FOR ', NPROCS,
     -          ' PROCESSES.'
	   call MPI_Finalize(ierr)
           stop
        endif

C	check for invalid combinations of dimensions and
C	process-grid dimensions

	error = .FALSE.

 	if (NX - (XPROCS-1)*NXlsize .le. 0) then
 		print*, 'invalid combination of NX, XPROCS'
 		error = .TRUE.
 	endif
 	if (NY - (YPROCS-1)*NYlsize .le. 0) then
 		print*, 'invalid combination of NY, YPROCS'
 		error = .TRUE.
 	endif
 	if (error) STOP

	if (my_id .lt. NPROCS) then
	   call pxpypz (my_id+1, ipx, ipy, ipz)
           call mesh_grid_proc(ipx,ipy,ipz)
        endif

	call MPI_Finalize(ierr)

#else
	ipx=1
	ipy=1
	ipz=1
        call mesh_grid_proc(ipx,ipy,ipz)
#endif

	end
C       ================================================================
C
C	main program for grid process
C


        subroutine mesh_grid_proc(a_iprocx, a_iprocy, a_iprocz)

      use mesh_parms
      use mesh_common

	integer:: a_iprocx, a_iprocy, a_iprocz
	integer, external:: iprocnum

	iprocx = a_iprocx
	iprocy = a_iprocy
	iprocz = a_iprocz

	iproc_idx = iprocnum(iprocx, iprocy, iprocz)

	NXlocal = NXlsize
	IF (iprocx .EQ. XPROCS) NXlocal = NX-(XPROCS-1)*NXlsize
	NYlocal = NYlsize
	IF (iprocy .EQ. YPROCS) NYlocal = NY-(YPROCS-1)*NYlsize
c
	call MHDMain
c
	return
	end