c
       program mwave301
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      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
c      common /para_info/ks,ks1,kgs,ke,ke1,kge,kss,irank,isize
       common /para_x/is,is1,igs,ie,ie1,ige,iss,irankx,isizex
       common /para_y/js,js1,jgs,je,je1,jge,jss,iranky,isizey
       common /para_z/ks,ks1,kgs,ke,ke1,kge,kss,irankz,isizez
c      common /para_a/irank,isize,npex,npey,npez
       common /para_a/irank,isize
c for mpi_gatherv
       parameter (npex=2,npey=2,npez=2)
       parameter (npe=npex*npey*npez,npexy=npex*npey)
       integer recvcount(npe),displs(npe),bound(3,npe)
       integer itable(-1:npex,-1:npey,-1:npez)
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)
c      parameter(nb=4,iip0=8,iiq0=1,iir0=1,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)/npez+1)
       parameter(nyy=(ny2-1)/npey+1)
       parameter(nxx=(nx2-1)/npex+1)
       parameter(nxx3=nxx+2,nyy3=nyy+2,nzz3=nzz+2)
c      parameter(nwxy=nxx3*nyy3,nwyz=nyy3*nzz3,nwzx=nzz3*nxx3)
c
       dimension f(nb,0:nxx+1,0:nyy+1,0:nzz+1),
     1           u(nb,0:nxx+1,0:nyy+1,0:nzz+1),
     2           v(nb,0:nxx+1,0:nyy+1,0:nzz+1)
       dimension ftemp1x(nb,nyy3,nzz3),ftemp2x(nb,nyy3,nzz3)
       dimension ftemp1y(nb,nxx3,nzz3),ftemp2y(nb,nxx3,nzz3)
       dimension ftemp1z(nb,nxx3,nyy3),ftemp2z(nb,nxx3,nyy3)
       dimension fgl(nxx,nyy,nzz)
c      dimension utemp1x(nyy3,nzz3),utemp2x(nyy3,nzz3)
c      dimension utemp1y(nxx3,nzz3),utemp2y(nxx3,nzz3)
c      dimension utemp1z(nxx3,nyy3),utemp2z(nxx3,nyy3)
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
c
       isizex=npex
       isizey=npey
       isizez=npez
c
       kk=nz2/isizez
       kmod=mod(nz2,isizez)
       jj=ny2/isizey
       jmod=mod(ny2,isizey)
       ii=nx2/isizex
       imod=mod(nx2,isizex)
c
c
c      ii1=irank+1
c      irankz=(ii1-1)/npexy
c      ii2=ii1-irankz*npexy
c      iranky=(ii2-1)/npex
c      ii3=ii2-iranky*npex
c      irankx=ii3-1
c
       do k=-1,isizez
       do j=-1,isizey
       do i=-1,isizex
       itable(i,j,k)=MPI_PROC_NULL
       end do
       end do
       end do
c
       irank1=0
       do k=0,isizez-1
       do j=0,isizey-1
       do i=0,isizex-1
       itable(i,j,k)=irank1
       if(irank.eq.irank1) then
       irankx=i
       iranky=j
       irankz=k
       end if
       irank1=irank1+1
       end do
       end do
       end do
c
       ii3=irankx+1
       ii2=ii3+iranky*isizex
       ii1=ii2+irankz*isizey
c
c
       ks=1
       kss=irankz*kk+min(kmod,irankz)
       if (irankz.lt.kmod) kk=kk+1
       ke=ks+kk-1
       kgs=ks+kss
       kge=ke+kss
       ks1=ks
       ke1=ke
       if (irankz.eq.0) ks1=2
       if (irankz.eq.isizez-1) ke1=ke-1
c
       js=1
       jss=iranky*jj+min(jmod,iranky)
       if (iranky.lt.jmod) jj=jj+1
       je=js+jj-1
       jgs=js+jss
       jge=je+jss
       js1=js
       je1=je
       if (iranky.eq.0) js1=2
       if (iranky.eq.isizey-1) je1=je-1
