c
       program ognca7
c
c      file name  /zeartha2/eartha14.f
c      eartha14.f  modified leap-frog scheme
c      ogtac.cntl(ogncb77) from ognc.cntl(ogncd75)
c      3d mhd simulation of 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)
       common /para_info/ks,ks1,ke,ke1,kss,irank,isize
c for mpi_gatherv
       parameter (npe=2)
       integer recvcount(npe),displs(npe)
       integer rcvfg(npe),topfg(npe)
CC MPI END
!hpf$  processors pe(npe)
c      parameter (nx=320,ny= 80,nz=160,nxp=100)
       parameter (nx=500,ny=100,nz=200,nxp=190)
c      parameter (nx=640,ny=160,nz=320,nxp=200)
c      parameter (nx=800,ny=200,nz=400,nxp=250)
c      parameter (nx=800,ny=200,nz=478,nxp=250)
c      parameter (nx=800,ny=200,nz=670,nxp=250)
c      parameter (nx=1000,ny= 500,nz=1118,nxp=250)
c      parameter (nx=1000,ny=1000,nz=1118,nxp=250)
c      parameter (nx=1000,ny=1000,nz=1258,nxp=250)
c      parameter (nx=1678,ny= 558,nz=1118,nxp=280)
c      parameter (nx=2238,ny= 558,nz=1118,nxp=820)
       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)
c      parameter (nb=8,nbb=11,n4=n3*nb,n5=n3*nbb)
       parameter (nb=8,nbb=11)
c      parameter (last=7680,nil=2,lll=18,ibig=0)
c      parameter (iiq0=8,iip0= 32,iis0=7680,thx=4.00)
c      parameter (last=1280,nil=2,lll=18,ibig=0)
c      parameter (last=2560,nil=2,lll=18,ibig=0)
c      parameter (last=5120,nil=2,lll=18,ibig=0)
c      parameter (last=512,nil=2,lll=18,ibig=0)
c      parameter (last=128,nil=2,lll=18,ibig=0)
       parameter (last= 32,nil=2,lll=18,ibig=0)
c      parameter (last=128,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=3840,thx=2.50)
c      parameter (iiq0=8,iip0= 32,iis0=3840,thx=2.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=512,thx=0.50)
c      parameter (iiq0=8,iip0= 32,iis0= 32,thx=0.50)
       parameter (iiq0=8,iip0=  8,iis0=  8,thx=0.50)
c      parameter (iiq0=8,iip0= 32,iis0=3840,thx=1.50)
c      parameter (iiq0=8,iip0= 32,iis0=1920,thx=2.00)
c      parameter (iiq0=8,iip0= 32,iis0=1920,thx=1.50)
c      parameter (iiq0=8,iip0= 32,iis0=1920,thx=1.00)
c      parameter (iiq0=8,iip0= 32,iis0=1920,thx=0.75)
c      parameter (iiq0=8,iip0= 32,iis0=1920,thx=0.50)
c      parameter (iiq0=8,iip0= 32,iis0=1920,thx=0.33333333)
c      parameter (iiq0=8,iip0= 32,iis0=1920,thx=0.25)
c      parameter (iiq0=8,iip0= 32,iis0=1920,thx=0.16666667)
c      parameter (rmuu=0.0020,eud0=0.0000,ro01=5.0e-4,pr01=3.56e-8)
c      parameter (rmuu=0.0020,eud0=0.0000,ro01=5.0e-4,pr01=3.56e-8)
       parameter (rmuu=0.0001,eud0=0.0000,ro01=5.0e-4,pr01=3.56e-8)
c      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)
c      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.002,rrat=0.1,rra1=0.1)
       parameter (eatt=0.0001,rrat=0.02,rra1=0.02)
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=270.750,ths0=0.0,kk8=k3*8)
       parameter (orid= 60.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)
       parameter (mfd=6,nfd=1441,iinob=last*840,ttno=512.0)
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
CC MPI START
       parameter(nzz=(nz2-1)/npe+1)
       dimension f(nx2,ny2,0:nzz+1,nb),u(nx2,ny2,0:nzz+1,nb),
     1           ff(nx2,ny2,0:nzz+1,nb),p(nx2,ny2,0:nzz+1,nbb),
     2           pp(nx2,ny2,0:nzz+1,3),
     3           fq(nx2,nb,0:nzz+1)
c for all_gather
c      dimension fg(nx2,ny2,nz2,nb)
       dimension fg(nx2,ny2,nz2),
     1           fqg(nx2,nb,nz2),fqg1(nx2,nb,nz2)
CC MPI END
       dimension fdd(mfd,nfd),gfdd(mfd,nfd)
       dimension hhx(n1),fbb(16),
     1           cj(10),cp(11),v(n2)
c
!hpf$  distribute f(*,*,block,*) onto pe
!hpf$  distribute u(*,*,block,*) onto pe
!hpf$  distribute pp(*,*,block,*) onto pe
!hpf$  distribute ff(*,*,block,*) onto pe
!hpf$  distribute p(*,*,block,*) onto pe
!hpf$  shadow     f(0,0,1:1,0)
!hpf$  shadow     u(0,0,1:1,0)
!hpf$  shadow     pp(0,0,1:1,0)
!hpf$  shadow     ff(0,0,1:1,0)
!hpf$  shadow     p(0,0,1:1,0)
!hpf$  asyncid id1
c      !xocl global gf,gu,gff,gpp,gfdd
c      equivalence (gf,f),(gff,ff),(gu,u),(gpp,pp)
c      equivalence (gfdd,fdd)
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
c      cp      xl   yl   zl   ra  vo0  p0   gra  07  vsw  ibr ibz
c      data cp/161.0,41.0,51.0,4.0,6.81,2.68,1.35,7.0,.044,9.0,-0.0/
c      data cp/128.4,32.4,64.4,4.0,0.00,2.68,1.35,7.0,.044,0.0,-0.0/
c      data cp/160.5,40.5,80.5,4.0,0.00,2.68,1.35,7.0,.044,0.0,-0.0/
       data cp/250.5,50.5,100.5,4.0,0.00,2.68,1.35,7.0,.044,0.0,-0.0/
c      data cp/223.9,55.9,111.9,3.5,0.00,2.68,1.35,7.0,.044,0.0,-0.0/
c      data cp/83.95,27.95,55.95,3.5,0.00,2.68,1.35,7.0,.044,0.0,-0.0/
c
c
       call clock(zt0)
CC MPI START
       call mpi_init(ier)
       call mpi_comm_rank(mpi_comm_world,irank,ier)
       call mpi_comm_size(mpi_comm_world,isize,ier)
c
       kk=nz2/isize
       kmod=mod(nz2,isize)
c
       ks=1
       kss=irank*kk+min(kmod,irank)
c
       if (irank.lt.kmod) kk=kk+1
       ke=ks+kk-1
       ks1=ks
       ke1=ke
       if (irank.eq.0) ks1=2
       if (irank.eq.isize-1) ke1=ke-1
c
       nword=(ke-ks+1)*nx2*ny2
       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
       nwfg=(ke-ks+1)*nx2*nb
       do i=1,isize
       rcvfg(i)=recvcount(i)*nb/ny2
       topfg(i)=displs(i)*nb/ny2
       end do
c
c
CC MPI END
c
c
c      rewind 12
c
CC MPI START
c      if(irank.eq.0)
c    1 open(9,file='ac990319.data',access='sequential',
c    2         form='unformatted')
c      if(irank.eq.0)
c    1 open(10,file='fb190840.data',access='sequential',
c    2         form='unformatted')
c
c      if(irank.eq.0)
c    1 open(11,file='/home/dpfs/usr6/a41456a/fa168001.data',
c    2         access='sequential',form='unformatted')
CC MPI END
c
c
       n6=n3*18
       n5=n3*nbb
       n4=n3*nb
c
c
cc      read(9) gfdd
       do 660 i=1,nfd
c      read(9,662) (gfdd(j,i),j=1,mfd)
       do 664 j=1,mfd
cc      fdd(j,i)=gfdd(j,i)
  664  continue
  662  format(1h ,6(1pe12.4))
  660  continue
c
c      do 300 iii=1,5
c      do 300 iii=1,30
c      do 300 iii=1,4
       do 300 iii=1,1
       bisz=1.5
c      bisz=3.0
c      bisz=0.0
       ntap=10+iii
c      cp(11)=bisz*float(iii-2)
       cp(11)=bisz
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
       do 20 m=1,nb
c
!hpf$  independent,new(i,j,k)
CC MPI START
       do 22 k=ks,ke
CC MPI END
!hpf$  on home(u(:,:,k,:)),local(u,f,i,j,m) begin
       do 2211 j=1,ny2
       do 2211 i=1,nx2
       u(i,j,k,m)=0.0
       f(i,j,k,m)=0.0
 2211  continue
!hpf$  end on
   22  continue
   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(4+nb)=-1.0
       fbb(6+nb)=-1.0
