c
       program mwave3
c
c      MPI (Massage Passing Interface)
cc     mwave3.f  main program of 3-dimensional wave equtaion
c      modified leap-frog method  2000.06.17
c      mwave701.f
c      implicit real*8 (a-h,o-z)
CC MPI START
       include 'mpif.h'
       integer istatus(mpi_status_size)
c      common /para_info/ks,ks1,ke,ke1,kss,irank,isize
       common /para_info/ks,ks1,kgs,ke,ke1,kge,kss,irank,isize
       common /para_a/nwxy
       common /para_table/itable
c for mpi_gatherv
       parameter (npe=8)
       integer recvcount(npe),displs(npe)
       integer itable(-1:npe)
CC MPI END
c      parameter(nx=30,ny=30,nz=30)
c      parameter(nx=100,ny=100,nz=100)
       parameter(nx=126,ny=126,nz=126)
c      parameter(nx=254,ny=254,nz=254)
c      parameter(nb=4,iip0=8,iiq0=4,iir0=4,last=4)
       parameter(nb=4,iip0=8,iiq0=1,iir0=1,last=4)
       parameter(nx1=nx+1,nx2=nx+2,ny1=ny+1,ny2=ny+2)
       parameter(nz1=nz+1,nz2=nz+2)
       parameter(n1=nx2,n2=n1*ny2,n3=n2*nz2,noinp=30)
       parameter(n4=n3*nb)
       parameter(thx=0.5,tam=0.5,vis0=0.10)
c
c      dimension f(nx2,ny2,nz2,nb),u(nx2,ny2,nz2,nb)
c      dimension gf(nx2,ny2,nz2,nb),gu(nx2,ny2,nz2,nb)
c      dimension v(nx2,ny2,nz2,nb)
CC MPI START
       parameter(nzz=(nz2-1)/npe+1)
       dimension f(nx2,ny2,0:nzz+1,nb),u(nx2,ny2,0:nzz+1,nb),
     1           v(nx2,ny2,0:nzz+1,nb)
c
       dimension ftemp1(nx2,ny2,nb),ftemp2(nx2,ny2,nb)
       dimension utemp1(nx2,ny2,nb),utemp2(nx2,ny2,nb)
       dimension vtemp1(nx2,ny2,nb),vtemp2(nx2,ny2,nb)
c
c for all_gather
       dimension fg(nx2,ny2,nz2)
CC MPI END
       dimension p(nb),ppin(10)
       real*8 zt0,zt1,zt2,zt
c      equivalence (gf,f),(gu,u)
c      common /blk/f
c
c           ppin xl   yl   zl   dxl  dyl  dzl  dn  dv   9  10
       data ppin/62.0,62.0,62.0,10.0,10.0,10.0,1.0,0.0, 9, 10/
c
CC MPI START
       call mpi_init(ier)
       call mpi_comm_rank(mpi_comm_world,irank,ier)
       call mpi_comm_size(mpi_comm_world,isize,ier)
c
       kk=nz2/isize
       kmod=mod(nz2,isize)
c
       do k=-1,isize
       itable(k)=MPI_PROC_NULL
       end do
c
       do k=0,isize-1
       itable(k)=k
       end do
c
       ks=1
       kss=irank*kk+min(kmod,irank)
c
       if (irank.lt.kmod) kk=kk+1
       ke=ks+kk-1
       kgs=ks+kss
       kge=ke+kss
       ks1=ks
       ke1=ke
       if (irank.eq.0) ks1=2
       if (irank.eq.isize-1) ke1=ke-1
c
       nword=(ke-ks+1)*nx2*ny2
       nwxy=n2*nb
       call mpi_gather(nword,1,mpi_integer,recvcount,
     *                 1,mpi_integer,0,mpi_comm_world,ier)
       displs(1)=0
       do i=2,isize
         displs(i)=displs(i-1)+recvcount(i-1)
       end do