c
       is=1
       iss=irankx*ii+min(imod,irankx)
       if (irankx.lt.imod) ii=ii+1
       ie=is+ii-1
       igs=is+iss
       ige=ie+iss
       is1=is
       ie1=ie
       if (irankx.eq.0) is1=2
       if (irankx.eq.isizex-1) ie1=ie-1
c
c
       nword=(ie-is+1)*(je-js+1)*(ke-ks+1)
       nwxy=(ie-is+3)*(je-js+3)
       nwyz=(je-js+3)*(ke-ks+3)
       nwzx=(ke-ks+3)*(ie-is+3)
       nword=nword*nb
       nwxy=nwxy*nb
       nwyz=nwyz*nb
       nwzx=nwzx*nb
c
       recvcount(ii1)=nword
       bound(1,ii1)=nwxy
       bound(2,ii1)=nwyz
       bound(3,ii1)=nwzx
c
       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 iia=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 20 m=1,nb
c      do m=1,nb
c
CC MPI START
c
       irightx = itable(isizex-1,iranky,irankz)
       ileftx  = itable(0,iranky,irankz)
       if (irankx.eq.isizex-1) then
       ftemp1x=f(:,ie-1,:,:)
c      ftemp2x=f(:,is,:,:)
       call mpi_send(ftemp1x,nwyz,mpi_real,ileftx,
     &               110,mpi_comm_world,ier)
       elseif (irankx.eq.0) then
       call mpi_recv(ftemp2x,nwyz,mpi_real,irightx,
     &               110,mpi_comm_world,istatus,ier)
       f(:,is,:,:)=ftemp2x
       end if
c
       if (irankx.eq.0) then
       ftemp1x=f(:,is+1,:,:)
       call mpi_send(ftemp1x,nwyz,mpi_real,irightx,
     &               115,mpi_comm_world,ier)
       elseif (irankx.eq.isizex-1) then
       call mpi_recv(ftemp2x,nwyz,mpi_real,ileftx,
     &               115,mpi_comm_world,istatus,ier)
       f(:,ie,:,:)=ftemp2x
       end if
c
c
       irighty = itable(irankx,isizey-1,irankz)
       ilefty  = itable(irankx,0,irankz)
       if (iranky.eq.isizey-1) then
       ftemp1y=f(:,:,je-1,:)
c      ftemp2y=f(:,:,js,:)
       call mpi_send(ftemp1y,nwzx,mpi_real,ilefty,
     &               120,mpi_comm_world,ier)
       elseif (iranky.eq.0) then
       call mpi_recv(ftemp2y,nwzx,mpi_real,irighty,
     &               120,mpi_comm_world,istatus,ier)
       f(:,:,js,:)=ftemp2y
       end if
c
       if (iranky.eq.0) then
       ftemp1y=f(:,:,js+1,:)
c      ftemp2y=f(:,:,je,:)
       call mpi_send(ftemp1y,nwzx,mpi_real,irighty,
     &               125,mpi_comm_world,ier)
       elseif (iranky.eq.isizey-1) then
       call mpi_recv(ftemp2y,nwzx,mpi_real,ilefty,
     &               125,mpi_comm_world,istatus,ier)
       f(:,:,je,:)=ftemp2y
       end if
c
c
       irightz = itable(irankx,iranky,isizez-1)
       ileftz  = itable(irankx,iranky,0)
       if (irankz.eq.isizez-1) then
       ftemp1z=f(:,:,:,ke-1)
c      ftemp2z=f(:,:,:,ks)
       call mpi_send(ftemp1z,nwxy,mpi_real,ileftz,
     &               130,mpi_comm_world,ier)
       elseif (irankz.eq.0) then
       call mpi_recv(ftemp2z,nwxy,mpi_real,irightz,
     &               130,mpi_comm_world,istatus,ier)
       f(:,:,:,ke+1)=ftemp2z
       end if
c
       if (irankz.eq.0) then
       ftemp1z=f(:,:,:,ks+1)
