program s_bem
c
      include "mpif.h"
      integer my_rank,nproc,ierr
c     
c     Initialize MPI and determine individual processor and
c     the size of the world
c     
      call init(nproc,my_rank)
      if(my_rank.eq.0) then
         call master(nproc)
      else
         call slave(my_rank)
      endif
      call MPI_Finalize(ierr)
      stop
      end
c---------------------------------------------------------------
      subroutine init(nproc,my_rank)
      include "mpif.h"
c
c     initialize the MPI
c
      integer my_rank,nproc
      character*30 nodename
      integer nchar
      logical flag
      call MPI_Init(ierr)
      call MPI_Comm_rank(MPI_COMM_WORLD,my_rank,ierr)
      call MPI_Comm_size(MPI_COMM_WORLD,nproc,ierr)
      call MPI_Get_processor_name(nodename,nchar,ierr)
      write(*,*) 'nodename=',nodename
      return
      end
c
c
      subroutine master(nproc)
c ---------------------------------------------------------
c Example fortran program illustrating the use of PVM 3
c
c     n: max. number of stations
c     l: max. number of offsets at each station
c     m: max. number of panels
c ---------------------------------------------------------
c     program 3dy1 
      include "mpif.h"
      parameter(nt=51, m=51, l=500)
      real*8 x(l),y(l),z(l),xc(l),yc(l),zc(l),h1(l),g1(l)
      real*8 h(l,l),g(l,l),xn(l,3),u(l),q(l),area(l),ar
      real*8 xc1,yc1,zc1,pi,sum1,sum2
      integer nsta,np,ir(nt,m),npt(nt),icn(l,3),nem(nt)
      integer ne,icn1(500),icn2(500),icn3(500)
      integer i, nproc,iwho
      integer nhost
      parameter (BUFSIZE=50000)
      character buf(50000)
      call input1(nsta,npt,ir,x,y,z,np,nt,m,l)
      call cnnct(nsta,npt,nem,icn,ir,ne,nt,m,l)
      do i=1,ne
         icn1(i)=icn(i,1)
         icn2(i)=icn(i,2)
         icn3(i)=icn(i,3)
      enddo
      nhost=nproc-1
      do ik=1,nhost
         idest=ik
         n=ik
         call gs_send(nproc,ne,np,icn1,icn2,icn3,x,y,z,n,idest)
      enddo
      ik=nhost
      jk=1
      do while(jk.le.ne)
         call gs_recv(xc1,yc1,zc1,xn,ar,h1,g1,iwho,n)
         jk=jk+1
         write(*,*) 'result from proc:',iwho
         xc(n)=xc1
         yc(n)=yc1
         zc(n)=zc1
         area(n)=ar
         do j=1,ne
            h(n,j)=h1(j)
            g(n,j)=g1(j)
         enddo
         ik=ik+1
         if (ik.le.ne) then 
            write(*,*) 'sending data of column ik=',ik,' to processor ',
     1           iwho
            idest=iwho
            n=ik
            call gs_send(nproc,ne,np,icn1,icn2,icn3,x,y,z,n,idest)
            print *,'sent ',iwho
         endif
      enddo
c
      do ik=1,nhost
         idest=ik
         n=-1
         call gs_send(nproc,ne,np,icn1,icn2,icn3,x,y,z,n,idest)
      enddo     
      pi=dacos(-1.d0)
      do i=1,ne 
         sum1=0.0d0       
         sum2=0.0d0                                         
         do j=1,ne        
           sum1=sum1+h(i,j)                                   
           sum2=sum2+g(i,j)                                
           if(i.eq.j) h(i,j)=2.d0*pi+h(i,j)
         enddo
      enddo
      call bndry(ne,xc,yc,zc,xn,h,g,l,u,q,area)
      return
      end