c
CC MPI END
c
       xl=ppin(1)
       yl=ppin(2)
       zl=ppin(3)
       dxl=ppin(4)
       dyl=ppin(5)
       dzl=ppin(6)
       dn=ppin(7)
       dv=ppin(8)
       vis=vis0
c
       hx=xl/float(nx1)
       hy=yl/float(ny1)
       hz=zl/float(nz1)
       t1=thx*hx
       dx1=0.5*t1/hx
       dy1=0.5*t1/hy
       dz1=0.5*t1/hz
       dx2=vis*0.50*(t1/hx)**2
       dy2=vis*0.50*(t1/hy)**2
       dz2=vis*0.50*(t1/hz)**2
c
c      initial parameters
c
CC MPI START
       if (irank.eq.0) then
       write(6,122) nx,ny,nz,last,nx2,ny2,nz2
       write(6,124) t1,hx,hy,hz,dx1,dy1,dz1
       end if
CC MPI END
  122  format(1h ,10i8)
  124  format(1h ,8(1pe10.3))
c
c      initial conditions
       call ainit1(f,ppin)
c
c      start of calculation
       call clock(zt0)
c
       do 500 ii=1,last
       call clock(zt1)
       do 300 iir=1,iir0
       do 200 iiq=1,iiq0
       do 100 iip=1,iip0
c
c      boundary conditions
c
c      do 22 k=2,nz1
CC MPI START
       do 22 k=ks1,ke1
CC MPI END
c
       do 221 m=1,nb
       do 221 j=2,ny1
       f(1,j,k,m)=f(nx1,j,k,m)
       f(nx2,j,k,m)=f(2,j,k,m)
  221  continue
c
       do 222 m=1,nb
       do 222 i=1,nx2
       f(i,1,k,m)=f(i,ny1,k,m)
       f(i,ny2,k,m)=f(i,2,k,m)
  222  continue
   22  continue
c
CC MPI START
c
       iright = isize-1
       ileft  = 0
c
       if (irank.eq.isize-1) then
       ftemp1=f(:,:,ke-1,:)
       call mpi_send(ftemp1,nwxy,mpi_real,ileft,
     &               110,mpi_comm_world,ier)
       elseif (irank.eq.0) then
       call mpi_recv(ftemp2,nwxy,mpi_real,iright,
     &               110,mpi_comm_world,istatus,ier)
       f(:,:,ks,:)=ftemp2
       end if
c
       if (irank.eq.0) then
       ftemp1=f(:,:,ks+1,:)
       call mpi_send(ftemp1,nwxy,mpi_real,iright,
     &               115,mpi_comm_world,ier)
       elseif (irank.eq.isize-1) then
       call mpi_recv(ftemp2,nwxy,mpi_real,ileft,
     &               115,mpi_comm_world,istatus,ier)
       f(:,:,ke,:)=ftemp2
       end if
c
CC MPI END
c
   20  continue
c
c
c      case of iip=1  2-step Lax-Wendroff method
       if(iip.eq.1) then
       t=0.5*t1
       dx=0.5*dx1
       dy=0.5*dy1
       dz=0.5*dz1
c      do 30 m=1,nb
!hpf$  reflect f
CC MPI START
c
       iright = itable(irank + 1)
       ileft  = itable(irank - 1)
c
       ftemp1=f(:,:,ks,:)
c      ftemp2=f(:,:,ke+1,:)
       call mpi_sendrecv(ftemp1,nwxy,mpi_real,ileft,120,
     &                   ftemp2,nwxy,mpi_real,iright,120,
     &                   mpi_comm_world,istatus,ier)
c      f(:,:,ks,:)=ftemp1
       f(:,:,ke+1,:)=ftemp2
c
CC MPI END
c      do 32 k=1,nz1
CC MPI START
       do 30 m=1,nb
       do 32 k=ks,ke1