c      ftemp2z=f(:,:,:,ke)
       call mpi_send(ftemp1z,nwxy,mpi_real,irightz,
     &               135,mpi_comm_world,ier)
       elseif (irankz.eq.isizez-1) then
       call mpi_recv(ftemp2z,nwxy,mpi_real,ileftz,
     &               135,mpi_comm_world,istatus,ier)
       f(:,:,:,ke)=ftemp2z
       end if
c
c
CC MPI END
c
c  20  continue
c      end do
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
CC MPI START
c
       irightx = itable(irankx+1,iranky,irankz)
       ileftx  = itable(irankx-1,iranky,irankz)
c
       irighty = itable(irankx,iranky+1,irankz)
       ilefty  = itable(irankx,iranky-1,irankz)
c
       irightz = itable(irankx,iranky,irankz+1)
       ileftz  = itable(irankx,iranky,irankz-1)
c
c      do m=1,nb
       ftemp1x=f(:,is,:,:)
c      ftemp2x=f(:,ie+1,:,:)
       ftemp1y=f(:,:,js,:)
c      ftemp2y=f(:,:,je+1,:)
       ftemp1z=f(:,:,:,ks)
c      ftemp2z=f(:,:,:,ke+1)
       call mpi_sendrecv(ftemp1x,nwyz,mpi_real,ileftx,200,
     &                   ftemp2x,nwyz,mpi_real,irightx,200,
     &                   mpi_comm_world,istatus,ier)
       call mpi_sendrecv(ftemp1y,nwzx,mpi_real,ilefty,210,
     &                   ftemp2y,nwzx,mpi_real,irighty,210,
     &                   mpi_comm_world,istatus,ier)
       call mpi_sendrecv(ftemp1z,nwxy,mpi_real,ileftz,220,
     &                   ftemp2z,nwxy,mpi_real,irightz,220,
     &                   mpi_comm_world,istatus,ier)
c      f(:,is,:,:)=ftemp1x
       f(:,ie+1,:,:)=ftemp2x
c      f(:,:,js,:)=ftemp1y
       f(:,:,je+1,:)=ftemp2y
c      f(:,:,:,ks)=ftemp1z
       f(:,:,:,ke+1)=ftemp2z
c      end do
c
CC MPI END
c
CC MPI START
c      do 30 m=1,nb
c      do 32 k=1,nz1
c      do 32 k=ks,ke1
c      do 321 j=1,ny1
c      do 321 j=js,je1
c      do 321 i=1,nx1
c      do 321 i=is,ie1
       do k=ks,ke1
       do j=js,je1
       do i=is,ie1
       do m=1,nb
CC MPI END
       u(m,i,j,k)=0.125*(f(m,i,j,k)+f(m,i+1,j,k)+
     1            f(m,i,j+1,k)+f(m,i+1,j+1,k)+
     2            f(m,i,j,k+1)+f(m,i+1,j,k+1)+
     3            f(m,i,j+1,k+1)+f(m,i+1,j+1,k+1))
       end do
       end do
       end do
       end do
c 321  continue
c  32  continue
c  30  continue
c
       else
       t=t1
       dx=dx1
       dy=dy1
       dz=dz1
       endif
c
c
c      first step
c
CC MPI START
c      do 40 k=1,nz1
c      do 40 k=ks,ke1
c      do 401 j=1,ny1
c      do 401 j=js,je1
c      do 401 i=1,nx1
c      do 401 i=is,ie1
       do k=ks,ke1
       do j=js,je1
       do i=is,ie1
CC MPI END
c
       p(1)=0.125*(f(1,i,j,k)+f(1,i+1,j,k)+
     1      f(1,i,j+1,k)+f(1,i+1,j+1,k)+
     2      f(1,i,j,k+1)+f(1,i+1,j,k+1)+
     3      f(1,i,j+1,k+1)+f(1,i+1,j+1,k+1))
       p(2)=0.125*(f(2,i,j,k)+f(2,i+1,j,k)+
     1      f(2,i,j+1,k)+f(2,i+1,j+1,k)+
     2      f(2,i,j,k+1)+f(2,i+1,j,k+1)+
     3      f(2,i,j+1,k+1)+f(2,i+1,j+1,k+1))
       p(3)=0.125*(f(3,i,j,k)+f(3,i+1,j,k)+
     1      f(3,i,j+1,k)+f(3,i+1,j+1,k)+
     2      f(3,i,j,k+1)+f(3,i+1,j,k+1)+
     3      f(3,i,j+1,k+1)+f(3,i+1,j+1,k+1))
       p(4)=0.125*(f(4,i,j,k)+f(4,i+1,j,k)+
     1      f(4,i,j+1,k)+f(4,i+1,j+1,k)+
     2      f(4,i,j,k+1)+f(4,i+1,j,k+1)+
     3      f(4,i,j+1,k+1)+f(4,i+1,j+1,k+1))
