c
       program mwave2dimtest
c
c      MPI (Massage Passing Interface)
c
CC MPI START
       include 'mpif.h'
       integer istatus(mpi_status_size)
       parameter (npex=1,npey=2,npez=4)
       parameter (npe=npey*npez)
c
       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
       common /para_table/itable
       integer recvcount(npe),displs(npe),bound(3,npe)
       integer itable(-1:npey,-1:npez)
c
CC MPI  END
c
c
            parameter(nx=126,ny=126,nz=126)
c
            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)
c
            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
CC MPI  START
            parameter(nzz=(nz2-1)/npez+1)
            parameter(nyy=(ny2-1)/npey+1)
c
       parameter(nyy3=nyy+2,nzz3=nzz+2)
c
       dimension f(nx2,0:nyy+1,0:nzz+1,nb),u(nx2,0:nyy+1,0:nzz+1,nb),
     1           v(nx2,0:nyy+1,0:nzz+1,nb)
c
c      dimension ftemp1x(nyy3,nzz3,nb),ftemp2x(nyy3,nzz3,nb)
       dimension ftemp1y(nx2,nzz3,nb),ftemp2y(nx2,nzz3,nb)
       dimension ftemp1z(nx2,nyy3,nb),ftemp2z(nx2,nyy3,nb)
c
CC MPI END
            dimension p(nb),ppin(10)
            real*8 zt0,zt1,zt2,zt
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
c          isizex=npex
            isizey=npey
            isizez=npez
c
            kk=nz2/isizez
            kmod=mod(nz2,isizez)
            jj=ny2/isizey
            jmod=mod(ny2,isizey)
c          ii=nx2/isizex
c          imod=mod(nx2,isizex)
c
c
           do k=-1,isizez
           do j=-1,isizey
           itable(j,k)=MPI_PROC_NULL
           end do
           end do
c
           irank1=0
           do k=0,isizez-1
           do j=0,isizey-1
c         do i=0,isizex-1
           itable(j,k)=irank1
           if(irank.eq.irank1) then
c         irankx=i
           iranky=j
           irankz=k
           end if
           irank1=irank1+1
           end do
           end do
c
c
c
c       ii3=irankx+1
c      ii2=ii3+iranky*isizex
c       ii1=ii2+irankz*isizey
c
         ii2=iranky+1
         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
c         is=1
c         iss=irankx*ii+min(imod,irankx)
c         if (irankx.lt.imod) ii=ii+1
c         ie=is+ii-1
c         igs=is+iss
c         ige=ie+iss
c         is1=is
c         ie1=ie
c         if (irankx.eq.0) is1=2
c         if (irankx.eq.isizex-1) ie1=ie-1
c
c
           nword=nx2*(je-js+1)*(ke-ks+1)
           nwxy=nx2*(je-js+3)
c           nwyz=(je-js+3)*(ke-ks+3)
           nwzx=(ke-ks+3)*nx2
           nword=nword*nb
           nwxy=nwxy*nb
c           nwyz=nwyz*nb
           nwzx=nwzx*nb
c
       recvcount(ii1)=nword
       bound(1,ii1)=nwxy
c       bound(2,ii1)=nwyz
       bound(3,ii1)=nwzx
c
c
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
c            initial parameters
c
c
CC MPI START
c
          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