c
c      initial condition
c      call equib7(f,pp,nxp,ro01,pr01,rra1,cj,cp)
       if(iii.eq.1)
     1 call equib8(f,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      do 410 ii=1,itap
c      do 174 m=1,nb
c      do 174 k=1,nz2
c      read(10) ((gf(i,j,k,m),i=1,nx2),j=1,ny2)
c 174  continue
c 410  continue
c 400  continue
c
c!hpf$  asynchronous(id1),nobuffer begin
c      f(1:nx2,1:ny2,1:nz2,1:nb)=gf(1:nx2,1:ny2,1:nz2,1:nb)
c!hpf$  end asynchronous
c!hpf$  asyncwait(id1)
c
c
       do 410 ii=1,itap
       do 174 m=1,nb
CC MPI START
       if(irank.eq.0) then
       do 1742 k=1,nz2
c      read(10) gf
c!hpf$  asynchronous(id1),nobuffer begin
c        f(1:nx2,1:ny2,k,m) = gf(1:nx2,1:ny2)
c!hpf$  end asynchronous
c!hpf$  asyncwait(id1)
cc       read(10) fg(1:nx2,1:ny2,k)
cc       read(10) f(1:nx2,1:ny2,k,m)
 1742  continue
       end if
cc         call mpi_scatterv(fg(1,1,1),recvcount,displs,
cc     *                    mpi_real,f(1,1,1,m),nword,mpi_real,0,
cc     *                    mpi_comm_world,ier)
CC MPI END
  174  continue
  410  continue
  400  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
       do 100 ii=1,last
c
       iiq=iiq+1
       iip=iip+1
       iis=iis+1
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
!hpf$  independent,new(i,j,k,m)
CC MPI START
       do 31 k=ks,ke
CC MPI END
!hpf$  on home(f(:,:,k,:)),local(f,i,j,m,x,xx,xx3,xx4,
!hpf$*         hhx) begin
       do 30 m=1,nb
CC MPI START
       if(k.eq.1 .and. irank.eq.0) then
CC MPI END
       do 311 j=2,ny1
       f(2,j,1,m)=f(1,j,2,m)
c      f(2,j,k,m)=f(1,j,k+1,m)
       do 311 i=3,nx1
       x=hhx(i)+xx4
       x=1.0-x/xx3
       x=amin1(x,1.0)
       x=1.0+x/3.0
       xx=1.0-x
       f(i,j,1,m)=x*f(i-1,j,2,m)+xx*f(i-2,j,3,m)
c      f(i,j,k,m)=x*f(i-1,j,k+1,m)+xx*f(i-2,j,k+2,m)
  311  continue
CC MPI START
       else if(k.eq.ke .and. irank.eq.isize-1) then
CC MPI END
       do 312 j=2,ny1
c      f(2,j,nz2,m)=f(1,j,nz1,m)
       f(2,j,k,m)=f(1,j,k-1,m)
       do 312 i=3,nx1
       x=hhx(i)+xx4
       x=1.0-x/xx3
       x=amin1(x,1.0)
       x=1.0+x/3.0
       xx=1.0-x
c      f(i,j,nz2,m)=x*f(i-1,j,nz1,m)+xx*f(i-2,j,nz,m)
       f(i,j,k,m)=x*f(i-1,j,k-1,m)+xx*f(i-2,j,k-2,m)
  312  continue
       end if
   30  continue
!hpf$  end on
   31  continue
c
c
c!hpf$  independent,new(k)
c      do 32 k=1,nz2
c      f(2:nx1,1,k,1:nb) = f(2:nx1,2,-k+nz3,1:nb)
c      end do
c
c!hpf$  asynchronous(id1),nobuffer begin
c      f(2:nx1,1,  1,m) = f(2:nx1,2,1+nz,m)
c      f(1:nx2,1:ny2,nz2,m) = f(1:nx2,1:ny2,   2,m)
c      f(2:nx1,1,1:nz2,1:nb) = f(2:nx1,2,nz2:1:-1,1:nb)
cc      f(2:nx1,1,nz2:1:-1,1:nb) = f(2:nx1,2,1:nz2,1:nb)
c      f(2:nx1,1,1:nz2,m) = f(2:nx1,2,nz2:1:-1,m)
c!hpf$  end asynchronous
c!hpf$  asyncwait(id1)
c
c      do 322 m=1,nb
CC MPI START
c        call mpi_gatherv(f(1,1,1,m),nword,mpi_real,fg(1,1,1),
c    *                    recvcount,displs,mpi_real,0,
c    *                    mpi_comm_world,ier)
c
c      if(irank.eq.0) then
c      do 2732 k=1,nz2
c      do 2732 i=2,nx1
c      fg(i,1,k)=fg(i,2,-k+nz3)
c2732  continue
c      end if
c
c        call mpi_scatterv(fg(1,1,1),recvcount,displs,
c    *                    mpi_real,f(1,1,1,m),nword,mpi_real,0,
c    *                    mpi_comm_world,ier)
CC MPI END
c 322  continue
c      call mpi_barrier(mpi_comm_world,ier)
c
c
CC MPI START
       do 1320 k=ks,ke
       do 1320 m=1,nb
       do 1320 i=1,nx2
       fq(i,m,k)=f(i,2,k,m)
 1320  continue
         call mpi_gatherv(fq(1,1,1),nwfg,mpi_real,fqg(1,1,1),
     *                    rcvfg,topfg,mpi_real,0,
     *                    mpi_comm_world,ier)
       call mpi_barrier(mpi_comm_world,ier)
       if(irank.eq.0) then
       do 1330 k=1,nz2
       do 1332 m=1,nb
       do 1332 i=1,nx2
       fqg1(i,m,k)=fqg(i,m,k)
 1332  continue
       do 1330 m=1,nb
       do 1330 i=1,nx2
       fqg(i,m,k)=fqg1(i,m,nz3-k)
 1330  continue
       end if
c      call mpi_barrier(mpi_comm_world,ier)
         call mpi_scatterv(fqg(1,1,1),rcvfg,topfg,
     *                    mpi_real,fq(1,1,1),nwfg,mpi_real,0,
     *                    mpi_comm_world,ier)
       call mpi_barrier(mpi_comm_world,ier)
       do 1340 k=ks,ke
       do 1340 m=1,nb
       do 1340 i=1,nx2
       f(i,1,k,m)=fq(i,m,k)
 1340  continue
CC MPI END
c
c
!hpf$  independent,new(i,j,k,m)
CC MPI START
       do 32 k=ks,ke
CC MPI END
!hpf$  on home(f(:,:,k,:)),local(f,u,fbb,i,j,m,x,xx,
!hpf$*         xx3,xx4,hhx,pr02,ro02) begin
       do 3311 m=1,nb
       f(2,ny2,k,m)=f(1,ny1,k,m)
c      f(2,1,k,m)=f(2,2,k,m)*fbb(m+nb)
       f(2,1,k,m)=f(2,1,k,m)*fbb(m+nb)
       do 3211 i=3,nx1
       x=hhx(i)+xx4
       x=1.0-x/xx3
       x=amin1(x,1.0)
       x=1.0+x/3.0
       xx=1.0-x
       f(i,ny2,k,m)=x*f(i-1,ny1,k,m)+xx*f(i-2,ny,k,m)
c      f(i,1,k,m)=f(i,2,k,m)*fbb(m+nb)
       f(i,1,k,m)=f(i,1,k,m)*fbb(m+nb)
 3211  continue
c
c      boundary condition at nx=nx2
       do 3311 j=1,ny2
       f(nx2,j,k,m)=f(nx1,j,k,m)
 3311  continue
c
c
       do 3411 j=1,ny2
       do 3411 i=1,nx2
       f(i,j,k,5)=amax1(f(i,j,k,5),pr02)
       f(i,j,k,1)=amax1(f(i,j,k,1),ro02)
       u(i,j,k,5)=amax1(u(i,j,k,5),pr02)
       u(i,j,k,1)=amax1(u(i,j,k,1),ro02)
 3411  continue
!hpf$  end on
   32  continue
c
c
!hpf$  reflect f
c
CC MPI START
       do m=1,nb
c        len=nx2*ny2
         if (irank.gt.0) then
           call mpi_send(f(1,1,ks,m),n2,mpi_real,irank-1,
     &                   100,mpi_comm_world,ier)
         end if
         if (irank.lt.isize-1) then
           call mpi_recv(f(1,1,ke+1,m),n2,mpi_real,irank+1,
     &                   100,mpi_comm_world,istatus,ier)
         end if
       end do
       call mpi_barrier(mpi_comm_world,ier)
CC MPI END
c
       if(ii.ne.1.and.iiq.ne.iiq0) go to 172
c      do 372 m=1,nb
!hpf$  independent,new(i,j,k,m)
CC MPI START
       do 371 k=ks,ke1
CC MPI END
!hpf$  on home(u(:,:,k,:)),local(f,u,i,j,m) begin
       do 3711 m=1,nb
       do 3711 j=1,ny1
       do 3711 i=1,nx1
       u(i,j,k,m)=0.125*(f(i,j,k,m)+f(i+1,j,k,m)
     1           +f(i,j+1,k,m)+f(i+1,j+1,k,m)
     2           +f(i,j,k+1,m)+f(i+1,j,k+1,m)
     3           +f(i,j+1,k+1,m)+f(i+1,j+1,k+1,m))
 3711  continue
!hpf$  end on
  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.ii.ne.1) go to 172
c      do 171 m=1,nb
c      do 171 k=1,nz2
c      write(ntap) ((gf(i,j,k,m),i=1,nx2),j=1,ny2)
c 171  continue
  172  continue
c
c
c
c      first step
c      step of k=k
!hpf$  independent,new(i,j,k)
CC MPI START
       do 90 k=ks,ke1
CC MPI END
!hpf$  on home(u(:,:,k,:)),local(f,u,p,pp,ff,i,j,m,
!hpf$*         x,y,z,hx,hy,hz,xx4,hhx,ar2,
!hpf$*         ro2,ro3,dx1,dy1,dz1,gam,eud,gra,ar1,
!hpf$*         uxd,x1,x2,x3,x5,y1,y5,ro02,pr02,vmax) begin
c
c      current
       do 38 j=1,ny1
       do 38 i=1,nx1
       p(i,j,k,11)=0.25*((f(i+1,j+1,k+1,7)+f(i+1,j,k+1,7)
     1   +f(i+1,j+1,k,7)+f(i+1,j,k,7)
     2   -f(i,j+1,k+1,7)-f(i,j,k+1,7)-f(i,j+1,k,7)-f(i,j,k,7))/hx
     3   -(f(i+1,j+1,k+1,6)-f(i+1,j,k+1,6)+f(i+1,j+1,k,6)-f(i+1,j,k,6)
     4   +f(i,j+1,k+1,6)-f(i,j,k+1,6)+f(i,j+1,k,6)-f(i,j,k,6))/hy)
       p(i,j,k,10)=0.25*((f(i+1,j+1,k+1,6)+f(i+1,j,k+1,6)
     1   -f(i+1,j+1,k,6)-f(i+1,j,k,6)
     2   +f(i,j+1,k+1,6)+f(i,j,k+1,6)-f(i,j+1,k,6)-f(i,j,k,6))/hz
     3   -(f(i+1,j+1,k+1,8)+f(i+1,j,k+1,8)+f(i+1,j+1,k,8)+f(i+1,j,k,8)
     4   -f(i,j+1,k+1,8)-f(i,j,k+1,8)-f(i,j+1,k,8)-f(i,j,k,8))/hx)
       p(i,j,k,9)=0.25*((f(i+1,j+1,k+1,8)-f(i+1,j,k+1,8)
     1   +f(i+1,j+1,k,8)-f(i+1,j,k,8)
     2   +f(i,j+1,k+1,8)-f(i,j,k+1,8)+f(i,j+1,k,8)-f(i,j,k,8))/hy
     3   -(f(i+1,j+1,k+1,7)+f(i+1,j,k+1,7)-f(i+1,j+1,k,7)-f(i+1,j,k,7)
     4   +f(i,j+1,k+1,7)+f(i,j,k+1,7)-f(i,j+1,k,7)-f(i,j,k,7))/hz)
   38  continue
       do 382 j=1,ny1
       do 382 i=1,nx1
       p(i,j,k,11)=p(i,j,k,11)-pp(i,j,k,3)
       p(i,j,k,10)=p(i,j,k,10)-pp(i,j,k,2)
       p(i,j,k,9)=p(i,j,k,9)-pp(i,j,k,1)
  382  continue
c
       do 381 m=1,nb
       do 381 j=1,ny1
       do 381 i=1,nx1
       p(i,j,k,m)=0.125*(f(i,j,k,m)+f(i+1,j,k,m)
     1      +f(i,j+1,k,m)+f(i+1,j+1,k,m)
     2      +f(i,j,k+1,m)+f(i+1,j,k+1,m)
     3      +f(i,j+1,k+1,m)+f(i+1,j+1,k+1,m))
c      u(i,j,k,m)=p(i,j,k,m)
  381  continue
c
c      first step
       xx4=0.5*hx*float(2*nxp-nx1-1)
       z=0.5*hz*float(2*(kss+k)-nz2)
       do 50 j=1,ny1
       y=0.5*hy*float(2*j-2)
       do 502 m=1,nb
       do 502 i=1,nx1
       ff(i,j,k,m)=u(i,j,k,m)
 502  continue
c
       do 501 i=1,nx1
       x=hhx(i)+xx4
       ro2=x*x+y*y+z*z
       ro2=amax1(ro2,ar2)
       ro3=ro2*sqrt(ro2)