c
       u(1,i,j,k)=u(1,i,j,k)
     1   -dx*(f(2,i+1,j+1,k+1)+f(2,i+1,j,k+1)
     2       -f(2,i,j+1,k+1)-f(2,i,j,k+1)
     3       +f(2,i+1,j+1,k)+f(2,i+1,j,k)
     4       -f(2,i,j+1,k)-f(2,i,j,k))
     5   -dy*(f(3,i+1,j+1,k+1)-f(3,i+1,j,k+1)
     6       +f(3,i,j+1,k+1)-f(3,i,j,k+1)
     7       +f(3,i+1,j+1,k)-f(3,i+1,j,k)
     8       +f(3,i,j+1,k)-f(3,i,j,k))
     1   -dz*(f(4,i+1,j+1,k+1)+f(4,i+1,j,k+1)
     2       +f(4,i,j+1,k+1)+f(4,i,j,k+1)
     3       -f(4,i+1,j+1,k)-f(4,i+1,j,k)
     4       -f(4,i,j+1,k)-f(4,i,j,k))
       u(2,i,j,k)=u(2,i,j,k)
     1   -dx*(f(1,i+1,j+1,k+1)+f(1,i+1,j,k+1)
     2       -f(1,i,j+1,k+1)-f(1,i,j,k+1)
     3       +f(1,i+1,j+1,k)+f(1,i+1,j,k)
     4       -f(1,i,j+1,k)-f(1,i,j,k))
       u(3,i,j,k)=u(3,i,j,k)
     1   -dy*(f(1,i+1,j+1,k+1)-f(1,i+1,j,k+1)
     2       +f(1,i,j+1,k+1)-f(1,i,j,k+1)
     3       +f(1,i+1,j+1,k)-f(1,i+1,j,k)
     4       +f(1,i,j+1,k)-f(1,i,j,k))
       u(4,i,j,k)=u(4,i,j,k)
     1   -dz*(f(1,i+1,j+1,k+1)+f(1,i+1,j,k+1)
     2       +f(1,i,j+1,k+1)+f(1,i,j,k+1)
     3       +f(1,i+1,j+1,k)+f(1,i+1,j,k)
     4       +f(1,i,j+1,k)+f(1,i,j,k))
c
       end do
       end do
       end do
c 401  continue
c  40  continue
c
c
c      preparation of second step
c      second step
c
c
c      do 61 m=1,nb
c
CC MPI START
c      do 62 k=1,nz2
c      do 62 k=ks,ke
c      do 621 j=1,ny2
c      do 621 j=js,je
c      do 621 i=1,nx2
c      do 621 i=is,ie
c      do 621 m=1,nb
       do k=ks,ke
       do j=js,je
       do i=is,ie
       do m=1,nb
CC MPI END
       v(m,i,j,k)=f(m,i,j,k)
       end do
       end do
       end do
       end do
c 621  continue
c  62  continue
c  61  continue
c
c
CC MPI START
c
       irightx = itable(irankx+1,iranky,irankz)
       ileftx  = itable(irankx-1,iranky,irankz)
c
       irighty = itable(irankx,iranky+1,irankz)
       ilefty  = itable(irankx,iranky-1,irankz)
c
       irightz = itable(irankx,iranky,irankz+1)
       ileftz  = itable(irankx,iranky,irankz-1)
c
c

