c
       program ogndd7
c
c      MPI (Massage Passing Interface)
c      file name  /zearthc2/earthc14.f
c      earthd14.f  modified leap-frog scheme
c      ogtac.cntl(ogncb77) from ognc.cntl(ogncd75)
c      3d mhd simulation of dipole tilt in 1/2 earth's magnetosphere
c      cartesian coordinate   finite resistivity  45 degree boudary
c      1990.10.22 modified  1986.06.16  by tatsuki ogino
c      implicit real*8 (a-h,o-z)
c
CC MPI START
       include 'mpif.h'
       integer istatus(mpi_status_size)
       parameter (npex=2,npey=2,npez=2)
       parameter (npe=npex*npey*npez,npexy=npex*npey)
c for mpi_gatherv
       integer recvcount(npe),displs(npe),bound(3,npe)
       integer itable(-1:npex,-1:npey,-1:npez)
       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
       common /para_table/itable
CC MPI END
c      parameter (nx=320,ny= 80,nz=160,nxp=100)
c      parameter (nx=500,ny=100,nz=200,nxp=150)
c      parameter (nx=318,ny=318,nz=318,nxp= 60)
c      parameter (nx=500,ny=200,nz=200,nxp=150)
       parameter (nx=510,ny=254,nz=254,nxp=155)
c      parameter (nx=500,ny=200,nz=400,nxp=150)
c      parameter (nx=640,ny=160,nz=320,nxp=200)
c      parameter (nx=800,ny=200,nz=400,nxp=250)
       parameter (n1=nx+2,n2=n1*(ny+2),n3=n2*(nz+2))
c      parameter (nb=8,nbb=11,n4=n3*nb,n5=n3*nbb,n6=n3*18)
       parameter (nb=8,mb=3,nbb=11)
c      parameter (last=2048,nil=2,lll=18,ibig=0)
c      parameter (last=8192,nil=2,lll=18,ibig=0)
c      parameter (last=4096,nil=2,lll=18,ibig=0)
c      parameter (last=1024,nil=2,lll=18,ibig=0)
       parameter (last=   8,nil=2,lll=18,ibig=0)
c      parameter (last=12288,nil=2,lll=18,ibig=0)
c      parameter (last= 448,nil=2,lll=18,ibig=0)
c      parameter (last= 32,nil=2,lll=18,ibig=0)
c      parameter (last= 32,nil=2,lll=18,ibig=0)
c      parameter (last=1920,nil=2,lll=18,ibig=0)
c      parameter (last=3840,nil=2,lll=18,ibig=0)
c      parameter (last=7680,nil=2,lll=18,ibig=0)
c      parameter (iiq0=8,iip0= 32,iis0=3840,thx=4.00)
c      parameter (iiq0=8,iip0= 32,iis0=3840,thx=3.50)
c      parameter (iiq0=8,iip0= 32,iis0=3840,thx=3.00)
c      parameter (iiq0=8,iip0= 32,iis0=1280,thx=2.00)
c      parameter (iiq0=8,iip0= 32,iis0=2560,thx=1.00)
c      parameter (iiq0=8,iip0= 32,iis0=5120,thx=0.50)
c      parameter (iiq0=8,iip0= 32,iis0=2048,thx=2.0833333)
c      parameter (iiq0=8,iip0= 32,iis0=8192,thx=0.5208333)
c      parameter (iiq0=8,iip0= 32,iis0=4096,thx=3.1266676)
c      parameter (iiq0=8,iip0= 32,iis0=1024,thx=2.0844451)
       parameter (iiq0=8,iip0=  8,iis0=   8,thx=2.0844451)
c      parameter (iiq0=8,iip0= 32,iis0=8192,thx=1.5633338)
c      parameter (rmuu=0.0020,eud0=0.0020,ro01=5.0e-4,pr01=3.56e-8)
       parameter (rmuu=0.0010,eud0=0.0000,ro01=5.0e-4,pr01=3.56e-8)
c      parameter (rmuu=0.0010,eud0=0.0010,ro01=5.0e-4,pr01=3.56e-8)
       parameter (aruu=30.0,pmu0=1.0,dfu0=1.0,xta=100.0,xtd=10.0)
       parameter (eatt=0.0010,rrat=0.1,rra1=0.1)
c      parameter (eatt=0.0020,rrat=0.2,rra1=0.2)
c      parameter (eatt=0.0020,rrat=0.1,rra1=0.1)
c      parameter (eatt=0.0020,rrat=0.3,rra1=0.3)
c      parameter (eatt=0.0010,rrat=0.2,rra1=0.2)
       parameter (nx1=nx+1,nx2=nx+2,ny1=ny+1,ny2=ny+2)
       parameter (nz1=nz+1,nz2=nz+2,nz3=nz+3)
       parameter (nin=1,itap=1)
       parameter (k2=n2*2,k3=n2*3,k4=n2*4,k5=n2*5,k6=n2*6)
       parameter (k7=n2*7,k8=n2*8,k9=n2*9,k10=n2*10,k11=n2*11)
c      parameter (orid= 45.375,ths0=0.0,kk8=k3*8)
c      parameter (orid= 90.375,ths0=0.0,kk8=k3*8)
c      parameter (orid=270.375,ths0=0.0,kk8=k3*8)
c      parameter (orid= 45.000,ths0=0.0,kk8=k3*8)
       parameter (orid= 90.000,ths0=0.0,kk8=k3*8)
c      parameter (orid=315.000,ths0=0.0,kk8=k3*8)
c      parameter (orid=  0.000,ths0=0.0,kk8=k3*8)
c      parameter (orid=270.1875,ths0=0.0,kk8=k3*8)
c      parameter (mfd=6,nfd=1441,iinob=last*12*0,ttno=128.0)
c      parameter (mfd=6,nfd=1441,iinob=last*12*2,ttno=256.0)
c      parameter (mfd=6,nfd=1441,iinob=last*24*4,ttno=256.0)
c      parameter (mfd=6,nfd=1441,iinob=last*6*19,ttno=256.0)
c      parameter (mfd=6,nfd=1441,iinob=last*6*20,ttno=512.0)
c      parameter (mfd=6,nfd=1441,iinob=last*780,ttno=512.0)
c      parameter (mfd=6,nfd=1441,iinob=last*0, ttno=102.4)
c      parameter (mfd=6,nfd=1441,iinob=last*0, ttno=204.8)
c      parameter (mfd=6,nfd=1441,iinob=last*25, ttno=819.2)
c      parameter (mfd=6,nfd=1441,iinob=last*140, ttno=409.6)
       parameter (mfd=6,nfd=1441,iinob=last*30, ttno=409.6)
c      parameter (mfd=6,nfd=1441,iinob=last*1410, ttno=448.0)
c      parameter (mfd=6,nfd=1441,iinob=last*30, ttno=448.0)
c      parameter (mfd=6,nfd=1441,iinob=last*40,ttno=102.4)
c      parameter (mfd=6,nfd=1441,iinob=last*12*1,ttno=512.0)
       double precision zt0,zt1,zt2,zt3,zt
       integer (kind=8):: n6,n5,n4
c
c      dimension f(nx2,ny2,nz2,nb),u(nx2,ny2,nz2,nb),
c    1           ff(nx2,ny2,nz2,nb),p(nx2,ny2,nz2,nbb),
c    2           pp(nx2,ny2,nz2,3),fdd(mfd,nfd)
c      dimension gf(nx2,ny2,nz2,nb),gu(nx2,ny2,nz2,nb),
c    1           gff(nx2,ny2,nz2,nb),gpp(nx2,ny2,nz2,3),
c    2           gfdd(mfd,nfd)
c      dimension gf(nx2,ny2,nz2,nb),gfdd(mfd,nfd)
c      dimension gf(nx2,ny2),gfdd(mfd,nfd)
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
       dimension f(nb,0:nxx+1,0:nyy+1,0:nzz+1),
     1           u(nb,0:nxx+1,0:nyy+1,0:nzz+1),
     1           ff(nb,0:nxx+1,0:nyy+1,0:nzz+1),
     1           p(nbb,0:nxx+1,0:nyy+1,0:nzz+1),
     2           pp(mb,0:nxx+1,0:nyy+1,0:nzz+1)
       dimension fgl(nxx,nyy,nzz)
       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 utemp1x(mb,nyy3,nzz3),utemp2x(mb,nyy3,nzz3)
       dimension utemp1y(mb,nxx3,nzz3),utemp2y(mb,nxx3,nzz3)
       dimension utemp1z(mb,nxx3,nyy3),utemp2z(mb,nxx3,nyy3)
c for all_gather
c      dimension fg(nx2,ny2,nz2,nb)
       dimension fg(nx2,ny2,nz2)
CC MPI END
       dimension gfdd(mfd,nfd)
       dimension hhx(n1),fbb(16),
     1           cj(10),cp(11),v(n2)
c
c      common /blk/f,pp
c      m  1ro 2vr 3vo 4vz 5pr 6br 7bo 8bz 9jr 10jo 11jz
c
c      cj      r1   dr  dz  b0  rm  a2    aid ar1 gam gra
       data cj/10.0,1.2,0.5,1.0,1.9,-0.01,0.7,0.8,0.9,1.0/
c      data cj/10.0,0.3,0.5,1.0,1.9,-0.01,0.7,0.8,0.9,1.0/
c
c      cp      xl    yl   zl   ra  vo0  p0   gra  til  vsw  ibr ibz
c      data cp/150.3,90.3,90.3,3.5,0.00,2.68,1.35, 30.,.044,0.0,-0.0/
c      data cp/150.3, 60.3, 60.3,3.5,0.0,2.68,1.35,30.,.044,0.0,-0.0/
       data cp/153.3, 76.5, 76.5,3.5,0.0,2.68,1.35,30.,.044,0.0,-0.0/
c      data cp/150.3,45.3,90.3,3.5,0.00,2.68,1.35,-30.,.044,0.0,-0.0/
c
c
CC MPI START
       call clock(zt0)
       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)
       mword=nword
       mwxy=nwxy*mb
       mwyz=nwyz*mb
       mwzx=nwzx*mb
c      nword=nword*nb
       nword=nword
       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
c for debug +++
c     write(6,*) '  irank  ,  ks  ,   ke  ',irank,ks,ke
c     if (irank.eq.0) write(6,*) ' value of recvcount ',
c    *                (recvcount(iia),iia=1,isize)
c     if (irank.eq.0) write(6,*) ' value of displs ',
c    *                (displs(iia),iia=1,isize)
c
CC MPI END
c
c      rewind 12
c
CC MPI START
       if(irank.eq.0) then
c      open(9,file='wn991021.data',access='sequential',
c    1         form='unformatted')
c      open(9,file='wn991022.data',access='sequential',
c    1         form='unformatted')
c      open(10,file='./data1/earthd/edd1d420.data',
c    1         access='sequential',form='unformatted')
c
c      open(11,file='ede3d301.data',
c    1         access='sequential',form='unformatted')
c      open(12,file='ede3d302.data',
c    1         access='sequential',form='unformatted')
c      open(13,file='ede3d303.data',
c    1         access='sequential',form='unformatted')
       end if
CC MPI END
c
c
       n4=n3*nb
       n5=n3*nbb
       n6=n3*18
c
c      read(9) gfdd
       do 660 i=1,nfd
c      read(9,662) (gfdd(j,i),j=1,mfd)
       do 664 j=1,mfd
c      fdd(j,i)=gfdd(j,i)
  664  continue
  662  format(1h ,6(1pe12.4))
  660  continue
c
c
c      til0=30.0
       til0=0.0
c      til0=90.0
       tild=5.0
       bis0=1.5       
c
c      do 300 iii=1,5
c      do 300 iii=1,10
       do 300 iii=1,3
c      do 300 iii=1,6
c
       ntap=10+iii
c      cp(11)=bisz*float(iii-2)
       cp(11)=bis0
c      cp(8)=til0+5.0*float(iii)
c      cp(8)=til0+10.0*float(iii)
       cp(8)=til0
c
       eat0=eatt
       rmu0=rmuu
       rmu0=rmuu
       aru=aruu
       pi=3.1415926
       gam=5.0/3.0
       gm1=2.0/(gam-1.0)
       gm2=0.5*(gam-1.0)
       cp(6)=10.0*(gam-1.0)*cp(7)/gam
       gra=cp(7)*1.0e-6
       po0=cp(6)*1.0e-7
       bis=cp(11)*1.0e-4
c      th=0.0+ths0*float(last)
c      th=orid*th*pi/180.0
       th=orid*pi/180.0
c      th=(orid+15.0*float(iii))*pi/180.0
       bisy=bis*cos(th)
       bisz=bis*sin(th)
       cp(11)=0.0
c      cp(11)=cp(11)
       vsw=cp(9)
       ar1=cp(4)
       ar2=ar1*ar1
       hx=cp(1)/float(nx1)
       hy=cp(2)/float(ny1)
       hz=cp(3)/float(nz1)
       t=0.5*hx*thx
       t1=0.5*t
       cj(8)=ar1
       cj(9)=gam
       cj(10)=gra
c
       dx1=0.25*t1/hx
       dy1=0.25*t1/hy
       dz1=0.25*t1/hz
       dx2=0.25*t/hx
       dy2=0.25*t/hy
       dz2=0.25*t/hz
       dx3=t1/(hx*hx)
       dy3=t1/(hy*hy)
       dz3=t1/(hz*hz)
       dx4=t/(hx*hx)
       dy4=t/(hy*hy)
       dz4=t/(hz*hz)
c
       rmu=t*rmu0*ro01
       pmu=t*rmu0*pmu0
       dfu=t*rmu0*dfu0
       eud=eud0
       ro02=rrat*ro01
       pr02=rrat*pr01
c
CC MPI START
       if(irank.eq.0)
     1 write (6,12) iii,last,nx,ny,nz,n1,n2,n3,n4,n5,n6,eat0,rmu0,aru,
     2 eud,rrat,hx,hy,hz,t,t1,ro01,pr01,gra,dx2,dy2,dz2,dx4,dy4,dz4,
     3 bis,orid,th,bisy,bisz,(cp(i),i=1,11),(cj(j),j=1,10)
CC MPI END
   12  format(1h ,5x,11i10/(1h ,5x,1p8e15.5))