c
       u(i,j,k,8)=u(i,j,k,8)
     1   +dx1*(f(i+1,j+1,k+1,4)*f(i+1,j+1,k+1,6)
     1   -f(i+1,j+1,k+1,2)*f(i+1,j+1,k+1,8)
     2   +f(i+1,j,k+1,4)*f(i+1,j,k+1,6)-f(i+1,j,k+1,2)*f(i+1,j,k+1,8)
     3   +f(i+1,j+1,k,4)*f(i+1,j+1,k,6)-f(i+1,j+1,k,2)*f(i+1,j+1,k,8)
     4   +f(i+1,j,k,4)*f(i+1,j,k,6)-f(i+1,j,k,2)*f(i+1,j,k,8)
     5   -f(i,j+1,k+1,4)*f(i,j+1,k+1,6)+f(i,j+1,k+1,2)*f(i,j+1,k+1,8)
     6   -f(i,j,k+1,4)*f(i,j,k+1,6)+f(i,j,k+1,2)*f(i,j,k+1,8)
     7   -f(i,j+1,k,4)*f(i,j+1,k,6)+f(i,j+1,k,2)*f(i,j+1,k,8)
     8   -f(i,j,k,4)*f(i,j,k,6)+f(i,j,k,2)*f(i,j,k,8))
     1   -dy1*(f(i+1,j+1,k+1,3)*f(i+1,j+1,k+1,8)
     1   -f(i+1,j+1,k+1,4)*f(i+1,j+1,k+1,7)
     2   -f(i+1,j,k+1,3)*f(i+1,j,k+1,8)+f(i+1,j,k+1,4)*f(i+1,j,k+1,7)
     3   +f(i+1,j+1,k,3)*f(i+1,j+1,k,8)-f(i+1,j+1,k,4)*f(i+1,j+1,k,7)
     4   -f(i+1,j,k,3)*f(i+1,j,k,8)+f(i+1,j,k,4)*f(i+1,j,k,7)
     5   +f(i,j+1,k+1,3)*f(i,j+1,k+1,8)-f(i,j+1,k+1,4)*f(i,j+1,k+1,7)
     6   -f(i,j,k+1,3)*f(i,j,k+1,8)+f(i,j,k+1,4)*f(i,j,k+1,7)
     7   +f(i,j+1,k,3)*f(i,j+1,k,8)-f(i,j+1,k,4)*f(i,j+1,k,7)
     8   -f(i,j,k,3)*f(i,j,k,8)+f(i,j,k,4)*f(i,j,k,7))
c
       u(i,j,k,7)=u(i,j,k,7)
     1   +dz1*(f(i+1,j+1,k+1,3)*f(i+1,j+1,k+1,8)
     1   -f(i+1,j+1,k+1,4)*f(i+1,j+1,k+1,7)
     2   +f(i+1,j,k+1,3)*f(i+1,j,k+1,8)-f(i+1,j,k+1,4)*f(i+1,j,k+1,7)
     3   -f(i+1,j+1,k,3)*f(i+1,j+1,k,8)+f(i+1,j+1,k,4)*f(i+1,j+1,k,7)
     4   -f(i+1,j,k,3)*f(i+1,j,k,8)+f(i+1,j,k,4)*f(i+1,j,k,7)
     5   +f(i,j+1,k+1,3)*f(i,j+1,k+1,8)-f(i,j+1,k+1,4)*f(i,j+1,k+1,7)
     6   +f(i,j,k+1,3)*f(i,j,k+1,8)-f(i,j,k+1,4)*f(i,j,k+1,7)
     7   -f(i,j+1,k,3)*f(i,j+1,k,8)+f(i,j+1,k,4)*f(i,j+1,k,7)
     8   -f(i,j,k,3)*f(i,j,k,8)+f(i,j,k,4)*f(i,j,k,7))
     1   -dx1*(f(i+1,j+1,k+1,2)*f(i+1,j+1,k+1,7)
     1   -f(i+1,j+1,k+1,3)*f(i+1,j+1,k+1,6)
     2   +f(i+1,j,k+1,2)*f(i+1,j,k+1,7)-f(i+1,j,k+1,3)*f(i+1,j,k+1,6)
     3   +f(i+1,j+1,k,2)*f(i+1,j+1,k,7)-f(i+1,j+1,k,3)*f(i+1,j+1,k,6)
     4   +f(i+1,j,k,2)*f(i+1,j,k,7)-f(i+1,j,k,3)*f(i+1,j,k,6)
     5   -f(i,j+1,k+1,2)*f(i,j+1,k+1,7)+f(i,j+1,k+1,3)*f(i,j+1,k+1,6)
     6   -f(i,j,k+1,2)*f(i,j,k+1,7)+f(i,j,k+1,3)*f(i,j,k+1,6)
     7   -f(i,j+1,k,2)*f(i,j+1,k,7)+f(i,j+1,k,3)*f(i,j+1,k,6)
     8   -f(i,j,k,2)*f(i,j,k,7)+f(i,j,k,3)*f(i,j,k,6))
c
       u(i,j,k,6)=u(i,j,k,6)
     1   +dy1*(f(i+1,j+1,k+1,2)*f(i+1,j+1,k+1,7)
     1   -f(i+1,j+1,k+1,3)*f(i+1,j+1,k+1,6)
     2   -f(i+1,j,k+1,2)*f(i+1,j,k+1,7)+f(i+1,j,k+1,3)*f(i+1,j,k+1,6)
     3   +f(i+1,j+1,k,2)*f(i+1,j+1,k,7)-f(i+1,j+1,k,3)*f(i+1,j+1,k,6)
     4   -f(i+1,j,k,2)*f(i+1,j,k,7)+f(i+1,j,k,3)*f(i+1,j,k,6)
     5   +f(i,j+1,k+1,2)*f(i,j+1,k+1,7)-f(i,j+1,k+1,3)*f(i,j+1,k+1,6)
     6   -f(i,j,k+1,2)*f(i,j,k+1,7)+f(i,j,k+1,3)*f(i,j,k+1,6)
     7   +f(i,j+1,k,2)*f(i,j+1,k,7)-f(i,j+1,k,3)*f(i,j+1,k,6)
     8   -f(i,j,k,2)*f(i,j,k,7)+f(i,j,k,3)*f(i,j,k,6))
     1   -dz1*(f(i+1,j+1,k+1,4)*f(i+1,j+1,k+1,6)
     1   -f(i+1,j+1,k+1,2)*f(i+1,j+1,k+1,8)
     2   +f(i+1,j,k+1,4)*f(i+1,j,k+1,6)-f(i+1,j,k+1,2)*f(i+1,j,k+1,8)
     3   -f(i+1,j+1,k,4)*f(i+1,j+1,k,6)+f(i+1,j+1,k,2)*f(i+1,j+1,k,8)
     4   -f(i+1,j,k,4)*f(i+1,j,k,6)+f(i+1,j,k,2)*f(i+1,j,k,8)
     5   +f(i,j+1,k+1,4)*f(i,j+1,k+1,6)-f(i,j+1,k+1,2)*f(i,j+1,k+1,8)
     6   +f(i,j,k+1,4)*f(i,j,k+1,6)-f(i,j,k+1,2)*f(i,j,k+1,8)
     7   -f(i,j+1,k,4)*f(i,j+1,k,6)+f(i,j+1,k,2)*f(i,j+1,k,8)
     8   -f(i,j,k,4)*f(i,j,k,6)+f(i,j,k,2)*f(i,j,k,8))
c
       x3=dx1*(f(i+1,j+1,k+1,2)+f(i+1,j,k+1,2)
     1   +f(i+1,j+1,k,2)+f(i+1,j,k,2)
     1   -f(i,j+1,k+1,2)-f(i,j,k+1,2)-f(i,j+1,k,2)-f(i,j,k,2))
     2   +dy1*(f(i+1,j+1,k+1,3)-f(i+1,j,k+1,3)
     2   +f(i+1,j+1,k,3)-f(i+1,j,k,3)
     3   +f(i,j+1,k+1,3)-f(i,j,k+1,3)+f(i,j+1,k,3)-f(i,j,k,3))
     4   +dz1*(f(i+1,j+1,k+1,4)+f(i+1,j,k+1,4)
     4   -f(i+1,j+1,k,4)-f(i+1,j,k,4)
     5   +f(i,j+1,k+1,4)+f(i,j,k+1,4)-f(i,j+1,k,4)-f(i,j,k,4))
c
       u(i,j,k,5)=u(i,j,k,5)-p(i,j,k,5)*gam*x3
     1   -dx1*p(i,j,k,2)*(f(i+1,j+1,k+1,5)+f(i+1,j,k+1,5)
     1   +f(i+1,j+1,k,5)+f(i+1,j,k,5)
     2   -f(i,j+1,k+1,5)-f(i,j,k+1,5)-f(i,j+1,k,5)-f(i,j,k,5))
     3   -dy1*p(i,j,k,3)*(f(i+1,j+1,k+1,5)-f(i+1,j,k+1,5)
     3   +f(i+1,j+1,k,5)-f(i+1,j,k,5)
     4   +f(i,j+1,k+1,5)-f(i,j,k+1,5)+f(i,j+1,k,5)-f(i,j,k,5))
     5   -dz1*p(i,j,k,4)*(f(i+1,j+1,k+1,5)+f(i+1,j,k+1,5)
     5   -f(i+1,j+1,k,5)-f(i+1,j,k,5)
     6   +f(i,j+1,k+1,5)+f(i,j,k+1,5)-f(i,j+1,k,5)-f(i,j,k,5))
c
c
       u(i,j,k,4)=u(i,j,k,4)+t1*((p(i,j,k,9)*p(i,j,k,7)
     1   -p(i,j,k,10)*p(i,j,k,6))/p(i,j,k,1)
     1   -0.0*eud*p(i,j,k,4)-gra*z/ro3)
     2   -dx1*p(i,j,k,2)*(f(i+1,j+1,k+1,4)+f(i+1,j,k+1,4)
     2   +f(i+1,j+1,k,4)+f(i+1,j,k,4)
     3   -f(i,j+1,k+1,4)-f(i,j,k+1,4)-f(i,j+1,k,4)-f(i,j,k,4))
     4   -dy1*p(i,j,k,3)*(f(i+1,j+1,k+1,4)-f(i+1,j,k+1,4)
     4   +f(i+1,j+1,k,4)-f(i+1,j,k,4)
     5   +f(i,j+1,k+1,4)-f(i,j,k+1,4)+f(i,j+1,k,4)-f(i,j,k,4))
     6   -dz1*p(i,j,k,4)*(f(i+1,j+1,k+1,4)+f(i+1,j,k+1,4)
     6   -f(i+1,j+1,k,4)-f(i+1,j,k,4)
     7   +f(i,j+1,k+1,4)+f(i,j,k+1,4)-f(i,j+1,k,4)-f(i,j,k,4))
     8   -dz1*(f(i+1,j+1,k+1,5)+f(i+1,j,k+1,5)
     8   -f(i+1,j+1,k,5)-f(i+1,j,k,5)
     9   +f(i,j+1,k+1,5)+f(i,j,k+1,5)-f(i,j+1,k,5)-f(i,j,k,5))
     9   /p(i,j,k,1)