c------------------------------------------------------------------
      subroutine gs_send(nproc,ne,np,icn1,icn2,icn3,x,y,z,n,idest)
      include "mpif.h"
      parameter(l=500)
      real*8 x(l),y(l),z(l),xc(l)
      real*8 xn(l,3)
      integer nproc,np,ne,icn1(500),icn2(500),icn3(500)
      parameter (BUFSIZE=50000)
      character buf(50000)
      itag1=10
      ipos=0
      call MPI_Pack(nproc,1,MPI_INTEGER,buf,BUFSIZE,ipos,
     *     MPI_COMM_WORLD,ierr)
      call MPI_Pack(ne,1,MPI_INTEGER,buf,BUFSIZE,ipos,
     *     MPI_COMM_WORLD,ierr)
      call MPI_Pack(np,1,MPI_INTEGER,buf,BUFSIZE,ipos,
     *     MPI_COMM_WORLD,ierr)
      call MPI_Pack(icn1,ne,MPI_INTEGER,buf,BUFSIZE,ipos,
     *     MPI_COMM_WORLD,ierr)
      call MPI_Pack(icn2,ne,MPI_INTEGER,buf,BUFSIZE,ipos,
     *     MPI_COMM_WORLD,ierr)
      call MPI_Pack(icn3,ne,MPI_INTEGER,buf,BUFSIZE,ipos,
     *     MPI_COMM_WORLD,ierr)
      call MPI_Pack(x,np,MPI_DOUBLE_PRECISION,buf,BUFSIZE,ipos,
     *     MPI_COMM_WORLD,ierr)
      call MPI_Pack(y,np,MPI_DOUBLE_PRECISION,buf,BUFSIZE,ipos,
     *     MPI_COMM_WORLD,ierr)
      call MPI_Pack(z,np,MPI_DOUBLE_PRECISION,buf,BUFSIZE,ipos,
     *     MPI_COMM_WORLD,ierr)
      call MPI_Pack(n,1,MPI_INTEGER,buf,BUFSIZE,ipos,
     *     MPI_COMM_WORLD,ierr)
      call MPI_Send(buf,ipos,MPI_PACKED,idest,itag1,
     *     MPI_COMM_WORLD,ierr)
      return
      end
c----------------------------------------------------------------
      subroutine gs_recv(xc1,yc1,zc1,xn,ar,h1,g1,iwho,n)
      include "mpif.h"
      parameter (l=500)
      real*8 xc1,yc1,zc1,xn(l,3),ar,h1(l),g1(l)
      integer iwho,n,status(MPI_STATUS_SIZE),nf
      parameter (BUFSIZE=50000)
      character buf(50000)
      ipos=0
      l3=3*l
      call MPI_Recv(buf,BUFSIZE,MPI_PACKED,MPI_ANY_SOURCE,
     1     MPI_ANY_TAG,MPI_COMM_WORLD,status,ierr)
      call MPI_Unpack(buf,BUFSIZE,ipos,xc1,1,MPI_DOUBLE_PRECISION,
     1     MPI_COMM_WORLD,ierr)
      call MPI_Unpack(buf,BUFSIZE,ipos,yc1,1,MPI_DOUBLE_PRECISION,
     1     MPI_COMM_WORLD,ierr)
      call MPI_Unpack(buf,BUFSIZE,ipos,zc1,1,MPI_DOUBLE_PRECISION,
     1     MPI_COMM_WORLD,ierr)
      call MPI_Unpack(buf,BUFSIZE,ipos,xn,l3,MPI_DOUBLE_PRECISION,
     1     MPI_COMM_WORLD,ierr)
      call MPI_Unpack(buf,BUFSIZE,ipos,ar,1,MPI_DOUBLE_PRECISION,
     1     MPI_COMM_WORLD,ierr)
      call MPI_Unpack(buf,BUFSIZE,ipos,nf,1,MPI_INTEGER,
     1     MPI_COMM_WORLD,ierr)
      call MPI_Unpack(buf,BUFSIZE,ipos,h1,nf,MPI_DOUBLE_PRECISION,
     1     MPI_COMM_WORLD,ierr)
      call MPI_Unpack(buf,BUFSIZE,ipos,g1,nf,MPI_DOUBLE_PRECISION,
     1     MPI_COMM_WORLD,ierr)
      call MPI_Unpack(buf,BUFSIZE,ipos,iwho,1,MPI_INTEGER,
     1     MPI_COMM_WORLD,ierr)
      call MPI_Unpack(buf,BUFSIZE,ipos,n,1,MPI_INTEGER,
     1     MPI_COMM_WORLD,ierr)
      return
      end