c
       if(iii.ge.2) go to 400
c
c
c      do 20 k=1,nz2
       do 20 k=ks,ke
c      do 20 j=1,ny2
       do 20 j=js,je
c      do 20 i=1,nx2
       do 20 i=is,ie
       do 20 m=1,nb
       u(m,i,j,k)=0.0
       f(m,i,j,k)=0.0
   20  continue
c
       do 221 i=1,n1
  221  hhx(i)=hx*float(i)
       do 222 i=1,nb
       fbb(i+nb)=1.0
  222  fbb(i)=1.0
       fbb(4)=-1.0
       fbb(6)=-1.0
       fbb(7)=-1.0
       fbb(3+nb)=-1.0
       fbb(7+nb)=-1.0
c
c      initial condition
c      if(iii.eq.1) call equib7(f,pp,nxp,ro01,pr01,rra1,cj,cp)
       if(iii.eq.1) call equib8(f,pp,nxp,ro01,pr01,rra1,cj,cp)
       if(iii.ge.2) call equib8a(u,pp,nxp,ro01,pr01,rra1,cj,cp)
c      go to 882
c      write(6,212) (f(i),i=1,n4)
  212  format(1h ,2x,1p10e12.4)
c 882  continue
c
c
c      initial condition
c
c      read(10) f
c      if(i.eq.1) write(ntap) f
       if(nin.eq.0) go to 400
       if(iii.ge.2) go to 400
c
c
       do 410 iia=1,itap
       do 174 m=1,nb
CC MPI START
       if(irank.eq.0) then
c
       do 1742 k=1,nz2
c      read(10) fg(1:nx2,1:ny2,k)
 1742  continue
       end if
c      call mpi_scatterv(fg(1,1,1),recvcount,displs,
c    1                    mpi_real,f(1,1,1,m),mword,mpi_real,0,
c    2                    mpi_comm_world,ier)
c      call mpi_scatterv(fg(1,1,1),recvcount,displs,
c    1                    mpi_real,fgl(1,1,1),mword,mpi_real,0,
c    2                    mpi_comm_world,ier)
c
       call mpi_barrier(mpi_comm_world,ier)
c      f(m,1:nxx,1:nyy,1:nzz)=fgl(1:nxx,1:nyy,1:nzz)
c
c      do k=ks,ke
c      do j=js,je
c      do i=is,ie
c      f(m,i,j,k)=fgl(i,j,k)
c      end do
c      end do
c      end do
c
       call mpi_barrier(mpi_comm_world,ier)
CC MPI END
  174  continue
  410  continue
  400  continue
c
c
c
c      do 280 k=2,nz1
       do 280 k=ks1,ke1
c
       xx4=0.5*hx*float(2*nxp-nx1-2)
       z=0.5*hz*float(2*(kss+k)-nz2-1)
c      do 280 j=2,ny1
       do 280 j=js1,je1
c      y=0.5*hy*float(2*j-3)
       y=0.5*hy*float(2*(jss+j)-ny2-1)
c
c      do 280 i=2,nx1
       do 280 i=is1,ie1
       x=hhx(iss+i)+xx4
       ro2=x*x+y*y+z*z
       x2=amax1(ro2/ar2-1.0,0.0)
       x2=aru*x2*x2
       x2=x2/(x2+1.0)
       f(8,i,j,k)=f(8,i,j,k)*x2+u(8,i,j,k)*(1.0-x2)
       f(7,i,j,k)=f(7,i,j,k)*x2+u(7,i,j,k)*(1.0-x2)
       f(6,i,j,k)=f(6,i,j,k)*x2+u(6,i,j,k)*(1.0-x2)
       f(5,i,j,k)=f(5,i,j,k)*x2+u(5,i,j,k)*(1.0-x2)
       f(4,i,j,k)=f(4,i,j,k)*x2+u(4,i,j,k)*(1.0-x2)
       f(3,i,j,k)=f(3,i,j,k)*x2+u(3,i,j,k)*(1.0-x2)
       f(2,i,j,k)=f(2,i,j,k)*x2+u(2,i,j,k)*(1.0-x2)
       f(1,i,j,k)=f(1,i,j,k)*x2+u(1,i,j,k)*(1.0-x2)
  280  continue
c
c
c
c      three dimensional cartesian model
       time=0.0
       iiq=0
       iip=0
       iis=0
       vmax=0.0
c
       call clock(zt0)
       call clock(zt1)
c
c
       do 100 iia=1,last
c
       tilx=float(iia+last*(iii-1))/float(last)
c      cp(8)=til0+tild*tilx
       cp(8)=til0
c
c
       iiq=iiq+1
       iip=iip+1
       iis=iis+1
c
CC MPI BOUNDARY CONDITION START
c      boundary condition at nz=1 and nz=nz2
c      boundary condition at ny=1 and ny=ny2
       xx4=0.5*hx*float(2*nxp-nx1-2)
       xx3=hhx(nx1)+xx4
c
c
c      boundary condition in z-direction
c      do 31 k=1,nz2
       do 31 k=ks,ke
       if(k.eq.1 .and. irankz.eq.0) then
c      do 311 j=2,ny1
       do 311 j=js1,je1
c      do 311 i=2,nx1
       do 311 i=is1,ie1
       do 311 m=1,nb
c      f(m,i,j,k)=x*f(m,i-1,j,k+1)+xx*f(m,i-2,j,k+2)
       f(m,i,j,k)=f(m,i-1,j,k+1)
  311  continue
       else if(k.eq.ke .and. irankz.eq.isizez-1) then
c      do 312 j=2,ny1
       do 312 j=js1,je1
c      do 312 i=2,nx1
       do 312 i=is1,ie1
       do 312 m=1,nb
c      f(m,i,j,k)=x*f(m,i-1,j,k-1)+xx*f(m,i-2,j,k-2)
       f(m,i,j,k)=f(m,i-1,j,k-1)
  312  continue
       end if
   31  continue
c
c
c      boundary condition in y-direction
c      do 32 j=1,ny2
       do 32 j=js,je
       if(j.eq.1 .and. iranky.eq.0) then
c      do 321 k=1,nz2
       do 321 k=ks,ke
c      do 321 i=2,nx1
       do 321 i=is1,ie1
       do 321 m=1,nb
c      f(m,i,j,k)=x*f(m,i-1,j+1,k)+xx*f(m,i-2,j+2,k)
       f(m,i,j,k)=f(m,i-1,j+1,k)
  321  continue
       else if(j.eq.je .and. iranky.eq.isizey-1) then
c      do 322 k=1,nz2
       do 322 k=ks,ke
c      do 322 i=2,nx1
       do 322 i=is1,ie1
       do 322 m=1,nb
c      f(m,i,j,k)=x*f(m,i-1,j-1,k)+xx*f(m,i-2,j-2,k)
       f(m,i,j,k)=f(m,i-1,j-1,k)
  322  continue
       end if
   32  continue
c
c
c      boundary condition in x-direction
c      boundary condition at nx=nx2
c      do 33 i=1,nx2
       do 33 i=is,ie
       if(i.eq.ie .and. irankx.eq.isizex-1) then
c      do 331 k=1,nz2
       do 331 k=ks,ke
c      do 331 j=1,ny2
       do 331 j=js,je
       do 331 m=1,nb
       f(m,i,j,k)=f(m,i-1,j,k)
  331  continue
       end if
   33  continue
c
c
c      do 34 k=1,nz2
       do 34 k=ks,ke
c      do 34 j=1,ny2
       do 34 j=js,je
c      do 34 i=1,nx2
       do 34 i=is,ie
       f(5,i,j,k)=amax1(f(5,i,j,k),pr02)
       f(1,i,j,k)=amax1(f(1,i,j,k),ro02)
       u(5,i,j,k)=amax1(u(5,i,j,k),pr02)
       u(1,i,j,k)=amax1(u(1,i,j,k),ro02)
   34  continue
c
CC MPI BOUNDARY CONDITION END
c
c
!hpf$  reflect f
c
CC MPI START
       irightx = itable(irankx+1,iranky,irankz)
       ileftx  = itable(irankx-1,iranky,irankz)
       irighty = itable(irankx,iranky+1,irankz)
       ilefty  = itable(irankx,iranky-1,irankz)
       irightz = itable(irankx,iranky,irankz+1)
       ileftz  = itable(irankx,iranky,irankz-1)
c
       ftemp1x=f(:,is,:,:)
       ftemp1y=f(:,:,js,:)
       ftemp1z=f(:,:,:,ks)
       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(:,ie+1,:,:)=ftemp2x
       f(:,:,je+1,:)=ftemp2y
       f(:,:,:,ke+1)=ftemp2z
CC MPI END
c
       if(iia.ne.1.and.iiq.ne.iiq0) go to 172
c      do 372 m=1,nb
c      do 371 k=1,nz1
       do 371 k=ks,ke1
c      do 371 j=1,ny1
       do 371 j=js,je1
c      do 371 i=1,nx1
       do 371 i=is,ie1
       do 371 m=1,nb
       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))
  371  continue
c
c 372  continue
       t1=0.5*t
       dx1=0.25*t1/hx
       dy1=0.25*t1/hy
       dz1=0.25*t1/hz
       dx3=t1/(hx*hx)
       dy3=t1/(hy*hy)
       dz3=t1/(hz*hz)
c
c      if(nin.ne.0.or.iia.ne.1) go to 172
c      do 171 m=1,nb
c      do 171 k=1,nz2
c      write(ntap) ((gf(m,i,j,k),i=1,nx2),j=1,ny2)
c 171  continue
  172  continue
c
c
c
c      first step
c      step of k=k
c      do 90 k=1,nz1
       do 90 k=ks,ke1
c
c      current
c      do 38 j=1,ny1
       do 38 j=js,je1
c      do 38 i=1,nx1
       do 38 i=is,ie1
       p(11,i,j,k)=0.25*((f(7,i+1,j+1,k+1)+f(7,i+1,j,k+1)
     1   +f(7,i+1,j+1,k)+f(7,i+1,j,k)
     2   -f(7,i,j+1,k+1)-f(7,i,j,k+1)-f(7,i,j+1,k)-f(7,i,j,k))/hx
     3   -(f(6,i+1,j+1,k+1)-f(6,i+1,j,k+1)+f(6,i+1,j+1,k)-f(6,i+1,j,k)
     4   +f(6,i,j+1,k+1)-f(6,i,j,k+1)+f(6,i,j+1,k)-f(6,i,j,k))/hy)
       p(10,i,j,k)=0.25*((f(6,i+1,j+1,k+1)+f(6,i+1,j,k+1)
     1   -f(6,i+1,j+1,k)-f(6,i+1,j,k)
     2   +f(6,i,j+1,k+1)+f(6,i,j,k+1)-f(6,i,j+1,k)-f(6,i,j,k))/hz
     3   -(f(8,i+1,j+1,k+1)+f(8,i+1,j,k+1)+f(8,i+1,j+1,k)+f(8,i+1,j,k)
     4   -f(8,i,j+1,k+1)-f(8,i,j,k+1)-f(8,i,j+1,k)-f(8,i,j,k))/hx)
       p(9,i,j,k)=0.25*((f(8,i+1,j+1,k+1)-f(8,i+1,j,k+1)
     1   +f(8,i+1,j+1,k)-f(8,i+1,j,k)
     2   +f(8,i,j+1,k+1)-f(8,i,j,k+1)+f(8,i,j+1,k)-f(8,i,j,k))/hy
     3   -(f(7,i+1,j+1,k+1)+f(7,i+1,j,k+1)-f(7,i+1,j+1,k)-f(7,i+1,j,k)
     4   +f(7,i,j+1,k+1)+f(7,i,j,k+1)-f(7,i,j+1,k)-f(7,i,j,k))/hz)
   38  continue
c      do 382 j=1,ny1
       do 382 j=js,je1
c      do 382 i=1,nx1
       do 382 i=is,ie1
       p(11,i,j,k)=p(11,i,j,k)-pp(3,i,j,k)
       p(10,i,j,k)=p(10,i,j,k)-pp(2,i,j,k)
       p(9,i,j,k)=p(9,i,j,k)-pp(1,i,j,k)
  382  continue
c
c      do 381 j=1,ny1
       do 381 j=js,je1
c      do 381 i=1,nx1
       do 381 i=is,ie1
       do 381 m=1,nb
       p(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))
c      u(m,i,j,k)=p(m,i,j,k)
  381  continue
c
c      first step
       xx4=0.5*hx*float(2*nxp-nx1-1)
       z=0.5*hz*float(2*(kss+k)-nz2)
c      do 50 j=1,ny1
       do 50 j=js,je1
c      y=0.5*hy*float(2*(jss+j)-2)
       y=0.5*hy*float(2*(jss+j)-ny2)
c      do 502 i=1,nx1
       do 502 i=is,ie1
       do 502 m=1,nb
       ff(m,i,j,k)=u(m,i,j,k)
 502  continue
c
c      do 501 i=1,nx1
       do 501 i=is,ie1
       x=hhx(iss+i)+xx4
       ro2=x*x+y*y+z*z
       ro2=amax1(ro2,ar2)
       ro3=ro2*sqrt(ro2)