CC MPI END
       do 321 j=1,ny1
       do 321 i=1,nx1
       u(i,j,k,m)=0.125*(f(i,j,k,m)+f(i+1,j,k,m)+
     1            f(i,j+1,k,m)+f(i+1,j+1,k,m)+
     2            f(i,j,k+1,m)+f(i+1,j,k+1,m)+
     3            f(i,j+1,k+1,m)+f(i+1,j+1,k+1,m))
  321  continue
   32  continue
   30  continue
c
       else
       t=t1
       dx=dx1
       dy=dy1
       dz=dz1
       endif
c
c
c      first step
c
c      !hpf$  reflect f
c      do 40 k=1,nz1
CC MPI START
       do 40 k=ks,ke1
CC MPI END
       do 401 j=1,ny1
       do 401 i=1,nx1
c
       p(1)=0.125*(f(i,j,k,1)+f(i+1,j,k,1)+
     1      f(i,j+1,k,1)+f(i+1,j+1,k,1)+
     2      f(i,j,k+1,1)+f(i+1,j,k+1,1)+
     3      f(i,j+1,k+1,1)+f(i+1,j+1,k+1,1))
       p(2)=0.125*(f(i,j,k,2)+f(i+1,j,k,2)+
     1      f(i,j+1,k,2)+f(i+1,j+1,k,2)+
     2      f(i,j,k+1,2)+f(i+1,j,k+1,2)+
     3      f(i,j+1,k+1,2)+f(i+1,j+1,k+1,2))
       p(3)=0.125*(f(i,j,k,3)+f(i+1,j,k,3)+
     1      f(i,j+1,k,3)+f(i+1,j+1,k,3)+
     2      f(i,j,k+1,3)+f(i+1,j,k+1,3)+
     3      f(i,j+1,k+1,3)+f(i+1,j+1,k+1,3))
       p(4)=0.125*(f(i,j,k,4)+f(i+1,j,k,4)+
     1      f(i,j+1,k,4)+f(i+1,j+1,k,4)+
     2      f(i,j,k+1,4)+f(i+1,j,k+1,4)+
     3      f(i,j+1,k+1,4)+f(i+1,j+1,k+1,4))
c
       u(i,j,k,1)=u(i,j,k,1)
     1   -dx*(f(i+1,j+1,k+1,2)+f(i+1,j,k+1,2)
     2       -f(i,j+1,k+1,2)-f(i,j,k+1,2)
     3       +f(i+1,j+1,k,2)+f(i+1,j,k,2)
     4       -f(i,j+1,k,2)-f(i,j,k,2))
     5   -dy*(f(i+1,j+1,k+1,3)-f(i+1,j,k+1,3)
     6       +f(i,j+1,k+1,3)-f(i,j,k+1,3)
     7       +f(i+1,j+1,k,3)-f(i+1,j,k,3)
     8       +f(i,j+1,k,3)-f(i,j,k,3))
     1   -dz*(f(i+1,j+1,k+1,4)+f(i+1,j,k+1,4)
     2       +f(i,j+1,k+1,4)+f(i,j,k+1,4)
     3       -f(i+1,j+1,k,4)-f(i+1,j,k,4)
     4       -f(i,j+1,k,4)-f(i,j,k,4))
       u(i,j,k,2)=u(i,j,k,2)
     1   -dx*(f(i+1,j+1,k+1,1)+f(i+1,j,k+1,1)
     2       -f(i,j+1,k+1,1)-f(i,j,k+1,1)
     3       +f(i+1,j+1,k,1)+f(i+1,j,k,1)
     4       -f(i,j+1,k,1)-f(i,j,k,1))
       u(i,j,k,3)=u(i,j,k,3)
     1   -dy*(f(i+1,j+1,k+1,1)-f(i+1,j,k+1,1)
     2       +f(i,j+1,k+1,1)-f(i,j,k+1,1)
     3       +f(i+1,j+1,k,1)-f(i+1,j,k,1)
     4       +f(i,j+1,k,1)-f(i,j,k,1))
       u(i,j,k,4)=u(i,j,k,4)
     1   -dz*(f(i+1,j+1,k+1,1)+f(i+1,j,k+1,1)
     2       +f(i,j+1,k+1,1)+f(i,j,k+1,1)
     3       +f(i+1,j+1,k,1)+f(i+1,j,k,1)
     4       +f(i,j+1,k,1)+f(i,j,k,1))
