c
       program hdiff3
c  
c      main program of 3-dimensional difussion equation
c      modified leap-frog method  2000.06.10
       implicit real*8 (a-h,o-z)
       parameter(npe=16)
!hpf$  processors pe(npe)
c      parameter(nx=30,ny=30,nz=30)
c      parameter(nx=100,ny=100,nz=100)
       parameter(nx=254,ny=254,nz=254)
       parameter(last=4,iip0=8)
       parameter(nx1=nx+1,nx2=nx+2,ny1=ny+1,ny2=ny+2)
       parameter(nz1=nz+1,nz2=nz+2)
       parameter(n1=nx2,n2=n1*ny2,n3=n2*nz2,noinp=30)
       parameter(thx=0.5,tam=0.5,dif0=0.10)
c
       dimension f(nx2,ny2,nz2),u(nx2,ny2,nz2)
c      dimension gf(nx2,ny2,nz2),gu(nx2,ny2,nz2)
       dimension rinp(noinp)
       dimension ppin(10)
c      integer iinp(noinp)
!hpf$  distribute f(*,*,block) onto pe
!hpf$  distribute u(*,*,block) onto pe
!hpf$  shadow     u(0,0,1:1)
!hpf$  asyncid id1
c      equivalence (gf,f),(gu,u)
       common /blk/f
c
c           ppin xl   yl   zl   dxl  dyl  dzl  dn  dv   9  10
       data ppin/62.0,62.0,62.0,10.0,10.0,10.0,1.0,0.0, 9, 10/
c
c
       xl=ppin(1)
       yl=ppin(2)
       zl=ppin(3)
       dxl=ppin(4)
       dyl=ppin(5)
       dzl=ppin(6)
       dn=ppin(7)
       dv=ppin(8)
       dif=dif0
c
       hx=xl/float(nx1)
       hy=yl/float(ny1)
       hz=zl/float(nz1)
       t1=thx*hx
       dx1=0.5*t1/hx
       dy1=0.5*t1/hy
       dz1=0.5*t1/hz
       dx2=dif*t1/(hx*hx)
       dy2=dif*t1/(hy*hy)
       dz2=dif*t1/(hz*hz)
c
       rinp(1)=xl
       rinp(2)=yl
       rinp(3)=zl
       rinp(4)=hx
       rinp(5)=hy
       rinp(6)=hz
       rinp(7)=t1
       rinp(8)=dx1
       rinp(9)=dy1
       rinp(10)=dz1
       rinp(11)=dx2
       rinp(12)=dy2
       rinp(13)=dz2
c
c      initial parameters
c
       write(6,122) nx,ny,nz,last,nx2,ny2,nz2
       write(6,124) t1,hx,hy,hz,dx1,dy1,dz1
  122  format(1h ,10i8)
  124  format(1h ,8(1pe10.3))
c
c      initial conditions
       call ainit1(ppin)
c
c      start of calculation
c
       call clock(zt0)
c
       do 500 ii=1,last
c
       call clock(zt1)
c
       do 100 iip=1,iip0
c
c      boundary conditions
c
c      periodic boundary in x-direction
!hpf$  independent,new(j,k)
       do 22 k=2,nz1
!hpf$  on home(f(:,:,k)),local(f,j) begin
       do 221 j=2,ny1
       f(1,j,k)=f(nx1,j,k)
       f(nx2,j,k)=f(2,j,k)
  221  continue
!hpf$  end on
   22  continue
c
c      periodic boundary in y-direction
!hpf$  independent,new(i,k)
       do 24 k=2,nz1
!hpf$  on home(f(:,:,k)),local(f,i) begin
       do 241 i=1,nx2
       f(i,1,k)=f(i,ny1,k)
       f(i,ny2,k)=f(i,2,k)
  241  continue
!hpf$  end on
   24  continue
c
c      periodic boundary in z-direction
c      do 26 k=1,1
c      do 26 j=1,ny2
c      do 26 i=1,nx2
c      f(i,j,1)=f(i,j,nz1)
c      f(i,j,nz2)=f(i,j,2)
c      f(i,j,k)=f(i,j,k+nz)
c  26  continue
c
c      do 28 k=nz2,nz2
c      do 28 j=1,ny2
c      do 28 i=1,nx2
c      f(i,j,1)=gf(i,j,nz1)
c      f(i,j,nz2)=gf(i,j,2)
c      f(i,j,k)=f(i,j,k-nz)
c  28  continue
c
!hpf$ asynchronous(id1),nobuffer begin
      f(1:nx2,1:ny2,  1) = f(1:nx2,1:ny2,1+nz)
      f(1:nx2,1:ny2,nz2) = f(1:nx2,1:ny2,   2)