c
       u(8,i,j,k)=u(8,i,j,k)
     1   +dx1*(f(4,i+1,j+1,k+1)*f(6,i+1,j+1,k+1)
     1   -f(2,i+1,j+1,k+1)*f(8,i+1,j+1,k+1)
     2   +f(4,i+1,j,k+1)*f(6,i+1,j,k+1)-f(2,i+1,j,k+1)*f(8,i+1,j,k+1)
     3   +f(4,i+1,j+1,k)*f(6,i+1,j+1,k)-f(2,i+1,j+1,k)*f(8,i+1,j+1,k)
     4   +f(4,i+1,j,k)*f(6,i+1,j,k)-f(2,i+1,j,k)*f(8,i+1,j,k)
     5   -f(4,i,j+1,k+1)*f(6,i,j+1,k+1)+f(2,i,j+1,k+1)*f(8,i,j+1,k+1)
     6   -f(4,i,j,k+1)*f(6,i,j,k+1)+f(2,i,j,k+1)*f(8,i,j,k+1)
     7   -f(4,i,j+1,k)*f(6,i,j+1,k)+f(2,i,j+1,k)*f(8,i,j+1,k)
     8   -f(4,i,j,k)*f(6,i,j,k)+f(2,i,j,k)*f(8,i,j,k))
     1   -dy1*(f(3,i+1,j+1,k+1)*f(8,i+1,j+1,k+1)
     1   -f(4,i+1,j+1,k+1)*f(7,i+1,j+1,k+1)
     2   -f(3,i+1,j,k+1)*f(8,i+1,j,k+1)+f(4,i+1,j,k+1)*f(7,i+1,j,k+1)
     3   +f(3,i+1,j+1,k)*f(8,i+1,j+1,k)-f(4,i+1,j+1,k)*f(7,i+1,j+1,k)
     4   -f(3,i+1,j,k)*f(8,i+1,j,k)+f(4,i+1,j,k)*f(7,i+1,j,k)
     5   +f(3,i,j+1,k+1)*f(8,i,j+1,k+1)-f(4,i,j+1,k+1)*f(7,i,j+1,k+1)
     6   -f(3,i,j,k+1)*f(8,i,j,k+1)+f(4,i,j,k+1)*f(7,i,j,k+1)
     7   +f(3,i,j+1,k)*f(8,i,j+1,k)-f(4,i,j+1,k)*f(7,i,j+1,k)
     8   -f(3,i,j,k)*f(8,i,j,k)+f(4,i,j,k)*f(7,i,j,k))
c
       u(7,i,j,k)=u(7,i,j,k)
     1   +dz1*(f(3,i+1,j+1,k+1)*f(8,i+1,j+1,k+1)
     1   -f(4,i+1,j+1,k+1)*f(7,i+1,j+1,k+1)
     2   +f(3,i+1,j,k+1)*f(8,i+1,j,k+1)-f(4,i+1,j,k+1)*f(7,i+1,j,k+1)
     3   -f(3,i+1,j+1,k)*f(8,i+1,j+1,k)+f(4,i+1,j+1,k)*f(7,i+1,j+1,k)
     4   -f(3,i+1,j,k)*f(8,i+1,j,k)+f(4,i+1,j,k)*f(7,i+1,j,k)
     5   +f(3,i,j+1,k+1)*f(8,i,j+1,k+1)-f(4,i,j+1,k+1)*f(7,i,j+1,k+1)
     6   +f(3,i,j,k+1)*f(8,i,j,k+1)-f(4,i,j,k+1)*f(7,i,j,k+1)
     7   -f(3,i,j+1,k)*f(8,i,j+1,k)+f(4,i,j+1,k)*f(7,i,j+1,k)
     8   -f(3,i,j,k)*f(8,i,j,k)+f(4,i,j,k)*f(7,i,j,k))
     1   -dx1*(f(2,i+1,j+1,k+1)*f(7,i+1,j+1,k+1)
     1   -f(3,i+1,j+1,k+1)*f(6,i+1,j+1,k+1)
     2   +f(2,i+1,j,k+1)*f(7,i+1,j,k+1)-f(3,i+1,j,k+1)*f(6,i+1,j,k+1)
     3   +f(2,i+1,j+1,k)*f(7,i+1,j+1,k)-f(3,i+1,j+1,k)*f(6,i+1,j+1,k)
     4   +f(2,i+1,j,k)*f(7,i+1,j,k)-f(3,i+1,j,k)*f(6,i+1,j,k)
     5   -f(2,i,j+1,k+1)*f(7,i,j+1,k+1)+f(3,i,j+1,k+1)*f(6,i,j+1,k+1)
     6   -f(2,i,j,k+1)*f(7,i,j,k+1)+f(3,i,j,k+1)*f(6,i,j,k+1)
     7   -f(2,i,j+1,k)*f(7,i,j+1,k)+f(3,i,j+1,k)*f(6,i,j+1,k)
     8   -f(2,i,j,k)*f(7,i,j,k)+f(3,i,j,k)*f(6,i,j,k))
c
       u(6,i,j,k)=u(6,i,j,k)
     1   +dy1*(f(2,i+1,j+1,k+1)*f(7,i+1,j+1,k+1)
     1   -f(3,i+1,j+1,k+1)*f(6,i+1,j+1,k+1)
     2   -f(2,i+1,j,k+1)*f(7,i+1,j,k+1)+f(3,i+1,j,k+1)*f(6,i+1,j,k+1)
     3   +f(2,i+1,j+1,k)*f(7,i+1,j+1,k)-f(3,i+1,j+1,k)*f(6,i+1,j+1,k)
     4   -f(2,i+1,j,k)*f(7,i+1,j,k)+f(3,i+1,j,k)*f(6,i+1,j,k)
     5   +f(2,i,j+1,k+1)*f(7,i,j+1,k+1)-f(3,i,j+1,k+1)*f(6,i,j+1,k+1)
     6   -f(2,i,j,k+1)*f(7,i,j,k+1)+f(3,i,j,k+1)*f(6,i,j,k+1)
     7   +f(2,i,j+1,k)*f(7,i,j+1,k)-f(3,i,j+1,k)*f(6,i,j+1,k)
     8   -f(2,i,j,k)*f(7,i,j,k)+f(3,i,j,k)*f(6,i,j,k))
     1   -dz1*(f(4,i+1,j+1,k+1)*f(6,i+1,j+1,k+1)
     1   -f(2,i+1,j+1,k+1)*f(8,i+1,j+1,k+1)
     2   +f(4,i+1,j,k+1)*f(6,i+1,j,k+1)-f(2,i+1,j,k+1)*f(8,i+1,j,k+1)
     3   -f(4,i+1,j+1,k)*f(6,i+1,j+1,k)+f(2,i+1,j+1,k)*f(8,i+1,j+1,k)
     4   -f(4,i+1,j,k)*f(6,i+1,j,k)+f(2,i+1,j,k)*f(8,i+1,j,k)
     5   +f(4,i,j+1,k+1)*f(6,i,j+1,k+1)-f(2,i,j+1,k+1)*f(8,i,j+1,k+1)
     6   +f(4,i,j,k+1)*f(6,i,j,k+1)-f(2,i,j,k+1)*f(8,i,j,k+1)
     7   -f(4,i,j+1,k)*f(6,i,j+1,k)+f(2,i,j+1,k)*f(8,i,j+1,k)
     8   -f(4,i,j,k)*f(6,i,j,k)+f(2,i,j,k)*f(8,i,j,k))
c
       x3=dx1*(f(2,i+1,j+1,k+1)+f(2,i+1,j,k+1)
     1   +f(2,i+1,j+1,k)+f(2,i+1,j,k)
     1   -f(2,i,j+1,k+1)-f(2,i,j,k+1)-f(2,i,j+1,k)-f(2,i,j,k))
     2   +dy1*(f(3,i+1,j+1,k+1)-f(3,i+1,j,k+1)
     2   +f(3,i+1,j+1,k)-f(3,i+1,j,k)
     3   +f(3,i,j+1,k+1)-f(3,i,j,k+1)+f(3,i,j+1,k)-f(3,i,j,k))
     4   +dz1*(f(4,i+1,j+1,k+1)+f(4,i+1,j,k+1)
     4   -f(4,i+1,j+1,k)-f(4,i+1,j,k)
     5   +f(4,i,j+1,k+1)+f(4,i,j,k+1)-f(4,i,j+1,k)-f(4,i,j,k))
c
       u(5,i,j,k)=u(5,i,j,k)-p(5,i,j,k)*gam*x3
     1   -dx1*p(2,i,j,k)*(f(5,i+1,j+1,k+1)+f(5,i+1,j,k+1)
     1   +f(5,i+1,j+1,k)+f(5,i+1,j,k)
     2   -f(5,i,j+1,k+1)-f(5,i,j,k+1)-f(5,i,j+1,k)-f(5,i,j,k))
     3   -dy1*p(3,i,j,k)*(f(5,i+1,j+1,k+1)-f(5,i+1,j,k+1)
     3   +f(5,i+1,j+1,k)-f(5,i+1,j,k)
     4   +f(5,i,j+1,k+1)-f(5,i,j,k+1)+f(5,i,j+1,k)-f(5,i,j,k))
     5   -dz1*p(4,i,j,k)*(f(5,i+1,j+1,k+1)+f(5,i+1,j,k+1)
     5   -f(5,i+1,j+1,k)-f(5,i+1,j,k)
     6   +f(5,i,j+1,k+1)+f(5,i,j,k+1)-f(5,i,j+1,k)-f(5,i,j,k))
c
c
       u(4,i,j,k)=u(4,i,j,k)+t1*((p(9,i,j,k)*p(7,i,j,k)
     1   -p(10,i,j,k)*p(6,i,j,k))/p(1,i,j,k)
     1   -0.0*eud*p(4,i,j,k)-gra*z/ro3)
     2   -dx1*p(2,i,j,k)*(f(4,i+1,j+1,k+1)+f(4,i+1,j,k+1)
     2   +f(4,i+1,j+1,k)+f(4,i+1,j,k)
     3   -f(4,i,j+1,k+1)-f(4,i,j,k+1)-f(4,i,j+1,k)-f(4,i,j,k))
     4   -dy1*p(3,i,j,k)*(f(4,i+1,j+1,k+1)-f(4,i+1,j,k+1)
     4   +f(4,i+1,j+1,k)-f(4,i+1,j,k)
     5   +f(4,i,j+1,k+1)-f(4,i,j,k+1)+f(4,i,j+1,k)-f(4,i,j,k))
     6   -dz1*p(4,i,j,k)*(f(4,i+1,j+1,k+1)+f(4,i+1,j,k+1)
     6   -f(4,i+1,j+1,k)-f(4,i+1,j,k)
     7   +f(4,i,j+1,k+1)+f(4,i,j,k+1)-f(4,i,j+1,k)-f(4,i,j,k))
     8   -dz1*(f(5,i+1,j+1,k+1)+f(5,i+1,j,k+1)
     8   -f(5,i+1,j+1,k)-f(5,i+1,j,k)
     9   +f(5,i,j+1,k+1)+f(5,i,j,k+1)-f(5,i,j+1,k)-f(5,i,j,k))
     9   /p(1,i,j,k)
c
       u(3,i,j,k)=u(3,i,j,k)+t1*((p(11,i,j,k)*p(6,i,j,k)
     1   -p(9,i,j,k)*p(8,i,j,k))/p(1,i,j,k)
     1   -0.0*eud*p(3,i,j,k)-gra*y/ro3)
     2   -dx1*p(2,i,j,k)*(f(3,i+1,j+1,k+1)+f(3,i+1,j,k+1)
     2   +f(3,i+1,j+1,k)+f(3,i+1,j,k)
     3   -f(3,i,j+1,k+1)-f(3,i,j,k+1)-f(3,i,j+1,k)-f(3,i,j,k))
     4   -dy1*p(3,i,j,k)*(f(3,i+1,j+1,k+1)-f(3,i+1,j,k+1)
     4   +f(3,i+1,j+1,k)-f(3,i+1,j,k)
     5   +f(3,i,j+1,k+1)-f(3,i,j,k+1)+f(3,i,j+1,k)-f(3,i,j,k))
     6   -dz1*p(4,i,j,k)*(f(3,i+1,j+1,k+1)+f(3,i+1,j,k+1)
     6   -f(3,i+1,j+1,k)-f(3,i+1,j,k)
     7   +f(3,i,j+1,k+1)+f(3,i,j,k+1)-f(3,i,j+1,k)-f(3,i,j,k))
     8   -dy1*(f(5,i+1,j+1,k+1)-f(5,i+1,j,k+1)
     8   +f(5,i+1,j+1,k)-f(5,i+1,j,k)
     9   +f(5,i,j+1,k+1)-f(5,i,j,k+1)+f(5,i,j+1,k)-f(5,i,j,k))
     9   /p(1,i,j,k)
c
c      uxd=amax1((x-4.0*ar1)/(4.0*ar1),0.0)
       uxd=amax1((x-xta)/xtd,0.0)
       uxd=amin1(uxd,1.0)
       uxd=amin1(uxd*p(2,i,j,k),0.0)
c
       u(2,i,j,k)=u(2,i,j,k)+t1*((p(10,i,j,k)*p(8,i,j,k)
     1   -p(11,i,j,k)*p(7,i,j,k))/p(1,i,j,k)
     1   -eud*uxd-gra*x/ro3)
     2   -dx1*p(2,i,j,k)*(f(2,i+1,j+1,k+1)+f(2,i+1,j,k+1)
     2   +f(2,i+1,j+1,k)+f(2,i+1,j,k)
     3   -f(2,i,j+1,k+1)-f(2,i,j,k+1)-f(2,i,j+1,k)-f(2,i,j,k))
     4   -dy1*p(3,i,j,k)*(f(2,i+1,j+1,k+1)-f(2,i+1,j,k+1)
     4   +f(2,i+1,j+1,k)-f(2,i+1,j,k)
     5   +f(2,i,j+1,k+1)-f(2,i,j,k+1)+f(2,i,j+1,k)-f(2,i,j,k))
     6   -dz1*p(4,i,j,k)*(f(2,i+1,j+1,k+1)+f(2,i+1,j,k+1)
     6   -f(2,i+1,j+1,k)-f(2,i+1,j,k)
     7   +f(2,i,j+1,k+1)+f(2,i,j,k+1)-f(2,i,j+1,k)-f(2,i,j,k))
     8   -dx1*(f(5,i+1,j+1,k+1)+f(5,i+1,j,k+1)
     8   +f(5,i+1,j+1,k)+f(5,i+1,j,k)
     9   -f(5,i,j+1,k+1)-f(5,i,j,k+1)-f(5,i,j+1,k)-f(5,i,j,k))
     9   /p(1,i,j,k)