c      do m=1,nb
c
       ftemp1x=u(:,ie,:,:)
       ftemp1y=u(:,:,je,:)
       ftemp1z=u(:,:,:,ke)
       call mpi_sendrecv(ftemp1x,nwyz,mpi_real,irightx,210,
     &                   ftemp2x,nwyz,mpi_real,ileftx,210,
     &                   mpi_comm_world,istatus,ier)
       call mpi_sendrecv(ftemp1y,nwzx,mpi_real,irighty,220,
     &                   ftemp2y,nwzx,mpi_real,ilefty,220,
     &                   mpi_comm_world,istatus,ier)
       call mpi_sendrecv(ftemp1z,nwxy,mpi_real,irightz,230,
     &                   ftemp2z,nwxy,mpi_real,ileftz,230,
     &                   mpi_comm_world,istatus,ier)
       u(:,is-1,:,:)=ftemp2x
       u(:,:,js-1,:)=ftemp2y
       u(:,:,:,ks-1)=ftemp2z
c      end do
c
       ftemp1x=v(:,ie,:,:)
       ftemp1y=v(:,:,je,:)
       ftemp1z=v(:,:,:,ke)
       call mpi_sendrecv(ftemp1x,nwyz,mpi_real,irightx,240,
     &                   ftemp2x,nwyz,mpi_real,ileftx,240,
     &                   mpi_comm_world,istatus,ier)
       call mpi_sendrecv(ftemp1y,nwzx,mpi_real,irighty,250,
     &                   ftemp2y,nwzx,mpi_real,ilefty,250,
     &                   mpi_comm_world,istatus,ier)
       call mpi_sendrecv(ftemp1z,nwxy,mpi_real,irightz,260,
     &                   ftemp2z,nwxy,mpi_real,ileftz,260,
     &                   mpi_comm_world,istatus,ier)
       v(:,is-1,:,:)=ftemp2x
       v(:,:,js-1,:)=ftemp2y
       v(:,:,:,ks-1)=ftemp2z
c
c
       ftemp1x=v(:,is,:,:)
       ftemp1y=v(:,:,js,:)
       ftemp1z=v(:,:,:,ks)
       call mpi_sendrecv(ftemp1x,nwyz,mpi_real,ileftx,270,
     &                   ftemp2x,nwyz,mpi_real,irightx,270,
     &                   mpi_comm_world,istatus,ier)
       call mpi_sendrecv(ftemp1y,nwzx,mpi_real,ilefty,280,
     &                   ftemp2y,nwzx,mpi_real,irighty,280,
     &                   mpi_comm_world,istatus,ier)
       call mpi_sendrecv(ftemp1z,nwxy,mpi_real,ileftz,290,
     &                   ftemp2z,nwxy,mpi_real,irightz,290,
     &                   mpi_comm_world,istatus,ier)
       v(:,ie+1,:,:)=ftemp2x
       v(:,:,je+1,:)=ftemp2y
       v(:,:,:,ke+1)=ftemp2z
c      end do
c
c      do i=1,npe
c      irank1=i-1
c      if(irank.eq.0) then
c      if(irank1.eq.irank) then
c      write(6,120) iia,i,irank,irankx,iranky,irankz,
c    *              isizex,isizey,isizez,irightx,ileftx,
c    *              irighty,ilefty,irightz,ileftz
c      end if
c      end do
  120  format(1h ,15i5)
CC MPI END
c
CC MPI START
c      do 60 k=2,nz1
c      do 60 k=ks1,ke1
c      do 601 j=2,ny1
c      do 601 j=js1,je1
c      do 601 i=2,nx1
c      do 601 i=is1,ie1
       do k=ks1,ke1
       do j=js1,je1
       do i=is1,ie1