c
  401  continue
   40  continue
c
c
c      preparation of second step
c      second step
c
c
c      do 61 m=1,nb
c
c      do 62 k=1,nz2
CC MPI START
       do 62 k=ks,ke
CC MPI END
       do 621 m=1,nb
       do 621 j=1,ny2
       do 621 i=1,nx2
       v(i,j,k,m)=f(i,j,k,m)
  621  continue
   62  continue
c  61  continue
c
!hpf$  reflect u
!hpf$  reflect v
CC MPI START
c
       iright = itable(irank + 1)
       ileft  = itable(irank - 1)
c
c      do m=1,nb
c
       utemp1=u(:,:,ke,:)
       call mpi_sendrecv(utemp1,nwxy,mpi_real,iright,130,
     &                   utemp2,nwxy,mpi_real,ileft,130,
     &                   mpi_comm_world,istatus,ier)
       u(:,:,ks-1,:)=utemp2
c
       utemp1=v(:,:,ke,:)
       call mpi_sendrecv(utemp1,nwxy,mpi_real,iright,140,
     &                   utemp2,nwxy,mpi_real,ileft,140,
     &                   mpi_comm_world,istatus,ier)
       v(:,:,ks-1,:)=utemp2
c
       utemp1=v(:,:,ks,:)
       call mpi_sendrecv(utemp1,nwxy,mpi_real,ileft,150,
     &                   utemp2,nwxy,mpi_real,iright,150,
     &                   mpi_comm_world,istatus,ier)
       v(:,:,ke+1,:)=utemp2
c
CC MPI END
c      do 60 k=2,nz1
CC MPI START
       do 60 k=ks1,ke1
CC MPI END
       do 601 j=2,ny1
       do 601 i=2,nx1
c
       p(1)=0.125*(u(i,j,k,1)+u(i-1,j,k,1)+
     1      u(i,j-1,k,1)+u(i-1,j-1,k,1)+
     2      u(i,j,k-1,1)+u(i-1,j,k-1,1)+
     3      u(i,j-1,k-1,1)+u(i-1,j-1,k-1,1))
       p(2)=0.125*(u(i,j,k,2)+u(i-1,j,k,2)+
     1      u(i,j-1,k,2)+u(i-1,j-1,k,2)+
     2      u(i,j,k-1,2)+u(i-1,j,k-1,2)+
     3      u(i,j-1,k-1,2)+u(i-1,j-1,k-1,2))
       p(3)=0.125*(u(i,j,k,3)+u(i-1,j,k,3)+
     1      u(i,j-1,k,3)+u(i-1,j-1,k,3)+
     2      u(i,j,k-1,3)+u(i-1,j,k-1,3)+
     3      u(i,j-1,k-1,3)+u(i-1,j-1,k-1,3))
       p(4)=0.125*(u(i,j,k,4)+u(i-1,j,k,4)+
     1      u(i,j-1,k,4)+u(i-1,j-1,k,4)+
     2      u(i,j,k-1,4)+u(i-1,j,k-1,4)+
     3      u(i,j-1,k-1,4)+u(i-1,j-1,k-1,4))