c
       u(i,j,k,3)=u(i,j,k,3)+t1*((p(i,j,k,11)*p(i,j,k,6)
     1   -p(i,j,k,9)*p(i,j,k,8))/p(i,j,k,1)
     1   -0.0*eud*p(i,j,k,3)-gra*y/ro3)
     2   -dx1*p(i,j,k,2)*(f(i+1,j+1,k+1,3)+f(i+1,j,k+1,3)
     2   +f(i+1,j+1,k,3)+f(i+1,j,k,3)
     3   -f(i,j+1,k+1,3)-f(i,j,k+1,3)-f(i,j+1,k,3)-f(i,j,k,3))
     4   -dy1*p(i,j,k,3)*(f(i+1,j+1,k+1,3)-f(i+1,j,k+1,3)
     4   +f(i+1,j+1,k,3)-f(i+1,j,k,3)
     5   +f(i,j+1,k+1,3)-f(i,j,k+1,3)+f(i,j+1,k,3)-f(i,j,k,3))
     6   -dz1*p(i,j,k,4)*(f(i+1,j+1,k+1,3)+f(i+1,j,k+1,3)
     6   -f(i+1,j+1,k,3)-f(i+1,j,k,3)
     7   +f(i,j+1,k+1,3)+f(i,j,k+1,3)-f(i,j+1,k,3)-f(i,j,k,3))
     8   -dy1*(f(i+1,j+1,k+1,5)-f(i+1,j,k+1,5)
     8   +f(i+1,j+1,k,5)-f(i+1,j,k,5)
     9   +f(i,j+1,k+1,5)-f(i,j,k+1,5)+f(i,j+1,k,5)-f(i,j,k,5))
     9   /p(i,j,k,1)
c
       uxd=amax1((x-4.0*ar1)/(4.0*ar1),0.0)
       uxd=amin1(uxd,1.0)
       uxd=amin1(uxd*p(i,j,k,2),0.0)
c
       u(i,j,k,2)=u(i,j,k,2)+t1*((p(i,j,k,10)*p(i,j,k,8)
     1   -p(i,j,k,11)*p(i,j,k,7))/p(i,j,k,1)
     1   -eud*uxd-gra*x/ro3)
     2   -dx1*p(i,j,k,2)*(f(i+1,j+1,k+1,2)+f(i+1,j,k+1,2)
     2   +f(i+1,j+1,k,2)+f(i+1,j,k,2)
     3   -f(i,j+1,k+1,2)-f(i,j,k+1,2)-f(i,j+1,k,2)-f(i,j,k,2))
     4   -dy1*p(i,j,k,3)*(f(i+1,j+1,k+1,2)-f(i+1,j,k+1,2)
     4   +f(i+1,j+1,k,2)-f(i+1,j,k,2)
     5   +f(i,j+1,k+1,2)-f(i,j,k+1,2)+f(i,j+1,k,2)-f(i,j,k,2))
     6   -dz1*p(i,j,k,4)*(f(i+1,j+1,k+1,2)+f(i+1,j,k+1,2)
     6   -f(i+1,j+1,k,2)-f(i+1,j,k,2)
     7   +f(i,j+1,k+1,2)+f(i,j,k+1,2)-f(i,j+1,k,2)-f(i,j,k,2))
     8   -dx1*(f(i+1,j+1,k+1,5)+f(i+1,j,k+1,5)
     8   +f(i+1,j+1,k,5)+f(i+1,j,k,5)
     9   -f(i,j+1,k+1,5)-f(i,j,k+1,5)-f(i,j+1,k,5)-f(i,j,k,5))
     9   /p(i,j,k,1)
c
       u(i,j,k,1)=u(i,j,k,1)
     1   -dx1*(f(i+1,j+1,k+1,2)*f(i+1,j+1,k+1,1)
     1   +f(i+1,j,k+1,2)*f(i+1,j,k+1,1)
     2   +f(i+1,j+1,k,2)*f(i+1,j+1,k,1)+f(i+1,j,k,2)*f(i+1,j,k,1)
     3   -f(i,j+1,k+1,2)*f(i,j+1,k+1,1)-f(i,j,k+1,2)*f(i,j,k+1,1)
     4   -f(i,j+1,k,2)*f(i,j+1,k,1)-f(i,j,k,2)*f(i,j,k,1))
     5   -dy1*(f(i+1,j+1,k+1,3)*f(i+1,j+1,k+1,1)
     5   -f(i+1,j,k+1,3)*f(i+1,j,k+1,1)
     6   +f(i+1,j+1,k,3)*f(i+1,j+1,k,1)-f(i+1,j,k,3)*f(i+1,j,k,1)
     7   +f(i,j+1,k+1,3)*f(i,j+1,k+1,1)-f(i,j,k+1,3)*f(i,j,k+1,1)
     8   +f(i,j+1,k,3)*f(i,j+1,k,1)-f(i,j,k,3)*f(i,j,k,1))
     1   -dz1*(f(i+1,j+1,k+1,4)*f(i+1,j+1,k+1,1)
     1   +f(i+1,j,k+1,4)*f(i+1,j,k+1,1)
     2   -f(i+1,j+1,k,4)*f(i+1,j+1,k,1)-f(i+1,j,k,4)*f(i+1,j,k,1)
     3   +f(i,j+1,k+1,4)*f(i,j+1,k+1,1)+f(i,j,k+1,4)*f(i,j,k+1,1)
     4   -f(i,j+1,k,4)*f(i,j+1,k,1)-f(i,j,k,4)*f(i,j,k,1))
  501  continue
c
       do 520 i=1,nx1
       x=hhx(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(i,j,k,8)=u(i,j,k,8)*x2+ff(i,j,k,8)*(1.0-x2)
       u(i,j,k,7)=u(i,j,k,7)*x2+ff(i,j,k,7)*(1.0-x2)
       u(i,j,k,6)=u(i,j,k,6)*x2+ff(i,j,k,6)*(1.0-x2)
       u(i,j,k,5)=u(i,j,k,5)*x2+ff(i,j,k,5)*(1.0-x2)
       u(i,j,k,4)=u(i,j,k,4)*x2+ff(i,j,k,4)*(1.0-x2)
       u(i,j,k,3)=u(i,j,k,3)*x2+ff(i,j,k,3)*(1.0-x2)
       u(i,j,k,2)=u(i,j,k,2)*x2+ff(i,j,k,2)*(1.0-x2)
       u(i,j,k,1)=u(i,j,k,1)*x2+ff(i,j,k,1)*(1.0-x2)
  520  continue
c
       do 524 i=1,nx1
       x5=u(i,j,k,5)
       x1=u(i,j,k,1)
       y5=amax1(x5,pr02)
       y1=amax1(x1,ro02)
c      x2=u(i,j,k,2)*u(i,j,k,2)+u(i,j,k,3)*u(i,j,k,3)+u(i,j,k,4)*u(i,j,k,4)
c      vmax=amax1(vmax,x2)
       x=amax1(x1,0.0)/y1
       u(i,j,k,5)=y5
       u(i,j,k,4)=x*u(i,j,k,4)
       u(i,j,k,3)=x*u(i,j,k,3)
       u(i,j,k,2)=x*u(i,j,k,2)
       u(i,j,k,1)=y1
  524  continue
   50  continue
!hpf$  end on
   90  continue
c
c
c
c      do 70 m=1,nb
!hpf$  independent,new(i,j,m)
CC MPI START
       do 72 k=ks,ke
CC MPI END
!hpf$  on home(ff(:,:,k,:)),local(f,ff,i,j,m) begin
       do 7211 m=1,nb
       do 7211 j=1,ny2
       do 7211 i=1,nx2
       ff(i,j,k,m)=f(i,j,k,m)
 7211  continue
!hpf$  end on
   72  continue
c  70  continue
c
!hpf$  reflect u
!hpf$  reflect ff
!hpf$  reflect pp
c
CC MPI START
       do m=1,nb
c        len=nx2*ny2
         if (irank.lt.isize-1) then
           call mpi_send(u(1,1,ke,m),n2,mpi_real,irank+1,
     &                   110,mpi_comm_world,ier)
           call mpi_send(ff(1,1,ke,m),n2,mpi_real,irank+1,
     &                   120,mpi_comm_world,ier)
         end if
         if (irank.gt.0) then
           call mpi_recv(u(1,1,ks-1,m),n2,mpi_real,irank-1,
     &                   110,mpi_comm_world,istatus,ier)
           call mpi_recv(ff(1,1,ks-1,m),n2,mpi_real,irank-1,
     &                   120,mpi_comm_world,istatus,ier)
         end if
         if (irank.gt.0) then
           call mpi_send(ff(1,1,ks,m),n2,mpi_real,irank-1,
     &                   125,mpi_comm_world,ier)
         end if
         if (irank.lt.isize-1) then
           call mpi_recv(ff(1,1,ke+1,m),n2,mpi_real,irank+1,
     &                   125,mpi_comm_world,istatus,ier)
         end if
       end do
       do m=1,3
c        len=nx2*ny2
         if (irank.lt.isize-1) then
           call mpi_send(pp(1,1,ke,m),n2,mpi_real,irank+1,
     &                   130,mpi_comm_world,ier)
         end if
         if (irank.gt.0) then
           call mpi_recv(pp(1,1,ks-1,m),n2,mpi_real,irank-1,
     &                   130,mpi_comm_world,istatus,ier)
         end if
       end do
       call mpi_barrier(mpi_comm_world,ier)
CC MPI END
c
c      second step
c
c
!hpf$  independent,new(i,j,k),reduction(max:vmax)
CC MPI START
       do 190 k=ks1,ke1
CC MPI END
!hpf$  on home(f(:,:,k,:)),local(f,u,p,pp,ff,i,j,m,
!hpf$*    x,y,z,hx,hy,hz,gam,pmu,
!hpf$*    xx4,ar1,ar2,ro2,ro3,hhx,x1,x2,x3,x5,eat,eat0,po0,
!hpf$*    dx2,dy2,dz2,dx4,dy4,dz4,eud,gra,rmu,uxd,
!hpf$*    ro02,pr02,y1,y5,vmax,aru) begin
c
c
c      current
       do 68 j=2,ny1
       do 68 i=2,nx1
       p(i,j,k,11)=0.25*((u(i,j,k,7)+u(i,j-1,k,7)
     1        +u(i,j,k-1,7)+u(i,j-1,k-1,7)
     2        -u(i-1,j,k,7)-u(i-1,j-1,k,7)
     3        -u(i-1,j,k-1,7)-u(i-1,j-1,k-1,7))/hx
     4        -(u(i,j,k,6)-u(i,j-1,k,6)
     5        +u(i,j,k-1,6)-u(i,j-1,k-1,6)
     6        +u(i-1,j,k,6)-u(i-1,j-1,k,6)
     7        +u(i-1,j,k-1,6)-u(i-1,j-1,k-1,6))/hy)
       p(i,j,k,10)=0.25*((u(i,j,k,6)+u(i,j-1,k,6)
     1        -u(i,j,k-1,6)-u(i,j-1,k-1,6)
     2        +u(i-1,j,k,6)+u(i-1,j-1,k,6)
     3        -u(i-1,j,k-1,6)-u(i-1,j-1,k-1,6))/hz
     4        -(u(i,j,k,8)+u(i,j-1,k,8)
     5        +u(i,j,k-1,8)+u(i,j-1,k-1,8)
     6        -u(i-1,j,k,8)-u(i-1,j-1,k,8)
     7        -u(i-1,j,k-1,8)-u(i-1,j-1,k-1,8))/hx)
       p(i,j,k,9)=0.25*((u(i,j,k,8)-u(i,j-1,k,8)
     1       +u(i,j,k-1,8)-u(i,j-1,k-1,8)
     2       +u(i-1,j,k,8)-u(i-1,j-1,k,8)
     3       +u(i-1,j,k-1,8)-u(i-1,j-1,k-1,8))/hy
     4       -(u(i,j,k,7)+u(i,j-1,k,7)
     5       -u(i,j,k-1,7)-u(i,j-1,k-1,7)
     6       +u(i-1,j,k,7)+u(i-1,j-1,k,7)
     7       -u(i-1,j,k-1,7)-u(i-1,j-1,k-1,7))/hz)
   68  continue
       do 682 j=2,ny1
       do 682 i=2,nx1
       p(i,j,k,11)=p(i,j,k,11)-0.125*(
     1          pp(i,j,k,3)+pp(i-1,j,k,3)
     2         +pp(i,j-1,k,3)+pp(i-1,j-1,k,3)
     3         +pp(i,j,k-1,3)+pp(i-1,j,k-1,3)
     4         +pp(i,j-1,k-1,3)+pp(i-1,j-1,k-1,3))
       p(i,j,k,10)=p(i,j,k,10)-0.125*(
     1          pp(i,j,k,2)+pp(i-1,j,k,2)
     2         +pp(i,j-1,k,2)+pp(i-1,j-1,k,2)
     3         +pp(i,j,k-1,2)+pp(i-1,j,k-1,2)
     4         +pp(i,j-1,k-1,2)+pp(i-1,j-1,k-1,2))
       p(i,j,k,9)=p(i,j,k,9)-0.125*(
     1          pp(i,j,k,1)+pp(i-1,j,k,1)
     2         +pp(i,j-1,k,1)+pp(i-1,j-1,k,1)
     3         +pp(i,j,k-1,1)+pp(i-1,j,k-1,1)
     4         +pp(i,j-1,k-1,1)+pp(i-1,j-1,k-1,1))
  682  continue
c
       do 681 m=1,nb
       do 681 j=2,ny1
       do 681 i=2,nx1
       p(i,j,k,m)=0.125*(u(i,j,k,m)+u(i-1,j,k,m)
     1      +u(i,j-1,k,m)+u(i-1,j-1,k,m)
     2      +u(i,j,k-1,m)+u(i-1,j,k-1,m)
     3      +u(i,j-1,k-1,m)+u(i-1,j-1,k-1,m))
  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)
       do 80 j=2,ny1
       y=0.5*hy*float(2*j-3)