CC MPI END
c
       p(1)=0.125*(u(1,i,j,k)+u(1,i-1,j,k)+
     1      u(1,i,j-1,k)+u(1,i-1,j-1,k)+
     2      u(1,i,j,k-1)+u(1,i-1,j,k-1)+
     3      u(1,i,j-1,k-1)+u(1,i-1,j-1,k-1))
       p(2)=0.125*(u(2,i,j,k)+u(2,i-1,j,k)+
     1      u(2,i,j-1,k)+u(2,i-1,j-1,k)+
     2      u(2,i,j,k-1)+u(2,i-1,j,k-1)+
     3      u(2,i,j-1,k-1)+u(2,i-1,j-1,k-1))
       p(3)=0.125*(u(3,i,j,k)+u(3,i-1,j,k)+
     1      u(3,i,j-1,k)+u(3,i-1,j-1,k)+
     2      u(3,i,j,k-1)+u(3,i-1,j,k-1)+
     3      u(3,i,j-1,k-1)+u(3,i-1,j-1,k-1))
       p(4)=0.125*(u(4,i,j,k)+u(4,i-1,j,k)+
     1      u(4,i,j-1,k)+u(4,i-1,j-1,k)+
     2      u(4,i,j,k-1)+u(4,i-1,j,k-1)+
     3      u(4,i,j-1,k-1)+u(4,i-1,j-1,k-1))
c
       f(1,i,j,k)=f(1,i,j,k)
     1   -dx1*(u(2,i,j,k)+u(2,i,j-1,k)
     2        -u(2,i-1,j,k)-u(2,i-1,j-1,k)
     3        +u(2,i,j,k-1)+u(2,i,j-1,k-1)
     4        -u(2,i-1,j,k-1)-u(2,i-1,j-1,k-1))
     5   -dy1*(u(3,i,j,k)-u(3,i,j-1,k)
     6        +u(3,i-1,j,k)-u(3,i-1,j-1,k)
     7        +u(3,i,j,k-1)-u(3,i,j-1,k-1)
     8        +u(3,i-1,j,k-1)-u(3,i-1,j-1,k-1))
     1   -dz1*(u(4,i,j,k)+u(4,i,j-1,k)
     2        +u(4,i-1,j,k)+u(4,i-1,j-1,k)
     3        -u(4,i,j,k-1)-u(4,i,j-1,k-1)
     4        -u(4,i-1,j,k-1)-u(4,i-1,j-1,k-1))
     5   +dx2*(v(1,i-1,j,k)-2.0*v(1,i,j,k)+v(1,i+1,j,k))
     6   +dy2*(v(1,i,j+1,k)-2.0*v(1,i,j,k)+v(1,i,j-1,k))
     7   +dz2*(v(1,i,j,k+1)-2.0*v(1,i,j,k)+v(1,i,j,k-1))
       f(2,i,j,k)=f(2,i,j,k)
     1   -dx1*(u(1,i,j,k)+u(1,i,j-1,k)
     2        -u(1,i-1,j,k)-u(1,i-1,j-1,k)
     3        +u(1,i,j,k-1)+u(1,i,j-1,k-1)
     4        -u(1,i-1,j,k-1)-u(1,i-1,j-1,k-1))
     5   +dx2*(v(2,i-1,j,k)-2.0*v(2,i,j,k)+v(2,i+1,j,k))
     6   +dy2*(v(2,i,j+1,k)-2.0*v(2,i,j,k)+v(2,i,j-1,k))
     7   +dz2*(v(2,i,j,k+1)-2.0*v(2,i,j,k)+v(2,i,j,k-1))
       f(3,i,j,k)=f(3,i,j,k)
     1   -dy1*(u(1,i,j,k)-u(1,i,j-1,k)
     2        +u(1,i-1,j,k)-u(1,i-1,j-1,k)
     3        +u(1,i,j,k-1)-u(1,i,j-1,k-1)
     4        +u(1,i-1,j,k-1)-u(1,i-1,j-1,k-1))
     5   +dx2*(v(3,i-1,j,k)-2.0*v(3,i,j,k)+v(3,i+1,j,k))
     6   +dy2*(v(3,i,j+1,k)-2.0*v(3,i,j,k)+v(3,i,j-1,k))
     7   +dz2*(v(3,i,j,k+1)-2.0*v(3,i,j,k)+v(3,i,j,k-1))
       f(4,i,j,k)=f(4,i,j,k)
     1   -dz1*(u(1,i,j,k)+u(1,i,j-1,k)
     2        +u(1,i-1,j,k)+u(1,i-1,j-1,k)
     3        -u(1,i,j,k-1)-u(1,i,j-1,k-1)
     4        -u(1,i-1,j,k-1)-u(1,i-1,j-1,k-1))
     5   +dx2*(v(4,i-1,j,k)-2.0*v(4,i,j,k)+v(4,i+1,j,k))
     6   +dy2*(v(4,i,j+1,k)-2.0*v(4,i,j,k)+v(4,i,j-1,k))
     7   +dz2*(v(4,i,j,k+1)-2.0*v(4,i,j,k)+v(4,i,j,k-1))