c
       u(1,i,j,k)=u(1,i,j,k)
     1   -dx1*(f(2,i+1,j+1,k+1)*f(1,i+1,j+1,k+1)
     1   +f(2,i+1,j,k+1)*f(1,i+1,j,k+1)
     2   +f(2,i+1,j+1,k)*f(1,i+1,j+1,k)+f(2,i+1,j,k)*f(1,i+1,j,k)
     3   -f(2,i,j+1,k+1)*f(1,i,j+1,k+1)-f(2,i,j,k+1)*f(1,i,j,k+1)
     4   -f(2,i,j+1,k)*f(1,i,j+1,k)-f(2,i,j,k)*f(1,i,j,k))
     5   -dy1*(f(3,i+1,j+1,k+1)*f(1,i+1,j+1,k+1)
     5   -f(3,i+1,j,k+1)*f(1,i+1,j,k+1)
     6   +f(3,i+1,j+1,k)*f(1,i+1,j+1,k)-f(3,i+1,j,k)*f(1,i+1,j,k)
     7   +f(3,i,j+1,k+1)*f(1,i,j+1,k+1)-f(3,i,j,k+1)*f(1,i,j,k+1)
     8   +f(3,i,j+1,k)*f(1,i,j+1,k)-f(3,i,j,k)*f(1,i,j,k))
     1   -dz1*(f(4,i+1,j+1,k+1)*f(1,i+1,j+1,k+1)
     1   +f(4,i+1,j,k+1)*f(1,i+1,j,k+1)
     2   -f(4,i+1,j+1,k)*f(1,i+1,j+1,k)-f(4,i+1,j,k)*f(1,i+1,j,k)
     3   +f(4,i,j+1,k+1)*f(1,i,j+1,k+1)+f(4,i,j,k+1)*f(1,i,j,k+1)
     4   -f(4,i,j+1,k)*f(1,i,j+1,k)-f(4,i,j,k)*f(1,i,j,k))
  501  continue
c
c      do 520 i=1,nx1
       do 520 i=is,ie1
       x=hhx(iss+i)+xx4
       ro2=x*x+y*y+z*z
       x2=amax1(ro2/ar2-1.0,0.0)
       x2=aru*x2*x2
       x2=x2/(x2+1.0)
       u(8,i,j,k)=u(8,i,j,k)*x2+ff(8,i,j,k)*(1.0-x2)
       u(7,i,j,k)=u(7,i,j,k)*x2+ff(7,i,j,k)*(1.0-x2)
       u(6,i,j,k)=u(6,i,j,k)*x2+ff(6,i,j,k)*(1.0-x2)
       u(5,i,j,k)=u(5,i,j,k)*x2+ff(5,i,j,k)*(1.0-x2)
       u(4,i,j,k)=u(4,i,j,k)*x2+ff(4,i,j,k)*(1.0-x2)
       u(3,i,j,k)=u(3,i,j,k)*x2+ff(3,i,j,k)*(1.0-x2)
       u(2,i,j,k)=u(2,i,j,k)*x2+ff(2,i,j,k)*(1.0-x2)
       u(1,i,j,k)=u(1,i,j,k)*x2+ff(1,i,j,k)*(1.0-x2)
  520  continue
c
c      do 524 i=1,nx1
       do 524 i=is,ie1
       x5=u(5,i,j,k)
       x1=u(1,i,j,k)
       y5=amax1(x5,pr02)
       y1=amax1(x1,ro02)
c      x2=u(2,i,j,k)*u(2,i,j,k)+u(3,i,j,k)*u(3,i,j,k)+u(4,i,j,k)*u(4,i,j,k)
c      vmax=amax1(vmax,x2)
       x=amax1(x1,0.0)/y1
       u(5,i,j,k)=y5
       u(4,i,j,k)=x*u(4,i,j,k)
       u(3,i,j,k)=x*u(3,i,j,k)
       u(2,i,j,k)=x*u(2,i,j,k)
       u(1,i,j,k)=y1
  524  continue
   50  continue
   90  continue
c
c
c
c      do 70 m=1,nb
c      do 72 k=1,nz2
       do 72 k=ks,ke
c      do 72 j=1,ny2
       do 72 j=js,je
c      do 72 i=1,nx2
       do 72 i=is,ie
       do 72 m=1,nb
       ff(m,i,j,k)=f(m,i,j,k)
   72  continue
c  70  continue
c
c
!hpf$  reflect u
!hpf$  reflect ff
!hpf$  reflect pp
c
CC MPI START
c
       irightx = itable(irankx+1,iranky,irankz)
       ileftx  = itable(irankx-1,iranky,irankz)
       irighty = itable(irankx,iranky+1,irankz)
       ilefty  = itable(irankx,iranky-1,irankz)
       irightz = itable(irankx,iranky,irankz+1)
       ileftz  = itable(irankx,iranky,irankz-1)
c
       ftemp1x=u(:,ie,:,:)
       ftemp1y=u(:,:,je,:)
       ftemp1z=u(:,:,:,ke)
       call mpi_sendrecv(ftemp1x,nwyz,mpi_real,irightx,310,
     &                   ftemp2x,nwyz,mpi_real,ileftx,310,
     &                   mpi_comm_world,istatus,ier)
       call mpi_sendrecv(ftemp1y,nwzx,mpi_real,irighty,320,
     &                   ftemp2y,nwzx,mpi_real,ilefty,320,
     &                   mpi_comm_world,istatus,ier)
       call mpi_sendrecv(ftemp1z,nwxy,mpi_real,irightz,330,
     &                   ftemp2z,nwxy,mpi_real,ileftz,330,
     &                   mpi_comm_world,istatus,ier)
       u(:,is-1,:,:)=ftemp2x
       u(:,:,js-1,:)=ftemp2y
       u(:,:,:,ks-1)=ftemp2z
c
       ftemp1x=ff(:,ie,:,:)
       ftemp1y=ff(:,:,je,:)
       ftemp1z=ff(:,:,:,ke)
       call mpi_sendrecv(ftemp1x,nwyz,mpi_real,irightx,410,
     &                   ftemp2x,nwyz,mpi_real,ileftx,410,
     &                   mpi_comm_world,istatus,ier)
       call mpi_sendrecv(ftemp1y,nwzx,mpi_real,irighty,420,
     &                   ftemp2y,nwzx,mpi_real,ilefty,420,
     &                   mpi_comm_world,istatus,ier)
       call mpi_sendrecv(ftemp1z,nwxy,mpi_real,irightz,430,
     &                   ftemp2z,nwxy,mpi_real,ileftz,430,
     &                   mpi_comm_world,istatus,ier)
       ff(:,is-1,:,:)=ftemp2x
       ff(:,:,js-1,:)=ftemp2y
       ff(:,:,:,ks-1)=ftemp2z
c
       ftemp1x=ff(:,is,:,:)
       ftemp1y=ff(:,:,js,:)
       ftemp1z=ff(:,:,:,ks)
       call mpi_sendrecv(ftemp1x,nwyz,mpi_real,ileftx,470,
     &                   ftemp2x,nwyz,mpi_real,irightx,470,
     &                   mpi_comm_world,istatus,ier)
       call mpi_sendrecv(ftemp1y,nwzx,mpi_real,ilefty,480,
     &                   ftemp2y,nwzx,mpi_real,irighty,480,
     &                   mpi_comm_world,istatus,ier)
       call mpi_sendrecv(ftemp1z,nwxy,mpi_real,ileftz,490,
     &                   ftemp2z,nwxy,mpi_real,irightz,490,
     &                   mpi_comm_world,istatus,ier)
       ff(:,ie+1,:,:)=ftemp2x
       ff(:,:,je+1,:)=ftemp2y
       ff(:,:,:,ke+1)=ftemp2z
c
       utemp1x=pp(:,ie,:,:)
       utemp1y=pp(:,:,je,:)
       utemp1z=pp(:,:,:,ke)
       call mpi_sendrecv(utemp1x,mwyz,mpi_real,irightx,510,
     &                   utemp2x,mwyz,mpi_real,ileftx,510,
     &                   mpi_comm_world,istatus,ier)
       call mpi_sendrecv(utemp1y,mwzx,mpi_real,irighty,520,
     &                   utemp2y,mwzx,mpi_real,ilefty,520,
     &                   mpi_comm_world,istatus,ier)
       call mpi_sendrecv(utemp1z,mwxy,mpi_real,irightz,530,
     &                   utemp2z,mwxy,mpi_real,ileftz,530,
     &                   mpi_comm_world,istatus,ier)
       pp(:,is-1,:,:)=utemp2x
       pp(:,:,js-1,:)=utemp2y
       pp(:,:,:,ks-1)=utemp2z
c
CC MPI END
c
c
c      second step
c
c
c      do 190 k=2,nz1
       do 190 k=ks1,ke1
c
c
c      current
c      do 68 j=2,ny1
       do 68 j=js1,je1
c      do 68 i=2,nx1
       do 68 i=is1,ie1
       p(11,i,j,k)=0.25*((u(7,i,j,k)+u(7,i,j-1,k)
     1        +u(7,i,j,k-1)+u(7,i,j-1,k-1)
     2        -u(7,i-1,j,k)-u(7,i-1,j-1,k)
     3        -u(7,i-1,j,k-1)-u(7,i-1,j-1,k-1))/hx
     4        -(u(6,i,j,k)-u(6,i,j-1,k)
     5        +u(6,i,j,k-1)-u(6,i,j-1,k-1)
     6        +u(6,i-1,j,k)-u(6,i-1,j-1,k)
     7        +u(6,i-1,j,k-1)-u(6,i-1,j-1,k-1))/hy)
       p(10,i,j,k)=0.25*((u(6,i,j,k)+u(6,i,j-1,k)
     1        -u(6,i,j,k-1)-u(6,i,j-1,k-1)
     2        +u(6,i-1,j,k)+u(6,i-1,j-1,k)
     3        -u(6,i-1,j,k-1)-u(6,i-1,j-1,k-1))/hz
     4        -(u(8,i,j,k)+u(8,i,j-1,k)
     5        +u(8,i,j,k-1)+u(8,i,j-1,k-1)
     6        -u(8,i-1,j,k)-u(8,i-1,j-1,k)
     7        -u(8,i-1,j,k-1)-u(8,i-1,j-1,k-1))/hx)
       p(9,i,j,k)=0.25*((u(8,i,j,k)-u(8,i,j-1,k)
     1       +u(8,i,j,k-1)-u(8,i,j-1,k-1)
     2       +u(8,i-1,j,k)-u(8,i-1,j-1,k)
     3       +u(8,i-1,j,k-1)-u(8,i-1,j-1,k-1))/hy
     4       -(u(7,i,j,k)+u(7,i,j-1,k)
     5       -u(7,i,j,k-1)-u(7,i,j-1,k-1)
     6       +u(7,i-1,j,k)+u(7,i-1,j-1,k)
     7       -u(7,i-1,j,k-1)-u(7,i-1,j-1,k-1))/hz)
   68  continue
c      do 682 j=2,ny1
       do 682 j=js1,je1
c      do 682 i=2,nx1
       do 682 i=is1,ie1
       p(11,i,j,k)=p(10,i,j,k)-0.125*(
     1          pp(3,i,j,k)+pp(3,i-1,j,k)
     2         +pp(3,i,j-1,k)+pp(3,i-1,j-1,k)
     3         +pp(3,i,j,k-1)+pp(3,i-1,j,k-1)
     4         +pp(3,i,j-1,k-1)+pp(3,i-1,j-1,k-1))
       p(10,i,j,k)=p(10,i,j,k)-0.125*(
     1          pp(2,i,j,k)+pp(2,i-1,j,k)
     2         +pp(2,i,j-1,k)+pp(2,i-1,j-1,k)
     3         +pp(2,i,j,k-1)+pp(2,i-1,j,k-1)
     4         +pp(2,i,j-1,k-1)+pp(2,i-1,j-1,k-1))
       p(9,i,j,k)=p(9,i,j,k)-0.125*(
     1          pp(1,i,j,k)+pp(1,i-1,j,k)
     2         +pp(1,i,j-1,k)+pp(1,i-1,j-1,k)
     3         +pp(1,i,j,k-1)+pp(1,i-1,j,k-1)
     4         +pp(1,i,j-1,k-1)+pp(1,i-1,j-1,k-1))
  682  continue
c
c      do 681 j=2,ny1
       do 681 j=js1,je1
c      do 681 i=2,nx1
       do 681 i=is1,ie1
       do 681 m=1,nb
       p(m,i,j,k)=0.125*(u(m,i,j,k)+u(m,i-1,j,k)
     1      +u(m,i,j-1,k)+u(m,i-1,j-1,k)
     2      +u(m,i,j,k-1)+u(m,i-1,j,k-1)
     3      +u(m,i,j-1,k-1)+u(m,i-1,j-1,k-1))
  681  continue
c
c      second step
       xx4=0.5*hx*float(2*nxp-nx1-2)
       z=0.5*hz*float(2*(kss+k)-nz2-1)
c      do 80 j=2,ny1
       do 80 j=js1,je1
c      y=0.5*hy*float(2*(jss+j)-3)
       y=0.5*hy*float(2*(jss+j)-ny2-1)
c
c      do 801 i=2,nx1
       do 801 i=is1,ie1
       x=hhx(iss+i)+xx4
       ro2=x*x+y*y+z*z
       ro2=amax1(ro2,ar2)
       ro3=ro2*sqrt(ro2)
c
       x1=ff(5,i,j,k)/(ff(1,i,j,k)*po0)
       x2=ro2/ar2-1.0
       x2=aru*x2*x2
       x2=x2/(x2+1.0)
       eat=eat0*x2