c
       do 801 i=2,nx1
       x=hhx(i)+xx4
       ro2=x*x+y*y+z*z
       ro2=amax1(ro2,ar2)
       ro3=ro2*sqrt(ro2)
c
       x1=ff(i,j,k,5)/(ff(i,j,k,1)*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(i+1,j,k,8)-2.0*ff(i,j,k,8)+ff(i-1,j,k,8))
     1        +dy4*(ff(i,j+1,k,8)-2.0*ff(i,j,k,8)+ff(i,j-1,k,8))
     2        +dz4*(ff(i,j,k+1,8)-2.0*ff(i,j,k,8)+ff(i,j,k-1,8)))
       f(i,j,k,8)=ff(i,j,k,8)+x1
     1   +dx2*(u(i,j,k,4)*u(i,j,k,6)-u(i,j,k,2)*u(i,j,k,8)
     2   +u(i,j-1,k,4)*u(i,j-1,k,6)-u(i,j-1,k,2)*u(i,j-1,k,8)
     3   +u(i,j,k-1,4)*u(i,j,k-1,6)-u(i,j,k-1,2)*u(i,j,k-1,8)
     4   +u(i,j-1,k-1,4)*u(i,j-1,k-1,6)-u(i,j-1,k-1,2)*u(i,j-1,k-1,8)
     5   -u(i-1,j,k,4)*u(i-1,j,k,6)+u(i-1,j,k,2)*u(i-1,j,k,8)
     6   -u(i-1,j-1,k,4)*u(i-1,j-1,k,6)+u(i-1,j-1,k,2)*u(i-1,j-1,k,8)
     7   -u(i-1,j,k-1,4)*u(i-1,j,k-1,6)+u(i-1,j,k-1,2)*u(i-1,j,k-1,8)
     8   -u(i-1,j-1,k-1,4)*u(i-1,j-1,k-1,6)
     8   +u(i-1,j-1,k-1,2)*u(i-1,j-1,k-1,8))
     1   -dy2*(u(i,j,k,3)*u(i,j,k,8)-u(i,j,k,4)*u(i,j,k,7)
     2   -u(i,j-1,k,3)*u(i,j-1,k,8)+u(i,j-1,k,4)*u(i,j-1,k,7)
     3   +u(i,j,k-1,3)*u(i,j,k-1,8)-u(i,j,k-1,4)*u(i,j,k-1,7)
     4   -u(i,j-1,k-1,3)*u(i,j-1,k-1,8)+u(i,j-1,k-1,4)*u(i,j-1,k-1,7)
     5   +u(i-1,j,k,3)*u(i-1,j,k,8)-u(i-1,j,k,4)*u(i-1,j,k,7)
     6   -u(i-1,j-1,k,3)*u(i-1,j-1,k,8)+u(i-1,j-1,k,4)*u(i-1,j-1,k,7)
     7   +u(i-1,j,k-1,3)*u(i-1,j,k-1,8)-u(i-1,j,k-1,4)*u(i-1,j,k-1,7)
     8   -u(i-1,j-1,k-1,3)*u(i-1,j-1,k-1,8)
     8   +u(i-1,j-1,k-1,4)*u(i-1,j-1,k-1,7))
c
       x1=eat*(dx4*(ff(i+1,j,k,7)-2.0*ff(i,j,k,7)+ff(i-1,j,k,7))
     1        +dy4*(ff(i,j+1,k,7)-2.0*ff(i,j,k,7)+ff(i,j-1,k,7))
     2        +dz4*(ff(i,j,k+1,7)-2.0*ff(i,j,k,7)+ff(i,j,k-1,7)))
       f(i,j,k,7)=ff(i,j,k,7)+x1
     1   +dz2*(u(i,j,k,3)*u(i,j,k,8)-u(i,j,k,4)*u(i,j,k,7)
     2   +u(i,j-1,k,3)*u(i,j-1,k,8)-u(i,j-1,k,4)*u(i,j-1,k,7)
     3   -u(i,j,k-1,3)*u(i,j,k-1,8)+u(i,j,k-1,4)*u(i,j,k-1,7)
     4   -u(i,j-1,k-1,3)*u(i,j-1,k-1,8)+u(i,j-1,k-1,4)*u(i,j-1,k-1,7)
     5   +u(i-1,j,k,3)*u(i-1,j,k,8)-u(i-1,j,k,4)*u(i-1,j,k,7)
     6   +u(i-1,j-1,k,3)*u(i-1,j-1,k,8)-u(i-1,j-1,k,4)*u(i-1,j-1,k,7)
     7   -u(i-1,j,k-1,3)*u(i-1,j,k-1,8)+u(i-1,j,k-1,4)*u(i-1,j,k-1,7)
     8   -u(i-1,j-1,k-1,3)*u(i-1,j-1,k-1,8)
     8   +u(i-1,j-1,k-1,4)*u(i-1,j-1,k-1,7))
     1   -dx2*(u(i,j,k,2)*u(i,j,k,7)-u(i,j,k,3)*u(i,j,k,6)
     2   +u(i,j-1,k,2)*u(i,j-1,k,7)-u(i,j-1,k,3)*u(i,j-1,k,6)
     3   +u(i,j,k-1,2)*u(i,j,k-1,7)-u(i,j,k-1,3)*u(i,j,k-1,6)
     4   +u(i,j-1,k-1,2)*u(i,j-1,k-1,7)-u(i,j-1,k-1,3)*u(i,j-1,k-1,6)
     5   -u(i-1,j,k,2)*u(i-1,j,k,7)+u(i-1,j,k,3)*u(i-1,j,k,6)
     6   -u(i-1,j-1,k,2)*u(i-1,j-1,k,7)+u(i-1,j-1,k,3)*u(i-1,j-1,k,6)
     7   -u(i-1,j,k-1,2)*u(i-1,j,k-1,7)+u(i-1,j,k-1,3)*u(i-1,j,k-1,6)
     8   -u(i-1,j-1,k-1,2)*u(i-1,j-1,k-1,7)
     8   +u(i-1,j-1,k-1,3)*u(i-1,j-1,k-1,6))
c
       x1=eat*(dx4*(ff(i+1,j,k,6)-2.0*ff(i,j,k,6)+ff(i-1,j,k,6))
     1        +dy4*(ff(i,j+1,k,6)-2.0*ff(i,j,k,6)+ff(i,j-1,k,6))
     2        +dz4*(ff(i,j,k+1,6)-2.0*ff(i,j,k,6)+ff(i,j,k-1,6)))
       f(i,j,k,6)=ff(i,j,k,6)+x1
     1   +dy2*(u(i,j,k,2)*u(i,j,k,7)-u(i,j,k,3)*u(i,j,k,6)
     2   -u(i,j-1,k,2)*u(i,j-1,k,7)+u(i,j-1,k,3)*u(i,j-1,k,6)
     3   +u(i,j,k-1,2)*u(i,j,k-1,7)-u(i,j,k-1,3)*u(i,j,k-1,6)
     4   -u(i,j-1,k-1,2)*u(i,j-1,k-1,7)+u(i,j-1,k-1,3)*u(i,j-1,k-1,6)
     5   +u(i-1,j,k,2)*u(i-1,j,k,7)-u(i-1,j,k,3)*u(i-1,j,k,6)
     6   -u(i-1,j-1,k,2)*u(i-1,j-1,k,7)+u(i-1,j-1,k,3)*u(i-1,j-1,k,6)
     7   +u(i-1,j,k-1,2)*u(i-1,j,k-1,7)-u(i-1,j,k-1,3)*u(i-1,j,k-1,6)
     8   -u(i-1,j-1,k-1,2)*u(i-1,j-1,k-1,7)
     8   +u(i-1,j-1,k-1,3)*u(i-1,j-1,k-1,6))
     1   -dz2*(u(i,j,k,4)*u(i,j,k,6)-u(i,j,k,2)*u(i,j,k,8)
     2   +u(i,j-1,k,4)*u(i,j-1,k,6)-u(i,j-1,k,2)*u(i,j-1,k,8)
     3   -u(i,j,k-1,4)*u(i,j,k-1,6)+u(i,j,k-1,2)*u(i,j,k-1,8)
     4   -u(i,j-1,k-1,4)*u(i,j-1,k-1,6)+u(i,j-1,k-1,2)*u(i,j-1,k-1,8)
     5   +u(i-1,j,k,4)*u(i-1,j,k,6)-u(i-1,j,k,2)*u(i-1,j,k,8)
     6   +u(i-1,j-1,k,4)*u(i-1,j-1,k,6)-u(i-1,j-1,k,2)*u(i-1,j-1,k,8)
     7   -u(i-1,j,k-1,4)*u(i-1,j,k-1,6)+u(i-1,j,k-1,2)*u(i-1,j,k-1,8)
     8   -u(i-1,j-1,k-1,4)*u(i-1,j-1,k-1,6)
     8   +u(i-1,j-1,k-1,2)*u(i-1,j-1,k-1,8))