c
       end do
       end do
       end do
c 601  continue
c  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) iia,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
       fgl(1:nxx,1:nyy,1:nzz)=f(m,1:nxx,1:nyy,1:nzz)
c
c      do k=ks,ke
c      do j=js,je
c      do i=is,ie
c      fgl(i,j,k)=f(m,i,j,k)
c      end do
c      end do
c      end do
c
       call mpi_barrier(mpi_comm_world,ier)
c        call mpi_gatherv(f(1,1,1,m),mword,mpi_real,fg(1,1,1),
c    1                    recvcount,displs,mpi_real,0,
c    2                    mpi_comm_world,ier)
         call mpi_gatherv(fgl(1,1,1),mword,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)
CC MPI START
       include 'mpif.h'
       integer istatus(mpi_status_size)
c      common /para_info/ks,ks1,ke,ke1,kss,irank,isize
c      common /para_info/ks,ks1,kgs,ke,ke1,kge,kss,irank,isize
       common /para_x/is,is1,igs,ie,ie1,ige,iss,irankx,isizex
       common /para_y/js,js1,jgs,je,je1,jge,jss,iranky,isizey
       common /para_z/ks,ks1,kgs,ke,ke1,kge,kss,irankz,isizez
c      common /para_a/irank,isize,npex,npey,npez
       common /para_a/irank,isize
c for mpi_gatherv
       parameter (npex=2,npey=2,npez=2)
       parameter (npe=npex*npey*npez,npexy=npex*npey)
       integer recvcount(npe),displs(npe),bound(3,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)
       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)/npez+1)
       parameter(nyy=(ny2-1)/npey+1)
       parameter(nxx=(nx2-1)/npex+1)
c
c      dimension f(nb,0:nxx+1,0:nyy+1,0:nzz+1),
c    1           u(nb,0:nxx+1,0:nyy+1,0:nzz+1),
c    2           v(nb,0:nxx+1,0:nyy+1,0:nzz+1)
       dimension f(nb,0:nxx+1,0:nyy+1,0:nzz+1)
c for all_gather
c      dimension fg(nx2,ny2,nz2)
CC MPI END
c      dimension gf(nx2,ny2,nz2,nb)
       dimension ppin(10)
c      equivalence (gf,f)
c      common /blk/f
c
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)
c
       hx=xl/float(nx1)
       hy=yl/float(ny1)
       hz=zl/float(nz1)
c
CC MPI START
c      do 10 k=1,nz2
c      do 10 k=ks,ke
       do k=ks,ke
       z=0.5*hz*(2*(k+kss)-nz2-1)
c      do 101 j=1,ny2
c      do 101 j=js,je
       do j=js,je
       y=0.5*hy*(2*(j+jss)-ny2-1)
c      do 101 i=1,nx2
c      do 101 i=is,ie
       do i=is,ie
CC MPI END
c
c      z=0.5*hz*(2*k-nz2-1)
c      y=0.5*hy*(2*j-ny2-1)
       x=0.5*hx*(2*(i+iss)-nx2-1)
       dn1=0.0
       dv1=0.0
       ax1=sqrt(x*x+y*y+z*z)
       if(ax1.le.dxl) dn1=dn
       f(1,i,j,k)=dn1
       f(2,i,j,k)=0.0
       f(3,i,j,k)=0.0
       f(4,i,j,k)=0.0
       end do
       end do
       end do
c 101  continue
c  10  continue
c
       return
       end