c
       f(i,j,k,1)=f(i,j,k,1)
     1   -dx1*(u(i,j,k,2)+u(i,j-1,k,2)
     2        -u(i-1,j,k,2)-u(i-1,j-1,k,2)
     3        +u(i,j,k-1,2)+u(i,j-1,k-1,2)
     4        -u(i-1,j,k-1,2)-u(i-1,j-1,k-1,2))
     5   -dy1*(u(i,j,k,3)-u(i,j-1,k,3)
     6        +u(i-1,j,k,3)-u(i-1,j-1,k,3)
     7        +u(i,j,k-1,3)-u(i,j-1,k-1,3)
     8        +u(i-1,j,k-1,3)-u(i-1,j-1,k-1,3))
     1   -dz1*(u(i,j,k,4)+u(i,j-1,k,4)
     2        +u(i-1,j,k,4)+u(i-1,j-1,k,4)
     3        -u(i,j,k-1,4)-u(i,j-1,k-1,4)
     4        -u(i-1,j,k-1,4)-u(i-1,j-1,k-1,4))
     5   +dx2*(v(i-1,j,k,1)-2.0*v(i,j,k,1)+v(i+1,j,k,1))
     6   +dy2*(v(i,j+1,k,1)-2.0*v(i,j,k,1)+v(i,j-1,k,1))
     7   +dz2*(v(i,j,k+1,1)-2.0*v(i,j,k,1)+v(i,j,k-1,1))
       f(i,j,k,2)=f(i,j,k,2)
     1   -dx1*(u(i,j,k,1)+u(i,j-1,k,1)
     2        -u(i-1,j,k,1)-u(i-1,j-1,k,1)
     3        +u(i,j,k-1,1)+u(i,j-1,k-1,1)
     4        -u(i-1,j,k-1,1)-u(i-1,j-1,k-1,1))
     5   +dx2*(v(i-1,j,k,2)-2.0*v(i,j,k,2)+v(i+1,j,k,2))
     6   +dy2*(v(i,j+1,k,2)-2.0*v(i,j,k,2)+v(i,j-1,k,2))
     7   +dz2*(v(i,j,k+1,2)-2.0*v(i,j,k,2)+v(i,j,k-1,2))
       f(i,j,k,3)=f(i,j,k,3)
     1   -dy1*(u(i,j,k,1)-u(i,j-1,k,1)
     2        +u(i-1,j,k,1)-u(i-1,j-1,k,1)
     3        +u(i,j,k-1,1)-u(i,j-1,k-1,1)
     4        +u(i-1,j,k-1,1)-u(i-1,j-1,k-1,1))
     5   +dx2*(v(i-1,j,k,3)-2.0*v(i,j,k,3)+v(i+1,j,k,3))
     6   +dy2*(v(i,j+1,k,3)-2.0*v(i,j,k,3)+v(i,j-1,k,3))
     7   +dz2*(v(i,j,k+1,3)-2.0*v(i,j,k,3)+v(i,j,k-1,3))
       f(i,j,k,4)=f(i,j,k,4)
     1   -dz1*(u(i,j,k,1)+u(i,j-1,k,1)
     2        +u(i-1,j,k,1)+u(i-1,j-1,k,1)
     3        -u(i,j,k-1,1)-u(i,j-1,k-1,1)
     4        -u(i-1,j,k-1,1)-u(i-1,j-1,k-1,1))
     5   +dx2*(v(i-1,j,k,4)-2.0*v(i,j,k,4)+v(i+1,j,k,4))
     6   +dy2*(v(i,j+1,k,4)-2.0*v(i,j,k,4)+v(i,j-1,k,4))
     7   +dz2*(v(i,j,k+1,4)-2.0*v(i,j,k,4)+v(i,j,k-1,4))
c
  601  continue
   60  continue
c
c      end of 1 time step advance
c
  100  continue
  200  continue
  300  continue
c
       call clock(zt2)
       zt1=zt1-zt0
       zt2=zt2-zt0
       zt=zt2-zt1
CC MPI START
       if (irank.eq.0)
     * write(6,402) ii,zt0,zt1,zt2,zt
CC MPI END
  402  format(1h , i6,1pe12.3,3(0pf12.5))
c
c      write the output data
c
c      write(ntap) f
       do 173 m=1,nb