c
       x3=dx2*(u(i,j,k,2)+u(i,j-1,k,2)+u(i,j,k-1,2)+u(i,j-1,k-1,2)
     1   -u(i-1,j,k,2)-u(i-1,j-1,k,2)-u(i-1,j,k-1,2)-u(i-1,j-1,k-1,2))
     2   +dy2*(u(i,j,k,3)-u(i,j-1,k,3)+u(i,j,k-1,3)-u(i,j-1,k-1,3)
     3   +u(i-1,j,k,3)-u(i-1,j-1,k,3)+u(i-1,j,k-1,3)-u(i-1,j-1,k-1,3))
     4   +dz2*(u(i,j,k,4)+u(i,j-1,k,4)-u(i,j,k-1,4)-u(i,j-1,k-1,4)
     5   +u(i-1,j,k,4)+u(i-1,j-1,k,4)-u(i-1,j,k-1,4)-u(i-1,j-1,k-1,4))
       f(i,j,k,5)=ff(i,j,k,5)-p(i,j,k,5)*gam*x3
     1   -dx2*p(i,j,k,2)*
     1   (u(i,j,k,5)+u(i,j-1,k,5)+u(i,j,k-1,5)+u(i,j-1,k-1,5)
     2   -u(i-1,j,k,5)-u(i-1,j-1,k,5)-u(i-1,j,k-1,5)-u(i-1,j-1,k-1,5))
     3   -dy2*p(i,j,k,3)*
     3   (u(i,j,k,5)-u(i,j-1,k,5)+u(i,j,k-1,5)-u(i,j-1,k-1,5)
     4   +u(i-1,j,k,5)-u(i-1,j-1,k,5)+u(i-1,j,k-1,5)-u(i-1,j-1,k-1,5))
     5   -dz2*p(i,j,k,4)*
     5   (u(i,j,k,5)+u(i,j-1,k,5)-u(i,j,k-1,5)-u(i,j-1,k-1,5)
     6   +u(i-1,j,k,5)+u(i-1,j-1,k,5)-u(i-1,j,k-1,5)-u(i-1,j-1,k-1,5))
     7   +pmu*(dx4*(ff(i+1,j,k,5)-2.0*ff(i,j,k,5)+ff(i-1,j,k,5))
     8   +dy4*(ff(i,j+1,k,5)-2.0*ff(i,j,k,5)+ff(i,j-1,k,5))
     9   +dz4*(ff(i,j,k+1,5)-2.0*ff(i,j,k,5)+ff(i,j,k-1,5)))
c
       f(i,j,k,4)=ff(i,j,k,4)+t*((p(i,j,k,9)*p(i,j,k,7)
     1   -p(i,j,k,10)*p(i,j,k,6))/p(i,j,k,1)
     1   -0.0*eud*p(i,j,k,4)-gra*z/ro3)
     2   -dx2*p(i,j,k,2)*
     2   (u(i,j,k,4)+u(i,j-1,k,4)+u(i,j,k-1,4)+u(i,j-1,k-1,4)
     3   -u(i-1,j,k,4)-u(i-1,j-1,k,4)-u(i-1,j,k-1,4)-u(i-1,j-1,k-1,4))
     4   -dy2*p(i,j,k,3)*
     4   (u(i,j,k,4)-u(i,j-1,k,4)+u(i,j,k-1,4)-u(i,j-1,k-1,4)
     5   +u(i-1,j,k,4)-u(i-1,j-1,k,4)+u(i-1,j,k-1,4)-u(i-1,j-1,k-1,4))
     6   -dz2*p(i,j,k,4)*
     6   (u(i,j,k,4)+u(i,j-1,k,4)-u(i,j,k-1,4)-u(i,j-1,k-1,4)
     7   +u(i-1,j,k,4)+u(i-1,j-1,k,4)-u(i-1,j,k-1,4)-u(i-1,j-1,k-1,4))
     8   -dz2*(u(i,j,k,5)+u(i,j-1,k,5)-u(i,j,k-1,5)-u(i,j-1,k-1,5)
     9   +u(i-1,j,k,5)+u(i-1,j-1,k,5)-u(i-1,j,k-1,5)-u(i-1,j-1,k-1,5))
     9   /p(i,j,k,1)
     1   +rmu*(dx4*(ff(i+1,j,k,4)-2.0*ff(i,j,k,4)+ff(i-1,j,k,4))
     2   +dy4*(ff(i,j+1,k,4)-2.0*ff(i,j,k,4)+ff(i,j-1,k,4))
     3   +dz4*(ff(i,j,k+1,4)-2.0*ff(i,j,k,4)+ff(i,j,k-1,4)))/ff(i,j,k,1)
c
       f(i,j,k,3)=ff(i,j,k,3)+t*((p(i,j,k,11)*p(i,j,k,6)
     1   -p(i,j,k,9)*p(i,j,k,8))/p(i,j,k,1)
     1   -0.0*eud*p(i,j,k,3)-gra*y/ro3)
     2   -dx2*p(i,j,k,2)*
     2   (u(i,j,k,3)+u(i,j-1,k,3)+u(i,j,k-1,3)+u(i,j-1,k-1,3)
     3   -u(i-1,j,k,3)-u(i-1,j-1,k,3)-u(i-1,j,k-1,3)-u(i-1,j-1,k-1,3))
     4   -dy2*p(i,j,k,3)*
     4   (u(i,j,k,3)-u(i,j-1,k,3)+u(i,j,k-1,3)-u(i,j-1,k-1,3)
     5   +u(i-1,j,k,3)-u(i-1,j-1,k,3)+u(i-1,j,k-1,3)-u(i-1,j-1,k-1,3))
     6   -dz2*p(i,j,k,4)*
     6   (u(i,j,k,3)+u(i,j-1,k,3)-u(i,j,k-1,3)-u(i,j-1,k-1,3)
     7   +u(i-1,j,k,3)+u(i-1,j-1,k,3)-u(i-1,j,k-1,3)-u(i-1,j-1,k-1,3))
     8   -dy2*(u(i,j,k,5)-u(i,j-1,k,5)+u(i,j,k-1,5)-u(i,j-1,k-1,5)
     9   +u(i-1,j,k,5)-u(i-1,j-1,k,5)+u(i-1,j,k-1,5)-u(i-1,j-1,k-1,5))
     9   /p(i,j,k,1)
     1   +rmu*(dx4*(ff(i+1,j,k,3)-2.0*ff(i,j,k,3)+ff(i-1,j,k,3))
     2   +dy4*(ff(i,j+1,k,3)-2.0*ff(i,j,k,3)+ff(i,j-1,k,3))
     3   +dz4*(ff(i,j,k+1,3)-2.0*ff(i,j,k,3)+ff(i,j,k-1,3)))/ff(i,j,k,1)
c
       uxd=amax1((x-4.0*ar1)/(4.0*ar1),0.0)
       uxd=amin1(uxd,1.0)
       uxd=amin1(uxd*p(i,j,k,2),0.0)
c
       f(i,j,k,2)=ff(i,j,k,2)+t*((p(i,j,k,10)*p(i,j,k,8)
     1   -p(i,j,k,11)*p(i,j,k,7))/p(i,j,k,1)
     1   -eud*uxd-gra*x/ro3)
     2   -dx2*p(i,j,k,2)*
     2   (u(i,j,k,2)+u(i,j-1,k,2)+u(i,j,k-1,2)+u(i,j-1,k-1,2)
     3   -u(i-1,j,k,2)-u(i-1,j-1,k,2)-u(i-1,j,k-1,2)-u(i-1,j-1,k-1,2))
     4   -dy2*p(i,j,k,3)*
     4   (u(i,j,k,2)-u(i,j-1,k,2)+u(i,j,k-1,2)-u(i,j-1,k-1,2)
     5   +u(i-1,j,k,2)-u(i-1,j-1,k,2)+u(i-1,j,k-1,2)-u(i-1,j-1,k-1,2))
     6   -dz2*p(i,j,k,4)*
     6   (u(i,j,k,2)+u(i,j-1,k,2)-u(i,j,k-1,2)-u(i,j-1,k-1,2)
     7   +u(i-1,j,k,2)+u(i-1,j-1,k,2)-u(i-1,j,k-1,2)-u(i-1,j-1,k-1,2))
     8   -dx2*(u(i,j,k,5)+u(i,j-1,k,5)+u(i,j,k-1,5)+u(i,j-1,k-1,5)
     9   -u(i-1,j,k,5)-u(i-1,j-1,k,5)-u(i-1,j,k-1,5)-u(i-1,j-1,k-1,5))
     9   /p(i,j,k,1)
     1   +rmu*(dx4*(ff(i+1,j,k,2)-2.0*ff(i,j,k,2)+ff(i-1,j,k,2))
     2   +dy4*(ff(i,j+1,k,2)-2.0*ff(i,j,k,2)+ff(i,j-1,k,2))
     3   +dz4*(ff(i,j,k+1,2)-2.0*ff(i,j,k,2)+ff(i,j,k-1,2)))/ff(i,j,k,1)
c
       f(i,j,k,1)=ff(i,j,k,1)
     1   -dx2*(u(i,j,k,2)*u(i,j,k,1)+u(i,j-1,k,2)*u(i,j-1,k,1)
     2   +u(i,j,k-1,2)*u(i,j,k-1,1)+u(i,j-1,k-1,2)*u(i,j-1,k-1,1)
     3   -u(i-1,j,k,2)*u(i-1,j,k,1)-u(i-1,j-1,k,2)*u(i-1,j-1,k,1)
     4   -u(i-1,j,k-1,2)*u(i-1,j,k-1,1)
     4   -u(i-1,j-1,k-1,2)*u(i-1,j-1,k-1,1))
     5   -dy2*(u(i,j,k,3)*u(i,j,k,1)-u(i,j-1,k,3)*u(i,j-1,k,1)
     6   +u(i,j,k-1,3)*u(i,j,k-1,1)-u(i,j-1,k-1,3)*u(i,j-1,k-1,1)
     7   +u(i-1,j,k,3)*u(i-1,j,k,1)-u(i-1,j-1,k,3)*u(i-1,j-1,k,1)
     8   +u(i-1,j,k-1,3)*u(i-1,j,k-1,1)
     8   -u(i-1,j-1,k-1,3)*u(i-1,j-1,k-1,1))
     1   -dz2*(u(i,j,k,4)*u(i,j,k,1)+u(i,j-1,k,4)*u(i,j-1,k,1)
     2   -u(i,j,k-1,4)*u(i,j,k-1,1)-u(i,j-1,k-1,4)*u(i,j-1,k-1,1)
     3   +u(i-1,j,k,4)*u(i-1,j,k,1)+u(i-1,j-1,k,4)*u(i-1,j-1,k,1)
     4   -u(i-1,j,k-1,4)*u(i-1,j,k-1,1)
     4   -u(i-1,j-1,k-1,4)*u(i-1,j-1,k-1,1))
     5   +dfu*(dx4*(ff(i+1,j,k,1)-2.0*ff(i,j,k,1)+ff(i-1,j,k,1))
     6   +dy4*(ff(i,j+1,k,1)-2.0*ff(i,j,k,1)+ff(i,j-1,k,1))
     7   +dz4*(ff(i,j,k+1,1)-2.0*ff(i,j,k,1)+ff(i,j,k-1,1)))
c
  801  continue