!hpf$ end asynchronous
!hpf$ asyncwait(id1)
c
c
c      do 60 k=1,nz2
c      do 601 j=1,ny2
c      do 601 i=1,nx2
c      u(i,j,k)=f(i,j,k)
c 601  continue
c  60  continue
c
!hpf$ asynchronous(id1),nobuffer begin
      u(1:nx2,1:ny2,1:nz2) = f(1:nx2,1:ny2,1:nz2)
!hpf$ end asynchronous
!hpf$ asyncwait(id1)
c
c
!hpf$  reflect u
!hpf$  independent,new(i,j,k)
       do 80 k=2,nz1
!hpf$  on home(f(:,:,k)),local(f,u,i,j,dx2,dy2,dz2) begin
       do 801 j=2,ny1
       do 801 i=2,nx1
c
       f(i,j,k)=f(i,j,k)
     1   +dx2*(u(i-1,j,k)-2.0*u(i,j,k)+u(i+1,j,k))
     2   +dy2*(u(i,j-1,k)-2.0*u(i,j,k)+u(i,j+1,k))
     3   +dz2*(u(i,j,k-1)-2.0*u(i,j,k)+u(i,j,k+1))
c
  801  continue
!hpf$  end on
   80  continue
c
c
  100  continue
c
c
       call clock(zt2)
       zt1=zt1-zt0
       zt2=zt2-zt0
       zt=zt2-zt1
       write(6,402) ii,zt0,zt1,zt2,zt
  402  format(1h , i6,1pe12.3,3(0pf12.5))
c
c      write the output data
c
  500  continue
    9  continue
c
       stop
       end
       subroutine clock(ti)
       implicit real*8 (a-h,o-z)
c      real*8 ti,ti1
       ti=1.0d0
       ti1=1.0d0
c      call gettod(ti1)
       call fjhpf_gettod(ti1)
       ti=1.0d-6*ti1
c      x=0.0
c      y=secnds(x)
c      ti=1.0d0*y
       return
       end
       subroutine clock2(ti)
       implicit real*8 (a-h,o-z)
c       integer ati,time
       dimension tt(2)
c       ati=time()
       x=0.0
c       y=secnds(x)
c       ti=1.00000*dfloat(ati)
       y=etime(tt)
       ti=1.000d0*tt(1)
       return
       end
       subroutine ainit1(ppin)
       implicit real*8 (a-h,o-z)
       parameter(npe=16)
!hpf$  processors pe(npe)
c      parameter(nx=30,ny=30,nz=30)
c      parameter(nx=100,ny=100,nz=100)
       parameter(nx=254,ny=254,nz=254)
       parameter(iip0=8,last=4)
       parameter(nx1=nx+1,nx2=nx+2,ny1=ny+1,ny2=ny+2)
       parameter(nz1=nz+1,nz2=nz+2)
       parameter(n1=nx2,n2=n1*ny2,n3=n2*nz2,noinp=30)
c
       dimension f(nx2,ny2,nz2)
c      dimension gf(nx2,ny2,nz2)
       dimension ppin(10)
!hpf$  distribute f(*,*,block) onto pe
c      equivalence (gf,f)
       common /blk/f
c
c
       xl=ppin(1)
       yl=ppin(2)
       zl=ppin(3)
       dxl=ppin(4)
       dyl=ppin(5)
       dzl=ppin(6)
       dn=ppin(7)
       dv=ppin(8)
c
       hx=xl/float(nx1)
       hy=yl/float(ny1)
       hz=zl/float(nz1)
c
!hpf$  independent,new(i,j,k)
       do 10 k=1,nz2
!hpf$  on home(f(:,:,k)),local(f,i,j,x,y,z,hx,hy,hz,
!hpf$*         dn1,dv1,ax1,dxl,dn) begin
       do 101 j=1,ny2
       do 101 i=1,nx2
       z=0.5*hz*(2*k-nz2-1)
       y=0.5*hy*(2*j-ny2-1)
       x=0.5*hx*(2*i-nx2-1)
       dn1=0.0
       dv1=0.0
       ax1=sqrt(x*x+y*y+z*z)
       if(ax1.le.dxl) dn1=dn
       f(i,j,k)=dn1
  101  continue
!hpf$  end on
   10  continue
c
       return
       end