c      eat=eat0*x2/sqrt(x1*x1*x1)
c
       x1=eat*(dx4*(ff(8,i+1,j,k)-2.0*ff(8,i,j,k)+ff(8,i-1,j,k))
     1        +dy4*(ff(8,i,j+1,k)-2.0*ff(8,i,j,k)+ff(8,i,j-1,k))
     2        +dz4*(ff(8,i,j,k+1)-2.0*ff(8,i,j,k)+ff(8,i,j,k-1)))
       f(8,i,j,k)=ff(8,i,j,k)+x1
     1   +dx2*(u(4,i,j,k)*u(6,i,j,k)-u(2,i,j,k)*u(8,i,j,k)
     2   +u(4,i,j-1,k)*u(6,i,j-1,k)-u(2,i,j-1,k)*u(8,i,j-1,k)
     3   +u(4,i,j,k-1)*u(6,i,j,k-1)-u(2,i,j,k-1)*u(8,i,j,k-1)
     4   +u(4,i,j-1,k-1)*u(6,i,j-1,k-1)-u(2,i,j-1,k-1)*u(8,i,j-1,k-1)
     5   -u(4,i-1,j,k)*u(6,i-1,j,k)+u(2,i-1,j,k)*u(8,i-1,j,k)
     6   -u(4,i-1,j-1,k)*u(6,i-1,j-1,k)+u(2,i-1,j-1,k)*u(8,i-1,j-1,k)
     7   -u(4,i-1,j,k-1)*u(6,i-1,j,k-1)+u(2,i-1,j,k-1)*u(8,i-1,j,k-1)
     8   -u(4,i-1,j-1,k-1)*u(6,i-1,j-1,k-1)
     8   +u(2,i-1,j-1,k-1)*u(8,i-1,j-1,k-1))
     1   -dy2*(u(3,i,j,k)*u(8,i,j,k)-u(4,i,j,k)*u(7,i,j,k)
     2   -u(3,i,j-1,k)*u(8,i,j-1,k)+u(4,i,j-1,k)*u(7,i,j-1,k)
     3   +u(3,i,j,k-1)*u(8,i,j,k-1)-u(4,i,j,k-1)*u(7,i,j,k-1)
     4   -u(3,i,j-1,k-1)*u(8,i,j-1,k-1)+u(4,i,j-1,k-1)*u(7,i,j-1,k-1)
     5   +u(3,i-1,j,k)*u(8,i-1,j,k)-u(4,i-1,j,k)*u(7,i-1,j,k)
     6   -u(3,i-1,j-1,k)*u(8,i-1,j-1,k)+u(4,i-1,j-1,k)*u(7,i-1,j-1,k)
     7   +u(3,i-1,j,k-1)*u(8,i-1,j,k-1)-u(4,i-1,j,k-1)*u(7,i-1,j,k-1)
     8   -u(3,i-1,j-1,k-1)*u(8,i-1,j-1,k-1)
     8   +u(4,i-1,j-1,k-1)*u(7,i-1,j-1,k-1))
c
       x1=eat*(dx4*(ff(7,i+1,j,k)-2.0*ff(7,i,j,k)+ff(7,i-1,j,k))
     1        +dy4*(ff(7,i,j+1,k)-2.0*ff(7,i,j,k)+ff(7,i,j-1,k))
     2        +dz4*(ff(7,i,j,k+1)-2.0*ff(7,i,j,k)+ff(7,i,j,k-1)))
       f(7,i,j,k)=ff(7,i,j,k)+x1
     1   +dz2*(u(3,i,j,k)*u(8,i,j,k)-u(4,i,j,k)*u(7,i,j,k)
     2   +u(3,i,j-1,k)*u(8,i,j-1,k)-u(4,i,j-1,k)*u(7,i,j-1,k)
     3   -u(3,i,j,k-1)*u(8,i,j,k-1)+u(4,i,j,k-1)*u(7,i,j,k-1)
     4   -u(3,i,j-1,k-1)*u(8,i,j-1,k-1)+u(4,i,j-1,k-1)*u(7,i,j-1,k-1)
     5   +u(3,i-1,j,k)*u(8,i-1,j,k)-u(4,i-1,j,k)*u(7,i-1,j,k)
     6   +u(3,i-1,j-1,k)*u(8,i-1,j-1,k)-u(4,i-1,j-1,k)*u(7,i-1,j-1,k)
     7   -u(3,i-1,j,k-1)*u(8,i-1,j,k-1)+u(4,i-1,j,k-1)*u(7,i-1,j,k-1)
     8   -u(3,i-1,j-1,k-1)*u(8,i-1,j-1,k-1)
     8   +u(4,i-1,j-1,k-1)*u(7,i-1,j-1,k-1))
     1   -dx2*(u(2,i,j,k)*u(7,i,j,k)-u(3,i,j,k)*u(6,i,j,k)
     2   +u(2,i,j-1,k)*u(7,i,j-1,k)-u(3,i,j-1,k)*u(6,i,j-1,k)
     3   +u(2,i,j,k-1)*u(7,i,j,k-1)-u(3,i,j,k-1)*u(6,i,j,k-1)
     4   +u(2,i,j-1,k-1)*u(7,i,j-1,k-1)-u(3,i,j-1,k-1)*u(6,i,j-1,k-1)
     5   -u(2,i-1,j,k)*u(7,i-1,j,k)+u(3,i-1,j,k)*u(6,i-1,j,k)
     6   -u(2,i-1,j-1,k)*u(7,i-1,j-1,k)+u(3,i-1,j-1,k)*u(6,i-1,j-1,k)
     7   -u(2,i-1,j,k-1)*u(7,i-1,j,k-1)+u(3,i-1,j,k-1)*u(6,i-1,j,k-1)
     8   -u(2,i-1,j-1,k-1)*u(7,i-1,j-1,k-1)
     8   +u(3,i-1,j-1,k-1)*u(6,i-1,j-1,k-1))
c
       x1=eat*(dx4*(ff(6,i+1,j,k)-2.0*ff(6,i,j,k)+ff(6,i-1,j,k))
     1        +dy4*(ff(6,i,j+1,k)-2.0*ff(6,i,j,k)+ff(6,i,j-1,k))
     2        +dz4*(ff(6,i,j,k+1)-2.0*ff(6,i,j,k)+ff(6,i,j,k-1)))
       f(6,i,j,k)=ff(6,i,j,k)+x1
     1   +dy2*(u(2,i,j,k)*u(7,i,j,k)-u(3,i,j,k)*u(6,i,j,k)
     2   -u(2,i,j-1,k)*u(7,i,j-1,k)+u(3,i,j-1,k)*u(6,i,j-1,k)
     3   +u(2,i,j,k-1)*u(7,i,j,k-1)-u(3,i,j,k-1)*u(6,i,j,k-1)
     4   -u(2,i,j-1,k-1)*u(7,i,j-1,k-1)+u(3,i,j-1,k-1)*u(6,i,j-1,k-1)
     5   +u(2,i-1,j,k)*u(7,i-1,j,k)-u(3,i-1,j,k)*u(6,i-1,j,k)
     6   -u(2,i-1,j-1,k)*u(7,i-1,j-1,k)+u(3,i-1,j-1,k)*u(6,i-1,j-1,k)
     7   +u(2,i-1,j,k-1)*u(7,i-1,j,k-1)-u(3,i-1,j,k-1)*u(6,i-1,j,k-1)
     8   -u(2,i-1,j-1,k-1)*u(7,i-1,j-1,k-1)
     8   +u(3,i-1,j-1,k-1)*u(6,i-1,j-1,k-1))
     1   -dz2*(u(4,i,j,k)*u(6,i,j,k)-u(2,i,j,k)*u(8,i,j,k)
     2   +u(4,i,j-1,k)*u(6,i,j-1,k)-u(2,i,j-1,k)*u(8,i,j-1,k)
     3   -u(4,i,j,k-1)*u(6,i,j,k-1)+u(2,i,j,k-1)*u(8,i,j,k-1)
     4   -u(4,i,j-1,k-1)*u(6,i,j-1,k-1)+u(2,i,j-1,k-1)*u(8,i,j-1,k-1)
     5   +u(4,i-1,j,k)*u(6,i-1,j,k)-u(2,i-1,j,k)*u(8,i-1,j,k)
     6   +u(4,i-1,j-1,k)*u(6,i-1,j-1,k)-u(2,i-1,j-1,k)*u(8,i-1,j-1,k)
     7   -u(4,i-1,j,k-1)*u(6,i-1,j,k-1)+u(2,i-1,j,k-1)*u(8,i-1,j,k-1)
     8   -u(4,i-1,j-1,k-1)*u(6,i-1,j-1,k-1)
     8   +u(2,i-1,j-1,k-1)*u(8,i-1,j-1,k-1))
c
       x3=dx2*(u(2,i,j,k)+u(2,i,j-1,k)+u(2,i,j,k-1)+u(2,i,j-1,k-1)
     1   -u(2,i-1,j,k)-u(2,i-1,j-1,k)-u(2,i-1,j,k-1)-u(2,i-1,j-1,k-1))
     2   +dy2*(u(3,i,j,k)-u(3,i,j-1,k)+u(3,i,j,k-1)-u(3,i,j-1,k-1)
     3   +u(3,i-1,j,k)-u(3,i-1,j-1,k)+u(3,i-1,j,k-1)-u(3,i-1,j-1,k-1))
     4   +dz2*(u(4,i,j,k)+u(4,i,j-1,k)-u(4,i,j,k-1)-u(4,i,j-1,k-1)
     5   +u(4,i-1,j,k)+u(4,i-1,j-1,k)-u(4,i-1,j,k-1)-u(4,i-1,j-1,k-1))
       f(5,i,j,k)=ff(5,i,j,k)-p(5,i,j,k)*gam*x3
     1   -dx2*p(2,i,j,k)*
     1   (u(5,i,j,k)+u(5,i,j-1,k)+u(5,i,j,k-1)+u(5,i,j-1,k-1)
     2   -u(5,i-1,j,k)-u(5,i-1,j-1,k)-u(5,i-1,j,k-1)-u(5,i-1,j-1,k-1))
     3   -dy2*p(3,i,j,k)*
     3   (u(5,i,j,k)-u(5,i,j-1,k)+u(5,i,j,k-1)-u(5,i,j-1,k-1)
     4   +u(5,i-1,j,k)-u(5,i-1,j-1,k)+u(5,i-1,j,k-1)-u(5,i-1,j-1,k-1))
     5   -dz2*p(4,i,j,k)*
     5   (u(5,i,j,k)+u(5,i,j-1,k)-u(5,i,j,k-1)-u(5,i,j-1,k-1)
     6   +u(5,i-1,j,k)+u(5,i-1,j-1,k)-u(5,i-1,j,k-1)-u(5,i-1,j-1,k-1))
     7   +pmu*(dx4*(ff(5,i+1,j,k)-2.0*ff(5,i,j,k)+ff(5,i-1,j,k))
     8   +dy4*(ff(5,i,j+1,k)-2.0*ff(5,i,j,k)+ff(5,i,j-1,k))
     9   +dz4*(ff(5,i,j,k+1)-2.0*ff(5,i,j,k)+ff(5,i,j,k-1)))
c
       f(4,i,j,k)=ff(4,i,j,k)+t*((p(9,i,j,k)*p(7,i,j,k)
     1   -p(10,i,j,k)*p(6,i,j,k))/p(1,i,j,k)
     1   -0.0*eud*p(4,i,j,k)-gra*z/ro3)
     2   -dx2*p(2,i,j,k)*
     2   (u(4,i,j,k)+u(4,i,j-1,k)+u(4,i,j,k-1)+u(4,i,j-1,k-1)
     3   -u(4,i-1,j,k)-u(4,i-1,j-1,k)-u(4,i-1,j,k-1)-u(4,i-1,j-1,k-1))
     4   -dy2*p(3,i,j,k)*
     4   (u(4,i,j,k)-u(4,i,j-1,k)+u(4,i,j,k-1)-u(4,i,j-1,k-1)
     5   +u(4,i-1,j,k)-u(4,i-1,j-1,k)+u(4,i-1,j,k-1)-u(4,i-1,j-1,k-1))
     6   -dz2*p(4,i,j,k)*
     6   (u(4,i,j,k)+u(4,i,j-1,k)-u(4,i,j,k-1)-u(4,i,j-1,k-1)
     7   +u(4,i-1,j,k)+u(4,i-1,j-1,k)-u(4,i-1,j,k-1)-u(4,i-1,j-1,k-1))
     8   -dz2*(u(5,i,j,k)+u(5,i,j-1,k)-u(5,i,j,k-1)-u(5,i,j-1,k-1)
     9   +u(5,i-1,j,k)+u(5,i-1,j-1,k)-u(5,i-1,j,k-1)-u(5,i-1,j-1,k-1))
     9   /p(1,i,j,k)
     1   +rmu*(dx4*(ff(4,i+1,j,k)-2.0*ff(4,i,j,k)+ff(4,i-1,j,k))
     2   +dy4*(ff(4,i,j+1,k)-2.0*ff(4,i,j,k)+ff(4,i,j-1,k))
     3   +dz4*(ff(4,i,j,k+1)-2.0*ff(4,i,j,k)+ff(4,i,j,k-1)))/ff(1,i,j,k)
c
       f(3,i,j,k)=ff(3,i,j,k)+t*((p(11,i,j,k)*p(6,i,j,k)
     1   -p(9,i,j,k)*p(8,i,j,k))/p(1,i,j,k)
     1   -0.0*eud*p(3,i,j,k)-gra*y/ro3)
     2   -dx2*p(2,i,j,k)*
     2   (u(3,i,j,k)+u(3,i,j-1,k)+u(3,i,j,k-1)+u(3,i,j-1,k-1)
     3   -u(3,i-1,j,k)-u(3,i-1,j-1,k)-u(3,i-1,j,k-1)-u(3,i-1,j-1,k-1))
     4   -dy2*p(3,i,j,k)*
     4   (u(3,i,j,k)-u(3,i,j-1,k)+u(3,i,j,k-1)-u(3,i,j-1,k-1)
     5   +u(3,i-1,j,k)-u(3,i-1,j-1,k)+u(3,i-1,j,k-1)-u(3,i-1,j-1,k-1))
     6   -dz2*p(4,i,j,k)*
     6   (u(3,i,j,k)+u(3,i,j-1,k)-u(3,i,j,k-1)-u(3,i,j-1,k-1)
     7   +u(3,i-1,j,k)+u(3,i-1,j-1,k)-u(3,i-1,j,k-1)-u(3,i-1,j-1,k-1))
     8   -dy2*(u(5,i,j,k)-u(5,i,j-1,k)+u(5,i,j,k-1)-u(5,i,j-1,k-1)
     9   +u(5,i-1,j,k)-u(5,i-1,j-1,k)+u(5,i-1,j,k-1)-u(5,i-1,j-1,k-1))
     9   /p(1,i,j,k)
     1   +rmu*(dx4*(ff(3,i+1,j,k)-2.0*ff(3,i,j,k)+ff(3,i-1,j,k))
     2   +dy4*(ff(3,i,j+1,k)-2.0*ff(3,i,j,k)+ff(3,i,j-1,k))
     3   +dz4*(ff(3,i,j,k+1)-2.0*ff(3,i,j,k)+ff(3,i,j,k-1)))/ff(1,i,j,k)