CC MPI START
       call mpi_barrier(mpi_comm_world,ier)
         call mpi_gatherv(f(1,1,1,m),nword,mpi_real,fg(1,1,1),
     1                    recvcount,displs,mpi_real,0,
     2                    mpi_comm_world,ier)
c
       call mpi_barrier(mpi_comm_world,ier)
       if(irank.eq.0) then
c      do 1732 k=1,nz2
c      write(ntap) fg(1:nx2,1:ny2,k)
c1732  continue
c
       inz=nz2/2
       iny=ny2/2
       i1=1+n1*(iny-1)+n2*(inz-1)
       i2=nx2+n1*(iny-1)+n2*(inz-1)
       write(6,180) (fg(i,iny,inz),i=i1,i2)
  180  format(1h ,10(1pe10.2))
c
       end if
CC MPI END
  173  continue
c
c
  500  continue
    9  continue
c
c
       call mpi_finalize(ierr)
       stop
       end
       subroutine clock(ti)
       include 'mpif.h'
c      double precision mpi_wtime
       real*8 ti,ti1
c      call gettod(ti1)
c      call fjhpf_gettod(ti1)
c      call xclock(ti1,5)
       ti1=mpi_wtime()
c      ti=1.0d-6*ti1
       ti=ti1
       return
       end
       subroutine ainit1(f,ppin)
c      implicit real*8 (a-h,o-z)
       parameter(npe=8)
c      parameter(nx=30,ny=30,nz=30)
c      parameter(nx=100,ny=100,nz=100)
       parameter(nx=126,ny=126,nz=126)
c      parameter(nx=254,ny=254,nz=254)
       parameter(nb=4,iip0=8,iiq0=4,iir0=4,last=4)
       parameter(nx1=nx+1,nx2=nx+2,ny1=ny+1,ny2=ny+2)
       parameter(nz1=nz+1,nz2=nz+2)
       parameter(n1=nx2,n2=n1*ny2,n3=n2*nz2,noinp=30)
       parameter(n4=n3*nb)
c
c      dimension f(nx2,ny2,nz2,nb)
CC MPI START
       parameter(nzz=(nz2-1)/npe+1)
       dimension f(nx2,ny2,0:nzz+1,nb)
CC MPI END
c      dimension gf(nx2,ny2,nz2,nb)
       dimension ppin(10)
       integer itable(-1:npe)
c      equivalence (gf,f)
CC MPI START
       include 'mpif.h'
       integer istatus(mpi_status_size)
c      common /para_info/ks,ks1,ke,ke1,kss,irank,isize
       common /para_info/ks,ks1,kgs,ke,ke1,kge,kss,irank,isize
       common /para_a/nwxy
       common /para_table/itable
CC MPI END
c      common /blk/f
c
c
       pi=3.1415926
       xl=ppin(1)
       yl=ppin(2)
       zl=ppin(3)
       dxl=ppin(4)
       dyl=ppin(5)
       dzl=ppin(6)
       dn=ppin(7)
       dv=ppin(8)
c
       hx=xl/float(nx1)
       hy=yl/float(ny1)
       hz=zl/float(nz1)
c
c      do 10 k=1,nz2
CC MPI START
       do 10 k=ks,ke
CC MPI END
       z=0.5*hz*(2*(k+kss)-nz2-1)
       do 101 j=1,ny2
       y=0.5*hy*(2*j-ny2-1)
       do 101 i=1,nx2
c      z=0.5*hz*(2*k-nz2-1)
c      y=0.5*hy*(2*j-ny2-1)
       x=0.5*hx*(2*i-nx2-1)
       dn1=0.0
       dv1=0.0
       ax1=sqrt(x*x+y*y+z*z)
       if(ax1.le.dxl) dn1=dn
       f(i,j,k,1)=dn1
       f(i,j,k,2)=0.0
c
       f(i,j,k,2)=0.01*sin(2.0*2.0*pi*x/xl)
c
       f(i,j,k,3)=0.0
       f(i,j,k,4)=0.0
  101  continue
   10  continue
c
       return
       end