c
       do 820 i=2,nx1
       x=hhx(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(i,j,k,8)=f(i,j,k,8)*x2+ff(i,j,k,8)*(1.0-x2)
       f(i,j,k,7)=f(i,j,k,7)*x2+ff(i,j,k,7)*(1.0-x2)
       f(i,j,k,6)=f(i,j,k,6)*x2+ff(i,j,k,6)*(1.0-x2)
       f(i,j,k,5)=f(i,j,k,5)*x2+ff(i,j,k,5)*(1.0-x2)
       f(i,j,k,4)=f(i,j,k,4)*x2+ff(i,j,k,4)*(1.0-x2)
       f(i,j,k,3)=f(i,j,k,3)*x2+ff(i,j,k,3)*(1.0-x2)
       f(i,j,k,2)=f(i,j,k,2)*x2+ff(i,j,k,2)*(1.0-x2)
       f(i,j,k,1)=f(i,j,k,1)*x2+ff(i,j,k,1)*(1.0-x2)
  820  continue
c
       do 824 i=2,nx1
       x5=f(i,j,k,5)
       x1=f(i,j,k,1)
       y5=amax1(x5,pr02)
       y1=amax1(x1,ro02)
       x2=f(i,j,k,2)*f(i,j,k,2)+f(i,j,k,3)*f(i,j,k,3)
     1    +f(i,j,k,4)*f(i,j,k,4)
       vmax=amax1(vmax,x2)
       x=amax1(x1,0.0)/y1
       f(i,j,k,5)=y5
       f(i,j,k,4)=x*f(i,j,k,4)
       f(i,j,k,3)=x*f(i,j,k,3)
       f(i,j,k,2)=x*f(i,j,k,2)
       f(i,j,k,1)=y1
  824  continue
   80  continue
c
!hpf$  end on
  190  continue
c      !xocl end spread max(vmax)
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(ii.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(ii)+ths0*float(last)
c      th=orid*th*pi/180.0
       th=orid*pi/180.0
       ith=ii+last*(iii-1)
       th1=180.0*float(ith)/float(last*4)
       th2=orid+th1
       if(th2.gt.360.0) th2=th2-360.0
c      th=th2*pi/180.0
       bisy=bis*cos(th)
       bisz=bis*sin(th)
c
       iino=ii+last*(iii-1)+iinob
       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
!hpf$  independent,new(i,j,k)
CC MPI START
       do 102 k=ks,ke
CC MPI END
!hpf$  on home(f(:,:,k,:)),local(f,i,j,m,nx11,x,x1,hx,hhx,
!hpf$*         bis,bisz,bisy,pro1,vswx,roo1) begin
       do 1021 j=1,ny2
       do 1021 i=1,nx2
       x=(2.0*hhx(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(i,j,k,8)=(1.0-x)*f(i,j,k,8)+bis*x
       f(i,j,k,8)=(1.0-x)*f(i,j,k,8)+bisz*x
       f(i,j,k,7)=(1.0-x)*f(i,j,k,7)+bisy*x
       f(i,j,k,6)=(1.0-x)*f(i,j,k,6)
       f(i,j,k,5)=(1.0-x)*f(i,j,k,5)+pro1*x
       f(i,j,k,4)=(1.0-x)*f(i,j,k,4)
       f(i,j,k,3)=(1.0-x)*f(i,j,k,3)
       f(i,j,k,2)=(1.0-x)*f(i,j,k,2)+vswx*x
       f(i,j,k,1)=(1.0-x)*f(i,j,k,1)+roo1*x
 1021  continue
!hpf$  end on
  102  continue
c
!hpf$  independent,new(i,j,k)
CC MPI START
       do 104 k=ks,ke1
CC MPI END
!hpf$  on home(u(:,:,k,:)),local(u,i,j,m,nx11,x,x1,hx,hhx,
!hpf$*         bis,bisz,bisy,pro1,vswx,roo1) begin
       do 1041 j=1,ny1
       do 1041 i=1,nx1
       x=(2.0*hhx(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(i,j,k,8)=(1.0-x)*u(i,j,k,8)+bis*x
       u(i,j,k,8)=(1.0-x)*u(i,j,k,8)+bisz*x
       u(i,j,k,7)=(1.0-x)*u(i,j,k,7)+bisy*x
       u(i,j,k,6)=(1.0-x)*u(i,j,k,6)
       u(i,j,k,5)=(1.0-x)*u(i,j,k,5)+pro1*x
       u(i,j,k,4)=(1.0-x)*u(i,j,k,4)
       u(i,j,k,3)=(1.0-x)*u(i,j,k,3)
       u(i,j,k,2)=(1.0-x)*u(i,j,k,2)+vswx*x
       u(i,j,k,1)=(1.0-x)*u(i,j,k,1)+roo1*x
 1041  continue
!hpf$  end on
  104  continue
c
c
  403  if(iis.ne.iis0) go to 100
       iis=0
  402  continue
       call clock(zt2)
c      if(ii.lt.1024) go to 100
c
c
c
c!hpf$  asynchronous(id1),nobuffer begin
c       gf(1:nx2,1:ny2,1:nz2,1:nb)=f(1:nx2,1:ny2,1:nz2,1:nb)
c!hpf$  end asynchronous
c!hpf$  asyncwait(id1)
c
c      write(ntap) f
       do m=1,nb
CC MPI START
         call mpi_gatherv(f(1,1,1,m),nword,mpi_real,fg(1,1,1),
     *                    recvcount,displs,mpi_real,0,
     *                    mpi_comm_world,ier)
       if(irank.eq.0) then
       do 1732 k=1,nz2
cc        write(ntap) fg(1:nx2,1:ny2,k)
 1732  continue
       end if
CC MPI END
       end do
c
c      do 173 k=1,nz2
c!hpf$  asynchronous(id1),nobuffer begin
c      gf(1:nx2,1:ny2)=f(1:nx2,1:ny2,k,m)
c!hpf$  end asynchronous
c!hpf$  asyncwait(id1)
c      write(ntap) gf
cc       write(ntap) f(1:nx2,1:ny2,k,m)
c 173  continue
c
CC MPI START
       if(irank.eq.0)
     * write(6,661) ii,last,iii,ith,nin,itap,th2,th,x2,vmax
CC MPI END
  661  format(1h ,1x,6i10/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) ii,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!xocl end parallel
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)
       parameter (npe=2)
!hpf$  processors pe(npe)
c      parameter (nx=320,ny= 80,nz=160,nb=8)
       parameter (nx=500,ny=100,nz=200,nb=8)
c      parameter (nx=800,ny=200,nz=478,nb=8)
c      parameter (nx=800,ny=200,nz=670,nb=8)
c      parameter (nx=1000,ny= 500,nz=1118,nb=8)
c      parameter (nx=1000,ny=1000,nz=1118,nb=8)
c      parameter (nx=1000,ny=1000,nz=1258,nb=8)
c      parameter (nx=1678,ny= 558,nz=1118,nb=8)
c      parameter (nx=2238,ny= 558,nz=1118,nb=8)
       parameter (nx1=nx+1,nx2=nx+2,ny1=ny+1,ny2=ny+2)
       parameter (nz1=nz+1,nz2=nz+2)
c
CC MPI START
       parameter(nzz=(nz2-1)/npe+1)
       dimension f(nx2,ny2,0:nzz+1,nb),pp(nx2,ny2,0:nzz+1,3)
CC MPI END
c      dimension gf(nx2,ny2,nz2,nb),gpp(nx2,ny2,nz2,3)
c      dimension gf(nx2,ny2,nz2,nb)
       dimension cj(1),cp(1)
!hpf$  distribute f(*,*,block,*) onto pe
!hpf$  distribute pp(*,*,block,*) onto pe
!hpf$  shadow     f(0,0,1:1,0)
!hpf$  shadow     pp(0,0,1:1,0)
c      !xocl global gf,gpp
c      equivalence (gf,f),(gpp,pp)
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
CC MPI START
       common /para_info/ks,ks1,ke,ke1,kss,irank,isize
CC MPI END
       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)
c
!hpf$  independent,new(i,j,k)
CC MPI START
       do 10 k=ks,ke
CC MPI END
!hpf$  on home(f(:,:,k,:)),local(f,i,j,m,nxp,x,y,z,
!hpf$*        dn1,dv1,ax1,hx,hy,hz,dn,vr1,vsw,ro1,
!hpf$*        ro2,ro3,x1,ro02,cp,cj,p0,ro1,pr02,px1,
!hpf$*        y1,xx,bis,ro01) begin
       do 101 j=1,ny2
       vr1=vsw
       do 101 i=1,nx2
       x=0.5*hx*float(2*i-nx2-1+2*nxp)
       y=0.5*hy*float(2*j-3)
       z=0.5*hz*float(2*(kss+k)-nz2-1)
       ro2=x*x+y*y+z*z
       ro1=sqrt(ro2)
       ro3=ro1*ro2
       f(i,j,k,1)=1.0/ro3
c      f(i,j,k,1)=ro1**(-x1)
       if(f(i,j,k,1).lt.ro02) f(i,j,k,1)=ro02
       f(i,j,k,5)=1.0e-7*cp(6)/ro2
c      f(i,j,k,5)=p0*ro1**(-x2)
       if(f(i,j,k,5).lt.pr02) f(i,j,k,5)=pr02
       f(i,j,k,6)=-2.0*b0*x*z/(ro2*ro2)
       f(i,j,k,7)=-2.0*b0*y*z/(ro2*ro2)
       f(i,j,k,8)=b0*(x*x+y*y-z*z)/(ro2*ro2)
       f(i,j,k,6)=-3.0*b0*x*z/(ro2*ro3)
       f(i,j,k,7)=-3.0*b0*y*z/(ro2*ro3)
       f(i,j,k,8)=b0*(x*x+y*y-2.0*z*z)/(ro2*ro3)
       px1=sqrt(f(i,j,k,6)*f(i,j,k,6)+f(i,j,k,7)*f(i,j,k,7)
     1         +f(i,j,k,8)*f(i,j,k,8))
       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(i,j,k,2)=vr1
       xx=vr1/vsw
       f(i,j,k,5)=f(i,j,k,5)+(pr01-pr02)*xx
       f(i,j,k,8)=f(i,j,k,8)+bis
  101  continue
!hpf$  end on
   10  continue
c
c
c      current
!hpf$  reflect f
CC MPI START
       do m=1,nb
c        len=nx2*ny2
         if (irank.gt.0) then
           call mpi_send(f(1,1,ks,m),n2,mpi_real,irank-1,
     &                   140,mpi_comm_world,ier)
         end if
         if (irank.lt.isize-1) then
           call mpi_recv(f(1,1,ke+1,m),n2,mpi_real,irank+1,
     &                   140,mpi_comm_world,istatus,ier)
         end if
       end do
       call mpi_barrier(mpi_comm_world,ier)
CC MPI END
c
!hpf$  independent,new(i,j,k)
CC MPI START
       do 38 k=ks,ke1
CC MPI END
!hpf$  on home(pp(:,:,k,:)),local(pp,f,i,j,m) begin
       do 381 j=1,ny1
       do 381 i=1,nx1
       pp(i,j,k,3)=0.25*((f(i+1,j+1,k+1,7)+f(i+1,j,k+1,7)
     1   +f(i+1,j+1,k,7)+f(i+1,j,k,7)
     2   -f(i,j+1,k+1,7)-f(i,j,k+1,7)-f(i,j+1,k,7)-f(i,j,k,7))/hx
     3   -(f(i+1,j+1,k+1,6)-f(i+1,j,k+1,6)+f(i+1,j+1,k,6)-f(i+1,j,k,6)
     4   +f(i,j+1,k+1,6)-f(i,j,k+1,6)+f(i,j+1,k,6)-f(i,j,k,6))/hy)
       pp(i,j,k,2)=0.25*((f(i+1,j+1,k+1,6)+f(i+1,j,k+1,6)
     1   -f(i+1,j+1,k,6)-f(i+1,j,k,6)
     2   +f(i,j+1,k+1,6)+f(i,j,k+1,6)-f(i,j+1,k,6)-f(i,j,k,6))/hz
     3   -(f(i+1,j+1,k+1,8)+f(i+1,j,k+1,8)+f(i+1,j+1,k,8)+f(i+1,j,k,8)
     4   -f(i,j+1,k+1,8)-f(i,j,k+1,8)-f(i,j+1,k,8)-f(i,j,k,8))/hx)
       pp(i,j,k,1)=0.25*((f(i+1,j+1,k+1,8)-f(i+1,j,k+1,8)
     1   +f(i+1,j+1,k,8)-f(i+1,j,k,8)
     2   +f(i,j+1,k+1,8)-f(i,j,k+1,8)+f(i,j+1,k,8)-f(i,j,k,8))/hy
     3   -(f(i+1,j+1,k+1,7)+f(i+1,j,k+1,7)-f(i+1,j+1,k,7)-f(i+1,j,k,7)
     4   +f(i,j+1,k+1,7)+f(i,j,k+1,7)-f(i,j+1,k,7)-f(i,j,k,7))/hz)
  381  continue
!hpf$  end on
   38  continue
       return
       end
       subroutine equib8(f,pp,nxp,ro01,pr01,rrat,cj,cp)
c      implicit real*8 (a-h,o-z)
       parameter (npe=2)
!hpf$  processors pe(npe)
c      parameter (nx=320,ny= 80,nz=160,nb=8)
       parameter (nx=500,ny=100,nz=200,nb=8)
c      parameter (nx=800,ny=200,nz=478,nb=8)
c      parameter (nx=800,ny=200,nz=670,nb=8)
c      parameter (nx=1000,ny= 500,nz=1118,nb=8)
c      parameter (nx=1000,ny=1000,nz=1118,nb=8)
c      parameter (nx=1000,ny=1000,nz=1258,nb=8)
c      parameter (nx=1678,ny= 558,nz=1118,nb=8)
c      parameter (nx=2238,ny= 558,nz=1118,nb=8)
       parameter (nx1=nx+1,nx2=nx+2,ny1=ny+1,ny2=ny+2)
       parameter (nz1=nz+1,nz2=nz+2)
c
CC MPI START
       parameter(nzz=(nz2-1)/npe+1)
       dimension f(nx2,ny2,0:nzz+1,nb),pp(nx2,ny2,0:nzz+1,3)
CC MPI END
c      dimension gf(nx2,ny2,nz2,nb),gpp(nx2,ny2,nz2,3)
c      dimension gf(nx2,ny2,nz2,nb)
       dimension cj(1),cp(1)
!hpf$  distribute f(*,*,block,*) onto pe
!hpf$  distribute pp(*,*,block,*) onto pe
!hpf$  shadow     f(0,0,1:1,0)
!hpf$  shadow     pp(0,0,1:1,0)
c      !xocl global gf,gpp
c      equivalence (gf,f),(gpp,pp)
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
CC MPI START
       include 'mpif.h'
       integer istatus(mpi_status_size)
       common /para_info/ks,ks1,ke,ke1,kss,irank,isize
CC MPI END
       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)
       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
!hpf$  independent,new(i,j,k)
CC MPI START
       do 10 k=ks,ke
CC MPI END
!hpf$  on home(f(:,:,k,:)),local(f,i,j,m,nxp,x,y,z,
!hpf$*         dn1,dv1,ax1,vr1,vsw,x1m,xm,ro1,ro2,ro3,
!hpf$*        hx,hy,hz,dn,ro1m,ro2m,ro3m,x1,ro02,pr02,
!hpf$*        cp,cj,p0,b0) begin
       do 101 j=1,ny2
       vr1=vsw
       do 101 i=1,nx2
       x=0.5*hx*float(2*i-nx2-1+2*nxp)
       y=0.5*hy*float(2*j-3)
       z=0.5*hz*float(2*(kss+k)-nz2-1)
       x1m=x-2.0*xm
       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(i,j,k,1)=1.0/ro3
c      f(i,j,k,1)=ro1**(-x1)
       if(f(i,j,k,1).lt.ro02) f(i,j,k,1)=ro02
       f(i,j,k,5)=1.0e-7*cp(6)/ro2
c      f(i,j,k,5)=p0*ro1**(-x2)
       if(f(i,j,k,5).lt.pr02) f(i,j,k,5)=pr02
       f(i,j,k,6)=-2.0*b0*x*z/(ro2*ro2)
       f(i,j,k,7)=-2.0*b0*y*z/(ro2*ro2)
       f(i,j,k,8)=b0*(x*x+y*y-z*z)/(ro2*ro2)
       f(i,j,k,6)=-3.0*b0*x*z/(ro2*ro3)
       f(i,j,k,7)=-3.0*b0*y*z/(ro2*ro3)
       f(i,j,k,8)=b0*(x*x+y*y-2.0*z*z)/(ro2*ro3)
  101  continue
!hpf$  end on
   10  continue
c
c
c      current
!hpf$  reflect f
CC MPI START
       do m=1,nb
c        len=nx2*ny2
         if (irank.gt.0) then
           call mpi_send(f(1,1,ks,m),n2,mpi_real,irank-1,
     &                   140,mpi_comm_world,ier)
         end if
         if (irank.lt.isize-1) then
           call mpi_recv(f(1,1,ke+1,m),n2,mpi_real,irank+1,
     &                   140,mpi_comm_world,istatus,ier)
         end if
       end do
       call mpi_barrier(mpi_comm_world,ier)
CC MPI END
c
!hpf$  independent,new(i,j,k)
CC MPI START
       do 38 k=ks,ke1
CC MPI END
!hpf$  on home(pp(:,:,k,:)),local(pp,f,i,j,m,hx,hy,hz) begin
       do 381 j=1,ny1
       do 381 i=1,nx1
       pp(i,j,k,3)=0.25*((f(i+1,j+1,k+1,7)+f(i+1,j,k+1,7)
     1   +f(i+1,j+1,k,7)+f(i+1,j,k,7)
     2   -f(i,j+1,k+1,7)-f(i,j,k+1,7)-f(i,j+1,k,7)-f(i,j,k,7))/hx
     3   -(f(i+1,j+1,k+1,6)-f(i+1,j,k+1,6)+f(i+1,j+1,k,6)-f(i+1,j,k,6)
     4   +f(i,j+1,k+1,6)-f(i,j,k+1,6)+f(i,j+1,k,6)-f(i,j,k,6))/hy)
       pp(i,j,k,2)=0.25*((f(i+1,j+1,k+1,6)+f(i+1,j,k+1,6)
     1   -f(i+1,j+1,k,6)-f(i+1,j,k,6)
     2   +f(i,j+1,k+1,6)+f(i,j,k+1,6)-f(i,j+1,k,6)-f(i,j,k,6))/hz
     3   -(f(i+1,j+1,k+1,8)+f(i+1,j,k+1,8)+f(i+1,j+1,k,8)+f(i+1,j,k,8)
     4   -f(i,j+1,k+1,8)-f(i,j,k+1,8)-f(i,j+1,k,8)-f(i,j,k,8))/hx)
       pp(i,j,k,1)=0.25*((f(i+1,j+1,k+1,8)-f(i+1,j,k+1,8)
     1   +f(i+1,j+1,k,8)-f(i+1,j,k,8)
     2   +f(i,j+1,k+1,8)-f(i,j,k+1,8)+f(i,j+1,k,8)-f(i,j,k,8))/hy
     3   -(f(i+1,j+1,k+1,7)+f(i+1,j,k+1,7)-f(i+1,j+1,k,7)-f(i+1,j,k,7)
     4   +f(i,j+1,k+1,7)+f(i,j,k+1,7)-f(i,j+1,k,7)-f(i,j,k,7))/hz)
  381  continue
!hpf$  end on
   38  continue
c
c
!hpf$  independent,new(i,j,k)
CC MPI START
       do 40 k=ks,ke
CC MPI END
!hpf$  on home(f(:,:,k,:)),local(f,i,j,m,nxp,x,y,z,
!hpf$*         dn1,dv1,ax1,vr1,vsw,x1m,xm,ro1,ro2,ro3,
!hpf$*        hx,hy,hz,dn,ro1m,ro2m,ro3m,x1,ro02,pr02,
!hpf$*        cp,cj,p0,b0,bx1,bx2,by1,by2,bz1,bz2,y1,
!hpf$*        alp1,bet,px1,ro01,pr01,xx,bis) begin
       do 401 j=1,ny2
c      vr1=vsw
       do 401 i=1,nx2
       x=0.5*hx*float(2*i-nx2-1+2*nxp)
       y=0.5*hy*float(2*j-3)
CC MPI START
       z=0.5*hz*float(2*(kss+k)-nz2-1)
C MPI END
       x1m=x-2.0*xm
       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(i,j,k,1)=1.0/ro3
c      f(i,j,k,1)=ro1**(-x1)
       if(f(i,j,k,1).lt.ro02) f(i,j,k,1)=ro02
       f(i,j,k,5)=1.0e-7*cp(6)/ro2
c      f(i,j,k,5)=p0*ro1**(-x2)
       if(f(i,j,k,5).lt.pr02) f(i,j,k,5)=pr02
       f(i,j,k,6)=-2.0*b0*x*z/(ro2*ro2)
       f(i,j,k,7)=-2.0*b0*y*z/(ro2*ro2)
       f(i,j,k,8)=b0*(x*x+y*y-z*z)/(ro2*ro2)
       bx1=-3.0*b0*x*z/(ro2*ro3)
       by1=-3.0*b0*y*z/(ro2*ro3)
       bz1=b0*(x*x+y*y-2.0*z*z)/(ro2*ro3)
       bx2=-3.0*b0*x1m*z/(ro2m*ro3m)
       by2=-3.0*b0*y*z/(ro2m*ro3m)
       bz2=b0*(x1m*x1m+y*y-2.0*z*z)/(ro2m*ro3m)
       f(i,j,k,6)=bx1+bx2
       f(i,j,k,7)=by1+by2
       f(i,j,k,8)=bz1+bz2
       if(x.lt.xm) f(i,j,k,6)=0.0
       if(x.lt.xm) f(i,j,k,7)=0.0
       if(x.lt.xm) f(i,j,k,8)=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(i,j,k,6)*f(i,j,k,6)+f(i,j,k,7)*f(i,j,k,7)
     1         +f(i,j,k,8)*f(i,j,k,8))
       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(i,j,k,2)=vr1
       xx=vr1/vsw
       f(i,j,k,5)=f(i,j,k,5)+(pr01-pr02)*xx
       f(i,j,k,8)=f(i,j,k,8)+bis*y1
  401  continue
!hpf$  end on
   40  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