c
       uxd=amax1((x-4.0*ar1)/(4.0*ar1),0.0)
       uxd=amax1((x-xta)/xtd,0.0)
       uxd=amin1(uxd,1.0)
       uxd=amin1(uxd*p(2,i,j,k),0.0)
c
       f(2,i,j,k)=ff(2,i,j,k)+t*((p(10,i,j,k)*p(8,i,j,k)
     1   -p(11,i,j,k)*p(7,i,j,k))/p(1,i,j,k)
     1   -eud*uxd-gra*x/ro3)
     2   -dx2*p(2,i,j,k)*
     2   (u(2,i,j,k)+u(2,i,j-1,k)+u(2,i,j,k-1)+u(2,i,j-1,k-1)
     3   -u(2,i-1,j,k)-u(2,i-1,j-1,k)-u(2,i-1,j,k-1)-u(2,i-1,j-1,k-1))
     4   -dy2*p(3,i,j,k)*
     4   (u(2,i,j,k)-u(2,i,j-1,k)+u(2,i,j,k-1)-u(2,i,j-1,k-1)
     5   +u(2,i-1,j,k)-u(2,i-1,j-1,k)+u(2,i-1,j,k-1)-u(2,i-1,j-1,k-1))
     6   -dz2*p(4,i,j,k)*
     6   (u(2,i,j,k)+u(2,i,j-1,k)-u(2,i,j,k-1)-u(2,i,j-1,k-1)
     7   +u(2,i-1,j,k)+u(2,i-1,j-1,k)-u(2,i-1,j,k-1)-u(2,i-1,j-1,k-1))
     8   -dx2*(u(5,i,j,k)+u(5,i,j-1,k)+u(5,i,j,k-1)+u(5,i,j-1,k-1)
     9   -u(5,i-1,j,k)-u(5,i-1,j-1,k)-u(5,i-1,j,k-1)-u(5,i-1,j-1,k-1))
     9   /p(1,i,j,k)
     1   +rmu*(dx4*(ff(2,i+1,j,k)-2.0*ff(2,i,j,k)+ff(2,i-1,j,k))
     2   +dy4*(ff(2,i,j+1,k)-2.0*ff(2,i,j,k)+ff(2,i,j-1,k))
     3   +dz4*(ff(2,i,j,k+1)-2.0*ff(2,i,j,k)+ff(2,i,j,k-1)))/ff(1,i,j,k)
c
       f(1,i,j,k)=ff(1,i,j,k)
     1   -dx2*(u(2,i,j,k)*u(1,i,j,k)+u(2,i,j-1,k)*u(1,i,j-1,k)
     2   +u(2,i,j,k-1)*u(1,i,j,k-1)+u(2,i,j-1,k-1)*u(1,i,j-1,k-1)
     3   -u(2,i-1,j,k)*u(1,i-1,j,k)-u(2,i-1,j-1,k)*u(1,i-1,j-1,k)
     4   -u(2,i-1,j,k-1)*u(1,i-1,j,k-1)
     4   -u(2,i-1,j-1,k-1)*u(1,i-1,j-1,k-1))
     5   -dy2*(u(3,i,j,k)*u(1,i,j,k)-u(3,i,j-1,k)*u(1,i,j-1,k)
     6   +u(3,i,j,k-1)*u(1,i,j,k-1)-u(3,i,j-1,k-1)*u(1,i,j-1,k-1)
     7   +u(3,i-1,j,k)*u(1,i-1,j,k)-u(3,i-1,j-1,k)*u(1,i-1,j-1,k)
     8   +u(3,i-1,j,k-1)*u(1,i-1,j,k-1)
     8   -u(3,i-1,j-1,k-1)*u(1,i-1,j-1,k-1))
     1   -dz2*(u(4,i,j,k)*u(1,i,j,k)+u(4,i,j-1,k)*u(1,i,j-1,k)
     2   -u(4,i,j,k-1)*u(1,i,j,k-1)-u(4,i,j-1,k-1)*u(1,i,j-1,k-1)
     3   +u(4,i-1,j,k)*u(1,i-1,j,k)+u(4,i-1,j-1,k)*u(1,i-1,j-1,k)
     4   -u(4,i-1,j,k-1)*u(1,i-1,j,k-1)
     4   -u(4,i-1,j-1,k-1)*u(1,i-1,j-1,k-1))
     5   +dfu*(dx4*(ff(1,i+1,j,k)-2.0*ff(1,i,j,k)+ff(1,i-1,j,k))
     6   +dy4*(ff(1,i,j+1,k)-2.0*ff(1,i,j,k)+ff(1,i,j-1,k))
     7   +dz4*(ff(1,i,j,k+1)-2.0*ff(1,i,j,k)+ff(1,i,j,k-1)))
c
  801  continue
c
c      do 820 i=2,nx1
       do 820 i=is1,ie1
       x=hhx(iss+i)+xx4
       ro2=x*x+y*y+z*z
       x2=amax1(ro2/ar2-1.0,0.0)
       x2=aru*x2*x2
       x2=x2/(x2+1.0)
       f(8,i,j,k)=f(8,i,j,k)*x2+ff(8,i,j,k)*(1.0-x2)
       f(7,i,j,k)=f(7,i,j,k)*x2+ff(7,i,j,k)*(1.0-x2)
       f(6,i,j,k)=f(6,i,j,k)*x2+ff(6,i,j,k)*(1.0-x2)
       f(5,i,j,k)=f(5,i,j,k)*x2+ff(5,i,j,k)*(1.0-x2)
       f(4,i,j,k)=f(4,i,j,k)*x2+ff(4,i,j,k)*(1.0-x2)
       f(3,i,j,k)=f(3,i,j,k)*x2+ff(3,i,j,k)*(1.0-x2)
       f(2,i,j,k)=f(2,i,j,k)*x2+ff(2,i,j,k)*(1.0-x2)
       f(1,i,j,k)=f(1,i,j,k)*x2+ff(1,i,j,k)*(1.0-x2)
  820  continue
c
c      do 824 i=2,nx1
       do 824 i=is1,ie1
       x5=f(5,i,j,k)
       x1=f(1,i,j,k)
       y5=amax1(x5,pr02)
       y1=amax1(x1,ro02)
       x2=f(2,i,j,k)*f(2,i,j,k)+f(3,i,j,k)*f(3,i,j,k)
     1    +f(4,i,j,k)*f(4,i,j,k)
       vmax=amax1(vmax,x2)
       x=amax1(x1,0.0)/y1
       f(5,i,j,k)=y5
       f(4,i,j,k)=x*f(4,i,j,k)
       f(3,i,j,k)=x*f(3,i,j,k)
       f(2,i,j,k)=x*f(2,i,j,k)
       f(1,i,j,k)=y1
  824  continue
   80  continue
c
  190  continue
c
CC MPI START
       call mpi_allreduce(vmax,vmax1,1,mpi_real,mpi_max,
     *                    mpi_comm_world,ier)
       vmax=vmax1
CC MPI END
c
c
       time=time+t
c
       if(iia.ne.1.and.iiq.ne.iiq0) go to 399
       t1=t
       dx1=0.25*t1/hx
       dy1=0.25*t1/hy
       dz1=0.25*t1/hz
       dx3=t1/(hx*hx)
       dy3=t1/(hy*hy)
       dz3=t1/(hz*hz)
  399  continue
c
       if(vmax.gt.1.0) go to 402
c
       if(iiq.ne.iiq0) go to 401
       iiq=0
c
  401  if(iip.ne.iip0) go to 403
       iip=0
c
c      th=float(iia)+ths0*float(last)
c      th=orid*th*pi/180.0
c      th=orid*pi/180.0
c      ith=iia+last*(iii-1)
c      th1=180.0*float(ith)/float(last*4)
c      th2=orid+th1
c      if(th2.gt.360.0) th2=th2-360.0
c      th=th2*pi/180.0
c      bisy=bis*cos(th)
c      bisz=bis*sin(th)
       ith=iia+last*(iii-1)
       th1=180.0*float(ith)/float(last*4)
       th2=orid+th1
c
c      iino=iia+last*(iii-1)+iinob
c      jjno=ifix(float(iino)/ttno)+1
       x2=2.0
       x1=2.0*x2
       nx11=nx1-2*nxp
c      roo1=fdd(2,jjno)
c      vswx=-fdd(3,jjno)
c      pro1=fdd(4,jjno)
c      bisy=-fdd(5,jjno)
c      bisz=fdd(6,jjno)
       roo1=ro01
       vswx=vsw
       pro1=pr01
c      bisy=0.0
c      bisz=bis
c
c      do 102 k=1,nz2
       do 102 k=ks,ke
c      do 102 j=1,ny2
       do 102 j=js,je
c      do 102 i=1,nx2
       do 102 i=is,ie
       x=(2.0*hhx(iss+i)/hx-float(2+nx11))/float(nx11)
       x=(x1+1.0)*x+x1
       x=-x
       x=amax1(x,0.0)
       x=amin1(x,1.0)
       x=x*x
c      f(8,i,j,k)=(1.0-x)*f(8,i,j,k)+bis*x
       f(8,i,j,k)=(1.0-x)*f(8,i,j,k)+bisz*x
       f(7,i,j,k)=(1.0-x)*f(7,i,j,k)+bisy*x
       f(6,i,j,k)=(1.0-x)*f(6,i,j,k)
       f(5,i,j,k)=(1.0-x)*f(5,i,j,k)+pro1*x
       f(4,i,j,k)=(1.0-x)*f(4,i,j,k)
       f(3,i,j,k)=(1.0-x)*f(3,i,j,k)
       f(2,i,j,k)=(1.0-x)*f(2,i,j,k)+vswx*x
       f(1,i,j,k)=(1.0-x)*f(1,i,j,k)+roo1*x
  102  continue
c
c      do 104 k=1,nz1
       do 104 k=ks,ke1
c      do 104 j=1,ny1
       do 104 j=js,je1
c      do 104 i=1,nx1
       do 104 i=is,ie1
       x=(2.0*hhx(iss+i)/hx-float(1+nx11))/float(nx11)
       x=(x1+1.0)*x+x1
       x=-x
       x=amax1(x,0.0)
       x=amin1(x,1.0)
       x=x*x
c      u(8,i,j,k)=(1.0-x)*u(8,i,j,k)+bis*x
       u(8,i,j,k)=(1.0-x)*u(8,i,j,k)+bisz*x
       u(7,i,j,k)=(1.0-x)*u(7,i,j,k)+bisy*x
       u(6,i,j,k)=(1.0-x)*u(6,i,j,k)
       u(5,i,j,k)=(1.0-x)*u(5,i,j,k)+pro1*x
       u(4,i,j,k)=(1.0-x)*u(4,i,j,k)
       u(3,i,j,k)=(1.0-x)*u(3,i,j,k)
       u(2,i,j,k)=(1.0-x)*u(2,i,j,k)+vswx*x
       u(1,i,j,k)=(1.0-x)*u(1,i,j,k)+roo1*x
  104  continue
c
c
c      inner boundary condition
       call equib8a(ff,pp,nxp,ro01,pr01,rra1,cj,cp)
c
c      do 390 k=2,nz1
       do 390 k=ks1,ke1
c
       xx4=0.5*hx*float(2*nxp-nx1-2)
       z=0.5*hz*float(2*(kss+k)-nz2-1)
c      do 380 j=2,ny1
       do 380 j=js1,je1
c      y=0.5*hy*float(2*(jss+j)-3)
       y=0.5*hy*float(2*(jss+j)-ny2-1)
c
c      do 380 i=2,nx1
       do 380 i=is1,ie1
       x=hhx(iss+i)+xx4
       ro2=x*x+y*y+z*z
       x2=amax1(ro2/ar2-1.0,0.0)
       x2=aru*x2*x2
       x2=x2/(x2+1.0)
       f(8,i,j,k)=f(8,i,j,k)*x2+ff(8,i,j,k)*(1.0-x2)
       f(7,i,j,k)=f(7,i,j,k)*x2+ff(7,i,j,k)*(1.0-x2)
       f(6,i,j,k)=f(6,i,j,k)*x2+ff(6,i,j,k)*(1.0-x2)
       f(5,i,j,k)=f(5,i,j,k)*x2+ff(5,i,j,k)*(1.0-x2)
       f(4,i,j,k)=f(4,i,j,k)*x2+ff(4,i,j,k)*(1.0-x2)
       f(3,i,j,k)=f(3,i,j,k)*x2+ff(3,i,j,k)*(1.0-x2)
       f(2,i,j,k)=f(2,i,j,k)*x2+ff(2,i,j,k)*(1.0-x2)
       f(1,i,j,k)=f(1,i,j,k)*x2+ff(1,i,j,k)*(1.0-x2)
  380  continue
c
  390  continue
c
c
  403  if(iis.ne.iis0) go to 100
       iis=0
  402  continue
       call clock(zt2)
c      if(iia.lt.1024) go to 100
c
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
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)
c        call mpi_gatherv(fgl(1,1,1),mword,mpi_real,fg(1,1,1),
c    1                    recvcount,displs,mpi_real,0,
c    2                    mpi_comm_world,ier)
c
c      call mpi_barrier(mpi_comm_world,ier)
       if(irank.eq.0) then
       do 1732 k=1,nz2
c      write(ntap) fg(1:nx2,1:ny2,k)
 1732  continue
       end if
CC MPI END
  173  continue
c
c
c
CC MPI START
       if(irank.eq.0)
     * write(6,661) iia,last,iii,ith,nin,ntap,itap,th2,th,x2,vmax
CC MPI END
  661  format(1h ,1x,7i10/1x,1p4e15.7)
c
       call clock(zt3)
       zt1=zt1-zt0
       zt2=zt2-zt0
       zt3=zt3-zt0
       zt=zt2-zt1
CC MPI START
       if(irank.eq.0)
     * write(6,404) iia,zt0,zt1,zt2,zt3,zt
CC MPI END
  404  format(1h , i6,1pe12.3,4(0pf12.5))
       zt1=zt3+zt0
c
       if(vmax.gt.1.0) go to 9
c
  100  continue
  300  continue
    9  continue
c      rewind 12
       call mpi_finalize(ierr)
       stop
       end
       subroutine equib7(f,pp,nxp,ro01,pr01,rrat,cj,cp)
