!***********************************************************
!
! Aug 2005: MPI version (Monte Carlo on N processors)
!
!************************************************************
      program MonteCarlo
      use mpi
!     Use the following include if the mpi module is not available
!     include "mpif.h"
      IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) 
      parameter ( nmax = 1024 , npars = 64 )
      dimension params(npars)
      dimension mypdfv(nmax,1000),mypdfs(nmax,100)
      dimension ipdfv(nmax,1000),ipdfs(nmax,100)
      data / iword / 8

      call MPI_INIT(ierr)
      call MPI_COMM_RANK(MPI_COMM_WORLD, myid, ierr)
      call MPI_COMM_SIZE(MPI_COMM_WORLD, numprocs, ierr)
!------------------------------------------------------------
! Read INPUT file for run parameters
!------------------------------------------------------------
      if ( myid .eq. 0) then
	open(3,file='INPUT',access='sequential',form='formatted')
	read(3,*)
	read(3,*)
	read(3,*)
	read(3,*)
	read(3,*)
	read(3,*) dt,tfinal,tstep1,tstep2,tstep3,tstep4,iseed0
	close(3)

	params(1) = dt
	params(2) = tfinal
	params(3) = tstep1
	params(4) = tstep2
	params(5) = tstep3
	params(6) = tstep4
      endif
!------------------------------------------------------------
! BroadCAST params from node 0
!------------------------------------------------------------
      call MPI_BCAST(params,npars,MPI_DOUBLE_PRECISION,0, &
	MPI_COMM_WORLD,ierr)

      if (myid .ne. 0 ) then
	dt     = params(1)
	tfinal = params(2)
	tstep1 = params(3)
	tstep2 = params(4)
	tstep3 = params(5)
	tstep4 = params(6)
      endif

      call MPI_BCAST(iseed0,1,MPI_INTEGER4,0, &
	MPI_COMM_WORLD,ierr)
!------------------------------------------------------------
!  Initialize OUTPUT files
!------------------------------------------------------------
      if (myid .eq. 0 ) then
	open(20,file='pdfv.dat',form='unformatted', &
		status='new',access='direct',recl=iword*nmax*1000)
	close(20)
	open(20,file='pdfs.dat',form='unformatted', &
		status='new',access='direct',recl=iword*nmax*100)
	close(20)
      endif
!------------------------------------------------------------
!  Set up problem locally
!------------------------------------------------------------
      iseed = myid + iseed0
      t = 0.0d0
!------------------------------------------------------------
      do iii=1,ntotal
!------------------------------------------------------------
	call rk2(v,nspike,ispike,tspike,t,dt,iseed)
	call computePDF(v,nspike,ispike,tspike,mypdfv,mypdfs)
!----------------------------------------finally advance TIME
	t = t + dt
!-------------------------------Checkpoint every 10000th step
	if (mod(iii,10000) .eq. 0) then
!------------------------------------------------------------
!  send data to node 0 and construct histogram/PDF
!------------------------------------------------------------
	  ncorr = 1000*nmax
	  call MPI_REDUCE(mypdfv,ipdfv,ncorr, &
		MPI_INTEGER4,MPI_SUM,0,MPI_COMM_WORLD,ierr)
	  ncorr = 100*nmax
	  call MPI_REDUCE(mypdfs,ipdfs,ncorr, &
		MPI_INTEGER4,MPI_SUM,0,MPI_COMM_WORLD,ierr)
!------------------------------------------------------------
	  if (myid .eq. 0) then
!------------------------------------------------------------
	    open(20,file='pdfv.dat',form='unformatted', &
		status='old',access='direct',recl=iword*nmax*1000)
	    write(20,rec=1) ipdfv
	    close(20)
	    open(20,file='pdfs.dat',form='unformatted', &
		status='old',access='direct',recl=iword*nmax*100)
	    write(20,rec=1) ipdfs
	    close(20)
!------------------------------------------------------------
	  endif
!-----------------------------------------------------------
	endif
!------------------------------------------------------------
      enddo
!------------------------------------End of TimeStepping Loop
      ncorr = 1000*nmax
      call MPI_REDUCE(mypdfv,ipdfv,ncorr, &
		MPI_INTEGER4,MPI_SUM,0,MPI_COMM_WORLD,ierr)
      ncorr = 100*nmax
      call MPI_REDUCE(mypdfs,ipdfs,ncorr, &
		MPI_INTEGER4,MPI_SUM,0,MPI_COMM_WORLD,ierr)
!------------------------------------------------------------
      if (myid .eq. 0) then
!------------------------------------------------------------
	open(20,file='pdfv.dat',form='unformatted', &
		status='old',access='direct',recl=iword*150/2)
	write(20,rec=1) ipdfv
	close(20)
	open(20,file='pdfs.dat',form='unformatted', &
		status='old',access='direct',recl=iword*nmax*100)
	write(20,rec=1) pdfs
	close(20)
      endif

      call MPI_FINALIZE(ierr)
!------------------------------------------------------------
 9999 stop
      end
!************************************************************