c
c
CC MPI START
c
c       irightx = itable(isizex-1,iranky,irankz)
c      ileftx  = itable(0,iranky,irankz)
c       if (irankx.eq.isizex-1) then
c       ftemp1x=f(:,ie-1,:,:)
c      ftemp2x=f(:,is,:,:)
c       call mpi_send(ftemp1x,nwyz,mpi_real,ileftx,
c     &               110,mpi_comm_world,ier)
c       elseif (irankx.eq.0) then
c       call mpi_recv(ftemp2x,nwyz,mpi_real,irightx,
c     &               110,mpi_comm_world,istatus,ier)
c       f(:,is,:,:)=ftemp2x
c       end if
c
c      if (irankx.eq.0) then
c       ftemp1x=f(:,is+1,:,:)
c       call mpi_send(ftemp1x,nwyz,mpi_real,irightx,
c     &               115,mpi_comm_world,ier)
c       elseif (irankx.eq.isizex-1) then
c       call mpi_recv(ftemp2x,nwyz,mpi_real,ileftx,
c     &               115,mpi_comm_world,istatus,ier)
c       f(:,ie,:,:)=ftemp2x
c       end if
c
c
       irighty = itable(isizey-1,irankz)
       ilefty  = itable(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(iranky,isizez-1)
       ileftz  = itable(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
       if(iip.eq.1) then
       t=0.5*t1
       dx=0.5*dx1
       dy=0.5*dy1
       dz=0.5*dz1
c
c
c
CC MPI START
c
c       irightx = itable(irankx+1,iranky,irankz)
c       ileftx  = itable(irankx-1,iranky,irankz)
c
         irighty = itable(iranky+1,irankz)
         ilefty  = itable(iranky-1,irankz)
c
         irightz = itable(iranky,irankz+1)
         ileftz  = itable(iranky,irankz-1)
c
c
c      do m=1,nb
c        ftemp1x=f(:,is,:,:)
c       ftemp2x=f(:,ie+1,:,:)
         ftemp1y=f(:,js,:,:)
c       ftemp2y=f(:,:,je+1,:)
         ftemp1z=f(:,:,ks,:)
c      ftemp2z=f(:,:,:,ke+1)
c
c        call mpi_sendrecv(ftemp1x,nwyz,mpi_real,ileftx,200,
c     &                   ftemp2x,nwyz,mpi_real,irightx,200,
c     &                   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
c       f(:,is,:,:)=ftemp1x
c       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
c
c
CC MPI START
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccc
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
c      do i=is,ie1
       do m=1,nb
CC MPI END
       u(m,i,j,k)=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))
       end do
       end do
       end do
c      end do
c 321  continue
c  32  continue
c  30  continue
c
       else
       t=t1
       dx=dx1
       dy=dy1
       dz=dz1
       endif
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
c        first step
c
c
CC MPI START
c
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
c        do i=is,ie1
c
CC MPI END
c
c
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
c         end do
           end do
           end do
c
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
c
CC MPI START
c
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
c
          do k=ks,ke
          do j=js,je
c        do i=is,ie
          do m=1,nb
CC MPI END
          v(i,j,k,m)=f(i,j,k,m)
          end do
          end do
          end do
c        end do
c  621  continue
c   62  continue
c   61  continue
c
c
CC MPI START
c
c       irightx = itable(irankx+1,iranky,irankz)
c       ileftx  = itable(irankx-1,iranky,irankz)
c
         irighty = itable(iranky+1,irankz)
         ilefty  = itable(iranky-1,irankz)
c
         irightz = itable(iranky,irankz+1)
         ileftz  = itable(iranky,irankz-1)
c
c
c      do m=1,nb
c
c       ftemp1x=u(:,ie,:,:)
       ftemp1y=u(:,je,:,:)
       ftemp1z=u(:,:,ke,:)
c       call mpi_sendrecv(ftemp1x,nwyz,mpi_real,irightx,210,
c     &                   ftemp2x,nwyz,mpi_real,ileftx,210,
c    &                   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)
c       u(:,is-1,:,:)=ftemp2x
       u(:,js-1,:,:)=ftemp2y
       u(:,:,ks-1,:)=ftemp2z
c      end do
c
c       ftemp1x=v(:,ie,:,:)
       ftemp1y=v(:,je,:,:)
       ftemp1z=v(:,:,ke,:)
c       call mpi_sendrecv(ftemp1x,nwyz,mpi_real,irightx,240,
c     &                   ftemp2x,nwyz,mpi_real,ileftx,240,
c     &                   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)
c       v(:,is-1,:,:)=ftemp2x
       v(:,js-1,:,:)=ftemp2y
       v(:,:,ks-1,:)=ftemp2z
c
c
c       ftemp1x=v(:,is,:,:)
       ftemp1y=v(:,js,:,:)
       ftemp1z=v(:,:,ks,:)
c       call mpi_sendrecv(ftemp1x,nwyz,mpi_real,ileftx,270,
c     &                   ftemp2x,nwyz,mpi_real,irightx,270,
c     &                   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)
c       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
c      do i=is1,ie1
CC MPI END
c
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
       end do
       end do
c       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
c
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
  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)
c
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
c       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=1,npey=2,npez=4)
       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)
c
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)
c      dimension f(nb,0:nxx+1,0:nyy+1,0:nzz+1)
        dimension f(nx2,0:nyy+1,0:nzz+1,nb)
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
c
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)
       do i=1,nx2
c      do 101 i=is,ie
c      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-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
       f(i,j,k,3)=0.0
       f(i,j,k,4)=0.0
       end do
       end do
       end do
c  10  continue
c
       return
       end