c      implicit real*8 (a-h,o-z)
CC MPI START
       include 'mpif.h'
       integer istatus(mpi_status_size)
       parameter (npex=2,npey=2,npez=2)
       parameter (npe=npex*npey*npez,npexy=npex*npey)
c for mpi_gatherv
c      integer recvcount(npe),displs(npe),bound(3,npe)
       integer itable(-1:npex,-1:npey,-1:npez)
       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
       common /para_table/itable
CC MPI END
c      parameter (nx=320,ny= 80,nz=160,nb=8)
c      parameter (nx=500,ny=100,nz=200,nb=8)
c      parameter (nx=318,ny=318,nz=318,nb=8,mb=3)
c      parameter (nx=500,ny=200,nz=200,nb=8,mb=3)
       parameter (nx=510,ny=254,nz=254,nb=8,mb=3)
c      parameter (nx=500,ny=300,nz=300,nb=8)
c      parameter (nx=500,ny=200,nz=400,nb=8)
       parameter (nx1=nx+1,nx2=nx+2,ny1=ny+1,ny2=ny+2)
       parameter (nz1=nz+1,nz2=nz+2)
c
c      dimension f(nx2,ny2,nz2,nb),pp(nx2,ny2,nz2,3)
c      dimension gf(nx2,ny2,nz2,nb),gpp(nx2,ny2,nz2,3)
c      dimension gf(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
       dimension f(nb,0:nxx+1,0:nyy+1,0:nzz+1),
     1           pp(mb,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)
CC MPI END
       dimension cj(1),cp(1)
c      common /blk/f,pp
c      cj 1-r1 2-dr 3-dz 4-b0 5-rm 6-a2 7-aid 8-ar1 9-gam 10-gra
c      cp 1-r0 2-z0 3-ra 4-vo0 5-p0 6-gra 7-07 8-vsw 9-09 10-10
c
c
       aid=cj(7)
       ar1=cj(8)
       gam=cj(9)
       gra=cj(10)
       p0=(gam-1.0)*gra/gam
       bis=cp(11)*1.0e-4
       x1=1.0/(gam-1.0)
       x2=gam*x1
       b0=cj(4)
       ro02=rrat*ro01
       pr02=pr01
c
       ic0=4
       icd=6
       pi=3.1415926
       n1=nx2
       n2=n1*ny2
       n3=n2*nz2
       hx=cp(1)/float(nx1)
       hy=cp(2)/float(ny1)
       hz=cp(3)/float(nz1)
       vsw=cp(9)
       tl=pi*cp(8)/180.0
       cotl=cos(tl)
       sitl=sin(tl)
c
c      do 10 k=1,nz2
       do 10 k=ks,ke
       z=0.5*hz*float(2*(kss+k)-nz2-1)
c      do 101 j=1,ny2
       do 101 j=js,je
c      vr1=vsw
c      y=0.5*hy*float(2*(jss+j)-3)
       y=0.5*hy*float(2*(jss+j)-ny2-1)
c      do 101 i=1,nx2
       do 101 i=is,ie
       vr1=vsw
       x=0.5*hx*float(2*(iss+i)-nx2-1+2*nxp)
c      y=0.5*hy*float(2*(jss+j)-3)
c      z=0.5*hz*float(2*(kss+k)-nz2-1)
       xt=x*cotl+z*sitl
       zt=-x*sitl+z*cotl
       ro2=x*x+y*y+z*z
       ro1=sqrt(ro2)
       ro3=ro1*ro2
       f(1,i,j,k)=1.0/ro3
c      f(1,i,j,k)=ro1**(-x1)
       if(f(1,i,j,k).lt.ro02) f(1,i,j,k)=ro02
       f(5,i,j,k)=1.0e-7*cp(6)/ro2
c      f(5,i,j,k)=p0*ro1**(-x2)
       if(f(5,i,j,k).lt.pr02) f(5,i,j,k)=pr02
c
       bx=-2.0*b0*xt*zt/(ro2*ro2)
       by=-2.0*b0*y*zt/(ro2*ro2)
       bz=b0*(xt*xt+y*y-zt*zt)/(ro2*ro2)
       bx=-3.0*b0*xt*zt/(ro2*ro3)
       by=-3.0*b0*y*zt/(ro2*ro3)
       bz=b0*(xt*xt+y*y-2.0*zt*zt)/(ro2*ro3)
       f(6,i,j,k)=bx*cotl-bz*sitl
       f(7,i,j,k)=by
       f(8,i,j,k)=bx*sitl+bz*cotl
       px1=sqrt(f(6,i,j,k)*f(6,i,j,k)+f(7,i,j,k)*f(7,i,j,k)
     1         +f(8,i,j,k)*f(8,i,j,k))
       if(abs(vsw).lt.1.0e-8) go to 101
       y1=(ro01-cj(1)*(px1/vsw)**2)/ro01
       if(y1.le.0.0) y1=0.0
       y1=vsw*sqrt(y1)
       vr1=amin1(vr1,y1)
       f(2,i,j,k)=vr1
       xx=vr1/vsw
       f(1,i,j,k)=f(1,i,j,k)+(ro01-ro02)*xx
       f(5,i,j,k)=f(5,i,j,k)+(pr01-pr02)*xx
       f(8,i,j,k)=f(8,i,j,k)+bis
  101  continue
   10  continue
c
c
c      current
!hpf$  reflect f
c
CC MPI START
c
       irightx = itable(irankx+1,iranky,irankz)
       ileftx  = itable(irankx-1,iranky,irankz)
       irighty = itable(irankx,iranky+1,irankz)
       ilefty  = itable(irankx,iranky-1,irankz)
       irightz = itable(irankx,iranky,irankz+1)
       ileftz  = itable(irankx,iranky,irankz-1)
c
       ftemp1x=f(:,is,:,:)
       ftemp1y=f(:,:,js,:)
       ftemp1z=f(:,:,:,ks)
       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(:,ie+1,:,:)=ftemp2x
       f(:,:,je+1,:)=ftemp2y
       f(:,:,:,ke+1)=ftemp2z
c
CC MPI END
c
c      do 38 k=1,nz1
       do 38 k=ks,ke1
c
c      do 381 j=1,ny1
       do 381 j=js,je1
c      do 381 i=1,nx1
       do 381 i=is,ie1
       pp(3,i,j,k)=0.25*((f(7,i+1,j+1,k+1)+f(7,i+1,j,k+1)
     1   +f(7,i+1,j+1,k)+f(7,i+1,j,k)
     2   -f(7,i,j+1,k+1)-f(7,i,j,k+1)-f(7,i,j+1,k)-f(7,i,j,k))/hx
     3   -(f(6,i+1,j+1,k+1)-f(6,i+1,j,k+1)+f(6,i+1,j+1,k)-f(6,i+1,j,k)
     4   +f(6,i,j+1,k+1)-f(6,i,j,k+1)+f(6,i,j+1,k)-f(6,i,j,k))/hy)
       pp(2,i,j,k)=0.25*((f(6,i+1,j+1,k+1)+f(6,i+1,j,k+1)
     1   -f(6,i+1,j+1,k)-f(6,i+1,j,k)
     2   +f(6,i,j+1,k+1)+f(6,i,j,k+1)-f(6,i,j+1,k)-f(6,i,j,k))/hz
     3   -(f(8,i+1,j+1,k+1)+f(8,i+1,j,k+1)+f(8,i+1,j+1,k)+f(8,i+1,j,k)
     4   -f(8,i,j+1,k+1)-f(8,i,j,k+1)-f(8,i,j+1,k)-f(8,i,j,k))/hx)
       pp(1,i,j,k)=0.25*((f(8,i+1,j+1,k+1)-f(8,i+1,j,k+1)
     1   +f(8,i+1,j+1,k)-f(8,i+1,j,k)
     2   +f(8,i,j+1,k+1)-f(8,i,j,k+1)+f(8,i,j+1,k)-f(8,i,j,k))/hy
     3   -(f(7,i+1,j+1,k+1)+f(7,i+1,j,k+1)-f(7,i+1,j+1,k)-f(7,i+1,j,k)
     4   +f(7,i,j+1,k+1)+f(7,i,j,k+1)-f(7,i,j+1,k)-f(7,i,j,k))/hz)
  381  continue
   38  continue
c
       return
       end
       subroutine equib8(f,pp,nxp,ro01,pr01,rrat,cj,cp)
c      implicit real*8 (a-h,o-z)
CC MPI START
       include 'mpif.h'
       integer istatus(mpi_status_size)
       parameter (npex=2,npey=2,npez=2)
       parameter (npe=npex*npey*npez,npexy=npex*npey)
c for mpi_gatherv
c      integer recvcount(npe),displs(npe),bound(3,npe)
       integer itable(-1:npex,-1:npey,-1:npez)
       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
       common /para_table/itable
CC MPI END
c      parameter (nx=320,ny= 80,nz=160,nb=8)
c      parameter (nx=500,ny=100,nz=200,nb=8)
c      parameter (nx=318,ny=318,nz=318,nb=8,mb=3)
c      parameter (nx=500,ny=200,nz=200,nb=8,mb=3)
       parameter (nx=510,ny=254,nz=254,nb=8,mb=3)
c      parameter (nx=500,ny=300,nz=300,nb=8)
c      parameter (nx=500,ny=200,nz=400,nb=8)
       parameter (nx1=nx+1,nx2=nx+2,ny1=ny+1,ny2=ny+2)
       parameter (nz1=nz+1,nz2=nz+2)
c
c      dimension f(nx2,ny2,nz2,nb),pp(nx2,ny2,nz2,3)
c      dimension gf(nx2,ny2,nz2,nb),gpp(nx2,ny2,nz2,3)
c      dimension gf(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
       dimension f(nb,0:nxx+1,0:nyy+1,0:nzz+1),
     1           pp(mb,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)
CC MPI END
       dimension cj(1),cp(1)
c      common /blk/f,pp
c      cj 1-r1 2-dr 3-dz 4-b0 5-rm 6-a2 7-aid 8-ar1 9-gam 10-gra
c      cp 1-r0 2-z0 3-ra 4-vo0 5-p0 6-gra 7-07 8-vsw 9-09 10-10
c
c
       aid=cj(7)
       ar1=cj(8)
       gam=cj(9)
       gra=cj(10)
       p0=(gam-1.0)*gra/gam
       bis=cp(11)*1.0e-4
       x1=1.0/(gam-1.0)
       x2=gam*x1
       b0=cj(4)
       ro02=rrat*ro01
       pr02=pr01
c
       ic0=4
       icd=6
       pi=3.1415926
       n1=nx2
       n2=n1*ny2
       n3=n2*nz2
       hx=cp(1)/float(nx1)
       hy=cp(2)/float(ny1)
       hz=cp(3)/float(nz1)
       vsw=cp(9)
       tl=pi*cp(8)/180.0
       cotl=cos(tl)
       sitl=sin(tl)
       alp=cj(1)
       bet=cj(2)
       alp1=1.0/sqrt(alp)
       xm=(2.0*alp*b0*b0/(ro01*vsw*vsw))**(1.0/6.0)
       xm=-xm
c
c      do 10 k=1,nz2
       do 10 k=ks,ke
       z=0.5*hz*float(2*(kss+k)-nz2-1)
c      do 101 j=1,ny2
       do 101 j=js,je
c      vr1=vsw
c      y=0.5*hy*float(2*(jss+j)-3)
       y=0.5*hy*float(2*(jss+j)-ny2-1)
c      do 101 i=1,nx2
       do 101 i=is,ie
       vr1=vsw
       x=0.5*hx*float(2*(iss+i)-nx2-1+2*nxp)
c      y=0.5*hy*float(2*(jss+j)-3)
c      z=0.5*hz*float(2*(kss+k)-nz2-1)
       xt=x*cotl+z*sitl
       zt=-x*sitl+z*cotl
       ro2=x*x+y*y+z*z
       ro1=sqrt(ro2)
       ro3=ro1*ro2
       f(1,i,j,k)=1.0/ro3
c      f(1,i,j,k)=ro1**(-x1)
       if(f(1,i,j,k).lt.ro02) f(1,i,j,k)=ro02
       f(5,i,j,k)=1.0e-7*cp(6)/ro2
c      f(5,i,j,k)=p0*ro1**(-x2)
       if(f(5,i,j,k).lt.pr02) f(5,i,j,k)=pr02
c
       bx=-2.0*b0*xt*zt/(ro2*ro2)
       by=-2.0*b0*y*zt/(ro2*ro2)
       bz=b0*(xt*xt+y*y-zt*zt)/(ro2*ro2)
       bx=-3.0*b0*xt*zt/(ro2*ro3)
       by=-3.0*b0*y*zt/(ro2*ro3)
       bz=b0*(xt*xt+y*y-2.0*zt*zt)/(ro2*ro3)
       f(6,i,j,k)=bx*cotl-bz*sitl
       f(7,i,j,k)=by
       f(8,i,j,k)=bx*sitl+bz*cotl
  101  continue
   10  continue
c
c      write(6,372) nx2,ny2,nz2,nxp,
c    1              b0,cotl,sitl,ro02,pr02,rrat,ro01,pr01
c 372  format(1h ,4i4/8(1pe10.3))
c
c      current
!hpf$  reflect f
c
CC MPI START
c
       irightx = itable(irankx+1,iranky,irankz)
       ileftx  = itable(irankx-1,iranky,irankz)
       irighty = itable(irankx,iranky+1,irankz)
       ilefty  = itable(irankx,iranky-1,irankz)
       irightz = itable(irankx,iranky,irankz+1)
       ileftz  = itable(irankx,iranky,irankz-1)
c
       ftemp1x=f(:,is,:,:)
       ftemp1y=f(:,:,js,:)
       ftemp1z=f(:,:,:,ks)
       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(:,ie+1,:,:)=ftemp2x
       f(:,:,je+1,:)=ftemp2y
       f(:,:,:,ke+1)=ftemp2z
c
CC MPI END
c
c      do 38 k=1,nz1
       do 38 k=ks,ke1
c      do 381 j=1,ny1
       do 381 j=js,je1
c      do 381 i=1,nx1
       do 381 i=is,ie1
       pp(3,i,j,k)=0.25*((f(7,i+1,j+1,k+1)+f(7,i+1,j,k+1)
     1   +f(7,i+1,j+1,k)+f(7,i+1,j,k)
     2   -f(7,i,j+1,k+1)-f(7,i,j,k+1)-f(7,i,j+1,k)-f(7,i,j,k))/hx
     3   -(f(6,i+1,j+1,k+1)-f(6,i+1,j,k+1)+f(6,i+1,j+1,k)-f(6,i+1,j,k)
     4   +f(6,i,j+1,k+1)-f(6,i,j,k+1)+f(6,i,j+1,k)-f(6,i,j,k))/hy)
       pp(2,i,j,k)=0.25*((f(6,i+1,j+1,k+1)+f(6,i+1,j,k+1)
     1   -f(6,i+1,j+1,k)-f(6,i+1,j,k)
     2   +f(6,i,j+1,k+1)+f(6,i,j,k+1)-f(6,i,j+1,k)-f(6,i,j,k))/hz
     3   -(f(8,i+1,j+1,k+1)+f(8,i+1,j,k+1)+f(8,i+1,j+1,k)+f(8,i+1,j,k)
     4   -f(8,i,j+1,k+1)-f(8,i,j,k+1)-f(8,i,j+1,k)-f(8,i,j,k))/hx)
       pp(1,i,j,k)=0.25*((f(8,i+1,j+1,k+1)-f(8,i+1,j,k+1)
     1   +f(8,i+1,j+1,k)-f(8,i+1,j,k)
     2   +f(8,i,j+1,k+1)-f(8,i,j,k+1)+f(8,i,j+1,k)-f(8,i,j,k))/hy
     3   -(f(7,i+1,j+1,k+1)+f(7,i+1,j,k+1)-f(7,i+1,j+1,k)-f(7,i+1,j,k)
     4   +f(7,i,j+1,k+1)+f(7,i,j,k+1)-f(7,i,j+1,k)-f(7,i,j,k))/hz)
  381  continue
   38  continue
c
c
c      do 40 k=1,nz2
       do 40 k=ks,ke
       z=0.5*hz*float(2*(kss+k)-nz2-1)
c      do 401 j=1,ny2
       do 401 j=js,je
c      vr1=vsw
c      y=0.5*hy*float(2*(jss+j)-3)
       y=0.5*hy*float(2*(jss+j)-ny2-1)
c      do 401 i=1,nx2
       do 401 i=is,ie
       x=0.5*hx*float(2*(iss+i)-nx2-1+2*nxp)
c      y=0.5*hy*float(2*(jss+j)-3)
c      z=0.5*hz*float(2*(kss+k)-nz2-1)
       x1m=x-2.0*xm
c      x1m=x+2.0*xm
       xt=x*cotl+z*sitl
       zt=-x*sitl+z*cotl
       xtm=x1m*cotl-z*sitl
       ztm=x1m*sitl+z*cotl
       ro2=x*x+y*y+z*z
       ro1=sqrt(ro2)
       ro3=ro1*ro2
       ro2m=x1m*x1m+y*y+z*z
       ro1m=sqrt(ro2m)
       ro3m=ro1m*ro2m
       f(1,i,j,k)=1.0/ro3
c      f(1,i,j,k)=ro1**(-x1)
       if(f(1,i,j,k).lt.ro02) f(1,i,j,k)=ro02
       f(5,i,j,k)=1.0e-7*cp(6)/ro2
c      f(5,i,j,k)=p0*ro1**(-x2)
       if(f(5,i,j,k).lt.pr02) f(5,i,j,k)=pr02
c
       bx1=-2.0*b0*xt*zt/(ro2*ro2)
       by1=-2.0*b0*y*zt/(ro2*ro2)
       bz1=b0*(xt*xt+y*y-zt*zt)/(ro2*ro2)
       bx1=-3.0*b0*xt*zt/(ro2*ro3)
       by1=-3.0*b0*y*zt/(ro2*ro3)
       bz1=b0*(xt*xt+y*y-2.0*zt*zt)/(ro2*ro3)
       bx2=-2.0*b0*xtm*ztm/(ro2m*ro2m)
       by2=-2.0*b0*y*ztm/(ro2m*ro2m)
       bz2=b0*(xtm*xtm+y*y-ztm*ztm)/(ro2m*ro2m)
       bx2=-3.0*b0*xtm*ztm/(ro2m*ro3m)
       by2=-3.0*b0*y*ztm/(ro2m*ro3m)
       bz2=b0*(xtm*xtm+y*y-2.0*ztm*ztm)/(ro2m*ro3m)
c
       bx11=bx1*cotl-bz1*sitl
       by11=by1
       bz11=bx1*sitl+bz1*cotl
       bx22=bx2*cotl+bz2*sitl
       by22=by2
       bz22=-bx2*sitl+bz2*cotl
c      bx1=-3.0*b0*x*z/(ro2*ro3)
c      by1=-3.0*b0*y*z/(ro2*ro3)
c      bz1=b0*(x*x+y*y-2.0*z*z)/(ro2*ro3)
c      bx2=-3.0*b0*x1m*z/(ro2m*ro3m)
c      by2=-3.0*b0*y*z/(ro2m*ro3m)
c      bz2=b0*(x1m*x1m+y*y-2.0*z*z)/(ro2m*ro3m)
       f(6,i,j,k)=bx11+bx22
       f(7,i,j,k)=by11+by22
       f(8,i,j,k)=bz11+bz22
       if(x.lt.xm) f(6,i,j,k)=0.0
       if(x.lt.xm) f(7,i,j,k)=0.0
       if(x.lt.xm) f(8,i,j,k)=0.0
       y1=alp1+(1.0-alp1)*(x-xm)/(bet-1.0)/xm
       if(x.gt.xm) y1=0.0
       y1=amax1(y1,0.0)
       y1=amin1(y1,1.0)
       vr1=vsw*y1
c
       px1=sqrt(f(6,i,j,k)*f(6,i,j,k)+f(7,i,j,k)*f(7,i,j,k)
     1         +f(8,i,j,k)*f(8,i,j,k))
       if(abs(vsw).lt.1.0e-8) go to 401
c      y1=(ro01-cj(1)*(px1/vsw)**2)/ro01
c      if(y1.le.0.0) y1=0.0
c      y1=vsw*sqrt(y1)
c      vr1=amin1(vr1,y1)
       f(2,i,j,k)=vr1
       xx=vr1/vsw
       f(5,i,j,k)=f(5,i,j,k)+(pr01-pr02)*xx
       f(8,i,j,k)=f(8,i,j,k)+bis*y1
  401  continue
   40  continue
c
       return
       end
       subroutine equib8a(u,pp,nxp,ro01,pr01,rrat,cj,cp)
c      implicit real*8 (a-h,o-z)
CC MPI START
       include 'mpif.h'
       integer istatus(mpi_status_size)
       parameter (npex=2,npey=2,npez=2)
       parameter (npe=npex*npey*npez,npexy=npex*npey)
c for mpi_gatherv
c      integer recvcount(npe),displs(npe),bound(3,npe)
       integer itable(-1:npex,-1:npey,-1:npez)
       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
       common /para_table/itable
CC MPI END
c      parameter (nx=320,ny= 80,nz=160,nb=8)
c      parameter (nx=500,ny=100,nz=200,nb=8)
c      parameter (nx=318,ny=318,nz=318,nb=8,mb=3)
c      parameter (nx=500,ny=200,nz=200,nb=8,mb=3)
       parameter (nx=510,ny=254,nz=254,nb=8,mb=3)
c      parameter (nx=500,ny=300,nz=300,nb=8)
c      parameter (nx=500,ny=200,nz=400,nb=8)
       parameter (nx1=nx+1,nx2=nx+2,ny1=ny+1,ny2=ny+2)
       parameter (nz1=nz+1,nz2=nz+2)
c
c      dimension u(nx2,ny2,nz2,nb),pp(nx2,ny2,nz2,3)
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
       dimension u(nb,0:nxx+1,0:nyy+1,0:nzz+1),
     1           pp(mb,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)
CC MPI END
       dimension cj(1),cp(1)
c      cj 1-r1 2-dr 3-dz 4-b0 5-rm 6-a2 7-aid 8-ar1 9-gam 10-gra
c      cp 1-r0 2-z0 3-ra 4-vo0 5-p0 6-gra 7-07 8-vsw 9-09 10-10
c
c
       aid=cj(7)
       ar1=cj(8)
       gam=cj(9)
       gra=cj(10)
       p0=(gam-1.0)*gra/gam
       bis=cp(11)*1.0e-4
       x1=1.0/(gam-1.0)
       x2=gam*x1
       b0=cj(4)
       ro02=rrat*ro01
       pr02=pr01
c
       ic0=4
       icd=6
       pi=3.1415926
       n1=nx2
       n2=n1*ny2
       n3=n2*nz2
       hx=cp(1)/float(nx1)
       hy=cp(2)/float(ny1)
       hz=cp(3)/float(nz1)
       vsw=cp(9)
       tl=pi*cp(8)/180.0
       cotl=cos(tl)
       sitl=sin(tl)
       alp=cj(1)
       bet=cj(2)
       alp1=1.0/sqrt(alp)
       xm=(2.0*alp*b0*b0/(ro01*vsw*vsw))**(1.0/6.0)
       xm=-xm
c
c      do 10 k=1,nz2
       do 10 k=ks,ke
       z=0.5*hz*float(2*(kss+k)-nz2-1)
c      do 101 j=1,ny2
       do 101 j=js,je
c      vr1=vsw
c      y=0.5*hy*float(2*(jss+j)-3)
       y=0.5*hy*float(2*(jss+j)-ny2-1)
c      do 101 i=1,nx2
       do 101 i=is,ie
       vr1=vsw
       x=0.5*hx*float(2*(iss+i)-nx2-1+2*nxp)
c      y=0.5*hy*float(2*(jss+j)-3)
c      z=0.5*hz*float(2*(kss+k)-nz2-1)
       xt=x*cotl+z*sitl
       zt=-x*sitl+z*cotl
       ro2=x*x+y*y+z*z
       ro1=sqrt(ro2)
       ro3=ro1*ro2
       u(1,i,j,k)=1.0/ro3
c      u(1,i,j,k)=ro1**(-x1)
       if(u(1,i,j,k).lt.ro02) u(1,i,j,k)=ro02
       u(5,i,j,k)=1.0e-7*cp(6)/ro2
c      u(5,i,j,k)=p0*ro1**(-x2)
       if(u(5,i,j,k).lt.pr02) u(5,i,j,k)=pr02
c
       bx=-2.0*b0*xt*zt/(ro2*ro2)
       by=-2.0*b0*y*zt/(ro2*ro2)
       bz=b0*(xt*xt+y*y-zt*zt)/(ro2*ro2)
       bx=-3.0*b0*xt*zt/(ro2*ro3)
       by=-3.0*b0*y*zt/(ro2*ro3)
       bz=b0*(xt*xt+y*y-2.0*zt*zt)/(ro2*ro3)
       u(6,i,j,k)=bx*cotl-bz*sitl
       u(7,i,j,k)=by
       u(8,i,j,k)=bx*sitl+bz*cotl
  101  continue
   10  continue
c
c      write(6,372) nx2,ny2,nz2,nxp,
c    1              b0,cotl,sitl,ro02,pr02,rrat,ro01,pr01
c 372  format(1h ,4i4/8(1pe10.3))
c
c      current
!hpf$  reflect u
c
CC MPI START
c
       irightx = itable(irankx+1,iranky,irankz)
       ileftx  = itable(irankx-1,iranky,irankz)
       irighty = itable(irankx,iranky+1,irankz)
       ilefty  = itable(irankx,iranky-1,irankz)
       irightz = itable(irankx,iranky,irankz+1)
       ileftz  = itable(irankx,iranky,irankz-1)
c
       ftemp1x=u(:,is,:,:)
       ftemp1y=u(:,:,js,:)
       ftemp1z=u(:,:,:,ks)
       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
       u(:,ie+1,:,:)=ftemp2x
       u(:,:,je+1,:)=ftemp2y
       u(:,:,:,ke+1)=ftemp2z
c
CC MPI END
c
c
c      do 38 k=1,nz1
       do 38 k=ks,ke1
c      do 381 j=1,ny1
       do 381 j=js,je1
c      do 381 i=1,nx1
       do 381 i=is,ie1
       pp(3,i,j,k)=0.25*((u(7,i+1,j+1,k+1)+u(7,i+1,j,k+1)
     1   +u(7,i+1,j+1,k)+u(7,i+1,j,k)
     2   -u(7,i,j+1,k+1)-u(7,i,j,k+1)-u(7,i,j+1,k)-u(7,i,j,k))/hx
     3   -(u(6,i+1,j+1,k+1)-u(6,i+1,j,k+1)+u(6,i+1,j+1,k)-u(6,i+1,j,k)
     4   +u(6,i,j+1,k+1)-u(6,i,j,k+1)+u(6,i,j+1,k)-u(6,i,j,k))/hy)
       pp(2,i,j,k)=0.25*((u(6,i+1,j+1,k+1)+u(6,i+1,j,k+1)
     1   -u(6,i+1,j+1,k)-u(6,i+1,j,k)
     2   +u(6,i,j+1,k+1)+u(6,i,j,k+1)-u(6,i,j+1,k)-u(6,i,j,k))/hz
     3   -(u(8,i+1,j+1,k+1)+u(8,i+1,j,k+1)+u(8,i+1,j+1,k)+u(8,i+1,j,k)
     4   -u(8,i,j+1,k+1)-u(8,i,j,k+1)-u(8,i,j+1,k)-u(8,i,j,k))/hx)
       pp(1,i,j,k)=0.25*((u(8,i+1,j+1,k+1)-u(8,i+1,j,k+1)
     1   +u(8,i+1,j+1,k)-u(8,i+1,j,k)
     2   +u(8,i,j+1,k+1)-u(8,i,j,k+1)+u(8,i,j+1,k)-u(8,i,j,k))/hy
     3   -(u(7,i+1,j+1,k+1)+u(7,i+1,j,k+1)-u(7,i+1,j+1,k)-u(7,i+1,j,k)
     4   +u(7,i,j+1,k+1)+u(7,i,j,k+1)-u(7,i,j+1,k)-u(7,i,j,k))/hz)
  381  continue
   38  continue
c
       return
       end
       subroutine clock(ti)
       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