c
       program adiff32
c  
c      main program of 3-dimensional difussion equation
c      modified leap-frog method  2000.06.10
       parameter(npe=16)
!xocl processor 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)
!xocl index partition ind=(pe,index=1:nz2,part=band)
!xocl index partition indo=(pe,index=1:nz2,part=band,overlap=(1,1))
c
       real*8 f(nx2,ny2,nz2),u(nx2,ny2,nz2)
       real*8 gf(nx2,ny2,nz2),gu(nx2,ny2,nz2)
       real*8 rinp(noinp)
       dimension ppin(10)
c      integer iinp(noinp)
       real*8 zt0,zt1,zt2,zt
!xocl local f(:,:,/indo),u(:,:,/indo)
!xocl global gf,gu
       equivalence (gf,f),(gu,u)
       common /blk/gf
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
!xocl parallel region
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
!xocl overlapfix(f) (id)
!xocl movewait (id)
c
c      boundary conditions
c
c      periodic boundary in x-direction
!xocl spread do /ind
       do 22 k=2,nz1
       do 22 j=2,ny1
       f(1,j,k)=f(nx1,j,k)
       f(nx2,j,k)=f(2,j,k)
   22 continue
!xocl end spread
c
c      periodic boundary in y-direction
!xocl spread do /ind
       do 24 k=2,nz1
       do 24 i=1,nx2
       f(i,1,k)=f(i,ny1,k)
       f(i,ny2,k)=f(i,2,k)
   24  continue
!xocl end spread
c
c      periodic boundary in z-direction
!xocl spread move /indo,:,:
       do 26 k=1,1
       do 26 j=1,ny2
       do 26 i=1,nx2
c      f(i,j,1)=f(i,j,nz1)
c      f(i,j,nz2)=f(i,j,2)
       f(i,j,k)=gf(i,j,k+nz)
   26  continue
!xocl end spread(id)
!xocl movewait(id)
c
!xocl spread move /indo,:,:
       do 28 k=nz2,nz2
       do 28 j=1,ny2
       do 28 i=1,nx2
c      f(i,j,1)=gf(i,j,nz1)
c      f(i,j,nz2)=gf(i,j,2)
       f(i,j,k)=gf(i,j,k-nz)
   28  continue
!xocl end spread(id)
!xocl movewait(id)
c
!xocl overlapfix(f) (id)
!xocl movewait (id)
c
c
!xocl spread do /ind
       do 60 k=1,nz2
       do 60 j=1,ny2
       do 60 i=1,nx2
       u(i,j,k)=f(i,j,k)
   60  continue
!xocl end spread
!xocl overlapfix(u) (id)
!xocl movewait (id)
c
c
!xocl spread do /ind
       do 80 k=2,nz1
       do 80 j=2,ny1
       do 80 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
   80  continue
!xocl end spread
!xocl overlapfix(f) (id)
!xocl movewait (id)
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
!xocl end parallel
c
       stop
       end
       subroutine clock(ti)
       real*8 ti
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)
       parameter(npe=16)
!xocl processor pe(npe)
!xocl subprocessor pes(npe)=pe(1: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)
!xocl index partition ind=(pes,index=1:nz2,part=band)
!xocl index partition indo=(pes,index=1:nz2,part=band,overlap=(1,1))
c
       real*8 f(nx2,ny2,nz2)
       real*8 gf(nx2,ny2,nz2)
       dimension ppin(10)
!xocl local f(:,:,/indo)
!xocl global gf
       equivalence (gf,f)
       common /blk/gf
c
!xocl overlapfix(f) (id)
!xocl movewait (id)
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
!xocl spread do /ind
       do 10 k=1,nz2
       do 10 j=1,ny2
       do 10 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
   10  continue
!xocl end spread
c
       return
       end