	program laplace
	implicit double precision (a-h,o-z)

      parameter (mp=12) ! mp (number of processors) must be the same 
                        ! as in execution >mpiexec -np 12 laplace

	parameter (mn=201,mnode=mn*mn,mt=4)
c         mn = number of nodes in i and j directions (odd #)
c         mnode = total number of nodes
c         mt = number of different tasks
      common /t/t(mn,mn),told(mn,mn),nbc(mn,mn)

      dimension nod14m(mn,mn),nim(mnode),njm(mnode)
c         nod14m(mn,mn) = global node numbering
c         nim(mnode) = i direction location for each node
c         njm(mnode) = j direction location for each node
      dimension ntaskt(mt),ntaskp(mt),ntaskn(mt)
c         'ntaskt(i)' = total number of nodes for task i
c         'ntaskp(i)' = number of node calculations for task i in each processor
c         'ntaskn(i)' = number of proc's (task i) that calculate 'ntaskp(i)' nodes*
c         *remaining processors calculate 'ntaskp(i)-1' nodes
      parameter (mtp=(mn-1)/2+1)
      parameter (mtpm=mtp*mtp)
      dimension nod14t(mtpm,mt)
c         nod14t(mtpm,mt) = nodes on each task
c         mtp = maximum number of nodes in i/j directions in tasks
c         mtpm = maximum number of node calculations for any task

      dimension ntasks(mp,mt),ntaske(mp,mt)
c         ntasks(np,nt) = start counter* on processor np, in task nt
c         ntaske(np,nt) = end counter* on processor np, in task nt
c         *counter is n entry in nod14t(n,nt) for task nt

      parameter (mtrecv=mtpm/mp+1)
      dimension trecv(mtrecv)
c         trecv(i) = calculated values sent and received between processors
c         mtrecv = maximum number of nodes calculates on a given task, on a 
c                  given processor

      real e,etime,et(2)

c         MPI

	include 'mpif.h'
	integer status(MPI_STATUS_SIZE)
      call mpi_init(ierr)
	call mpi_comm_rank(MPI_COMM_WORLD,my_rank,ierr)
	call mpi_comm_size(MPI_COMM_WORLD,nsize,ierr)

      open(unit=1,file='t.dat',form='formatted',SHARED)

      if (my_rank.eq.0) then

c         Global node numbering

      nod=0
      do 2 j=1,mn
      do 2 i=1,mn
              nod=nod+1
              nod14m(i,j)=nod
              nim(nod)=i
              njm(nod)=j
2     continue

c         Tasks
c	The four tasks are based on the checkerboard
c	methodology. Any CFD node will eventually be
c	dependent (actual software developed for my work)
c	on all eight surrounding nodes (in 2D) so here 
c	instead of 2 tasks there are 4. In 3D there are more.

      npi=mn
      npj=mn
      nei=(mn-1)/2
      nej=(mn-1)/2

      do 3 nt=1,mt
      do 3 ntpm=1,mtpm
      nod14t(ntpm,nt)=0
3     continue

      nt=1 ! task 1
      n=0
      do 4 j=1,npj,2
      do 4 i=1,npi,2
      n=n+1
      nod14t(n,nt)=nod14m(i,j)
4     continue
      ntaskt(nt)=n

      nt=2 ! task 2
      n=0
      do 5 j=1,npj,2
      do 5 i=2,npi-1,2
      n=n+1
      nod14t(n,nt)=nod14m(i,j)
5     continue
      ntaskt(nt)=n

      nt=3 ! task 3
      n=0
      do 6 j=2,npj-1,2
      do 6 i=1,npi,2
      n=n+1
      nod14t(n,nt)=nod14m(i,j)
6     continue
      ntaskt(nt)=n

      nt=4 ! task 4
      n=0
      do 7 j=2,npj-1,2
      do 7 i=2,npi-1,2
      n=n+1
      nod14t(n,nt)=nod14m(i,j)
7     continue
      ntaskt(nt)=n

      do 8 nt=1,mt
      ntp=ntaskt(nt)/mp
      if (ntp*mp.eq.ntaskt(nt)) then
          ntaskp(nt)=ntp
          ntaskn(nt)=mp
      else
          ntaskp(nt)=ntp+1
          ntaskn(nt)=ntaskt(nt)-ntp*mp
      endif
8     continue

      do 11 nt=1,mt
      ncs=1
      ntn=ntaskn(nt)
      ntp=ntaskp(nt)
          do 9 np=1,ntn
              ntasks(np,nt)=ncs
              nce=ncs+ntp-1
              ntaske(np,nt)=nce
              ncs=nce+1
9         continue
              if (ntn.eq.mp) go to 11
              ntp=ntp-1
          do 10 np=ntn+1,mp
              ntasks(np,nt)=ncs
              nce=ncs+ntp-1
              ntaske(np,nt)=nce
              ncs=nce+1
10         continue
11    continue

      do 12 nt=1,mt
          if (ntaske(mp,nt).ne.ntaskt(nt)) then
              write(*,*) '# of nodes calculated in processors'
              write(*,*) 'do not correspond in task ',nt
              stop
          endif
12      continue

c         Initialize

      error=dlog10(1.0d-06)
      nc=1

	dx=1.0d0/dble(mn-1)
	do 15 i=1,mn
	do 15 j=1,mn
	t(i,j)=0.0d0
      nbc(i,j)=0
15    continue

c         Boundary Conditions

      ny=mn/10
	do 20 j=2,ny
      t(1,j)=dble(j-1)/dble(ny)
      nbc(1,j)=1
20    continue

      do 30 j=ny+1,mn
      t(1,j)=1.0d0
      nbc(1,j)=1
30    continue

      do 40 i=2,mn
      t(i,mn)=1.0d0
      nbc(i,mn)=1
40    continue

      do 41 i=1,mn
      nbc(i,1)=1
41    continue

      do 42 j=2,mn-1
      nbc(mn,j)=2
42    continue

c         to(i,j)

      do 45 i=1,mn
      do 45 j=1,mn
      told(i,j)=t(i,j)
45    continue

      e=etime(et)

      endif

c         Broadcast user input

      call MPI_Bcast(ntasks,mp*mt,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
      call MPI_Bcast(ntaske,mp*mt,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
      call MPI_Bcast(nod14t,mtpm*mt,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
      call MPI_Bcast(nim,mnode,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
      call MPI_Bcast(njm,mnode,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
      call MPI_Bcast(nbc,mnode,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
      call MPI_Bcast(error,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)

c         Laplace calculation - tasks divided on all processors

1     call MPI_Bcast(t,mnode,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
      do 70 nt=1,mt

      if (my_rank.eq.0) then
         np=mp
      else
         np=my_rank
      endif

      nts=ntasks(np,nt)
      nte=ntaske(np,nt)
      mc=0
      do 55 m=nts,nte
          mc=mc+1
          nod14=nod14t(m,nt)
          i=nim(nod14)
          j=njm(nod14)
              nbcv=nbc(i,j)
                  if (nbcv.eq.0) then
                      te=t(i+1,j)
                      tw=t(i-1,j)
                      tn=t(i,j+1)
                      ts=t(i,j-1)
                      trecvq=0.25d0*(te+tw+tn+ts)
                      go to 50
                  endif
                  if (nbcv.eq.1) then
                      trecvq=t(i,j)
                      go to 50
                  endif
                  if (nbcv.eq.2) then
                      trecvq=t(i-1,j)
                  endif
50    t(i,j)=trecvq
      if (my_rank.eq.0) go to 55
      trecv(mc)=trecvq
55    continue
      if (my_rank.ne.0) then
      call mpi_send(trecv,mtrecv,MPI_DOUBLE_PRECISION,0,1, ! send to master
     *		              MPI_COMM_WORLD,ierr)
      endif

      IF (my_rank.eq.0) THEN
          do 60 nproc=1,nsize-1
      call mpi_recv(trecv,mtrecv,MPI_DOUBLE_PRECISION,nproc,1, ! master receive
     *		              MPI_COMM_WORLD,status,ierr)
          nts=ntasks(nproc,nt)
          nte=ntaske(nproc,nt)
          mc=0
          do 60 m=nts,nte
          mc=mc+1
          nod14=nod14t(m,nt)
          i=nim(nod14)
          j=njm(nod14)
          t(i,j)=trecv(mc)
60        continue
      endif

70    continue

      if (my_rank.eq.0) then

      res=0.0d0
      do 92 i=2,mn
      do 92 j=2,mn-1
      tnewv=t(i,j)
      toldv=told(i,j)
      res=res+dabs(tnewv-toldv)
      told(i,j)=t(i,j)
92    continue

      if(res.eq.0.0d0) res=1.0d-99
      res=dlog10(res)

      write(*,94) nc,res
94    format(2x,i7,2x,f8.3)

      nc=nc+1

      if (res.le.error) then

      e=etime(et)

c         Output

      write(*,*) e,et(1),et(2)

      WRITE(1,*) mn,mn
      Y=-DX
      DO 110 J=1,MN
          Y=Y+DX
      X=-DX
      DO 110 I=1,MN
          X=X+DX
      WRITE(1,102) X,Y,T(I,J)
110   CONTINUE

102   FORMAT(3(2X,D15.8))

      endif

      endif

c         continue to next iteration...

      call MPI_Bcast(res,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
      if (res.gt.error) go to 1

c         ... else ... done.

      call mpi_finalize(ierr)

      end program laplace
