c
       program pwave2
c
cc     pwave2.f  main program of 2-dimensional wave equtaion
c      modified leap-frog method  1991.12.14
       parameter(npe=8)
!xocl processor pe(npe)
       parameter(nx=100,ny=100)
c      parameter(nb=3,iip0=8,iiq0=4,iir0=4,last=4)
       parameter(nb=3,iip0=8,iiq0=1,iir0=1,last=4)
       parameter(nx1=nx+1,nx2=nx+2,ny1=ny+1,ny2=ny+2)
       parameter(n1=nx2,n2=n1*ny2,n3=n2*nb,noinp=30)
       parameter(thx=0.5,tam=0.5,vis0=0.10)
!xocl index partition ind=(pe,index=1:ny2,part=band)
!xocl index partition indo=(pe,index=1:ny2,part=band,overlap=(1,1))
c
       real*8 f(nx2,ny2,nb),u(nx2,ny2,nb)
       real*8 gf(nx2,ny2,nb),gu(nx2,ny2,nb)
       real*8 v(nx2,ny2,nb),p(nb)
       real*8 zt0,zt1,zt2,zt
       dimension ppin(10)
!xocl local f(:,/indo,:),u(:,/indo,:),v(:,/indo,:)
!xocl global gf,gu
       equivalence (gf,f),(gu,u)
       common /blk/gf
c
c           ppin xl    yl    dxl  dyl  dn  dv   7   8   9  10
       data ppin/101.0,101.0,10.0,10.0,1.0,0.0, 7,  8,  9, 10/
c
!xocl parallel region
c
       xl=ppin(1)
       yl=ppin(2)
       dxl=ppin(3)
       dyl=ppin(4)
       dn=ppin(5)
       dv=ppin(6)
       vis=vis0
c
       hx=xl/float(nx1)
       hy=yl/float(ny1)
       t1=thx*hx
       dx1=0.5*t1/hx
       dy1=0.5*t1/hy
       dx2=vis*0.50*(t1/hx)**2
       dy2=vis*0.50*(t1/hy)**2
c
c      initial parameters
c
       write(6,122) nx,ny,last,nx2,ny2
       write(6,124) t1,dx1,dy1,dx2,dy2
  122  format(1h ,10i8)
  124  format(1h ,8(1pe10.2))
c
c      initial conditions
       call ainit1(ppin)
c
c      start of calculation
       call clock(zt0)
c
       do 500 ii=1,last
       call clock(zt1)
       do 300 iir=1,iir0
       do 200 iiq=1,iiq0
       do 100 iip=1,iip0
!xocl overlapfix(f) (id)
!xocl movewait (id)
c
c      boundary conditions
c
       do 20 m=1,nb
!xocl spread do /ind
       do 22 j=2,ny1
       f(1,j,m)=f(nx1,j,m)
       f(nx2,j,m)=f(2,j,m)
   22  continue
!xocl end spread
!xocl spread move /indo,:
       do 24 j=1,1
       do 24 i=1,nx2
       f(i,j,m)=gf(i,j+ny,m)
   24  continue
!xocl end spread(id)
!xocl movewait(id)
!xocl spread move /indo,:
       do 26 j=ny2,ny2
       do 26 i=1,nx2
       f(i,j,m)=gf(i,j-ny,m)
   26  continue
!xocl end spread(id)
!xocl movewait(id)
   20  continue
c
!xocl overlapfix(f) (id)
!xocl movewait (id)
c
c      case of iip=1  2-step Lax-Wendroff method
       if(iip.eq.1) then
       t=0.5*t1
       dx=0.5*dx1
       dy=0.5*dy1
       do 30 m=1,nb
!xocl spread do /ind
       do 32 j=1,ny1
       do 32 i=1,nx1
       u(i,j,m)=0.25*(f(i,j,m)+f(i+1,j,m)+
     1                f(i,j+1,m)+f(i+1,j+1,m))
   32  continue
!xocl end spread
   30  continue
c
       else
       t=t1
       dx=dx1
       dy=dy1
       endif
c
!xocl overlapfix(f,u) (id)
!xocl movewait (id)
c
c      first step
c
!xocl spread do /ind
       do 40 j=1,ny1
       do 40 i=1,nx1
c
       p(1)=0.25*(f(i,j,1)+f(i+1,j,1)+
     1            f(i,j+1,1)+f(i+1,j+1,1))
       p(2)=0.25*(f(i,j,2)+f(i+1,j,2)+
     1            f(i,j+1,2)+f(i+1,j+1,2))
       p(3)=0.25*(f(i,j,3)+f(i+1,j,3)+
     1            f(i,j+1,3)+f(i+1,j+1,3))
c
       u(i,j,1)=u(i,j,1)
     1   -dx*(f(i+1,j+1,2)+f(i+1,j,2)-f(i,j+1,2)-f(i,j,2))
     2   -dy*(f(i+1,j+1,3)-f(i+1,j,3)+f(i,j+1,3)-f(i,j,3))
       u(i,j,2)=u(i,j,2)
     1   -dx*(f(i+1,j+1,1)+f(i+1,j,1)-f(i,j+1,1)-f(i,j,1))
       u(i,j,3)=u(i,j,3)
     1   -dy*(f(i+1,j+1,1)-f(i+1,j,1)+f(i,j+1,1)-f(i,j,1))
c
   40  continue
!xocl end spread
c
!xocl overlapfix(u) (id)
!xocl movewait (id)
c
c      preparation of second step
c      second step
c
!xocl spread do /ind
       do 62 j=1,ny2
c
       do 62 m=1,nb
       do 62 i=1,nx2
       v(i,j,m)=f(i,j,m)
   62  continue
!xocl end spread
c
!xocl overlapfix(v) (id)
!xocl movewait (id)
c
c
!xocl spread do /ind
       do 60 j=2,ny1
       do 60 i=2,nx1
c
       p(1)=0.25*(u(i,j,1)+u(i-1,j,1)+
     1            u(i,j-1,1)+u(i-1,j-1,1))
       p(2)=0.25*(u(i,j,2)+u(i-1,j,2)+
     1            u(i,j-1,2)+u(i-1,j-1,2))
       p(3)=0.25*(u(i,j,3)+u(i-1,j,3)+
     1            u(i,j-1,3)+u(i-1,j-1,3))
c
       f(i,j,1)=f(i,j,1)
     1   -dx1*(u(i,j,2)+u(i,j-1,2)-u(i-1,j,2)-u(i-1,j-1,2))
     2   -dy1*(u(i,j,3)-u(i,j-1,3)+u(i-1,j,3)-u(i-1,j-1,3))
     3   +dx2*(v(i-1,j,1)-2.0*v(i,j,1)+v(i+1,j,1))
     4   +dy2*(v(i,j+1,1)-2.0*v(i,j,1)+v(i,j-1,1))
       f(i,j,2)=f(i,j,2)
     1   -dx1*(u(i,j,1)+u(i,j-1,1)-u(i-1,j,1)-u(i-1,j-1,1))
     2   +dx2*(v(i-1,j,2)-2.0*v(i,j,2)+v(i+1,j,2))
     3   +dy2*(v(i,j+1,2)-2.0*v(i,j,2)+v(i,j-1,2))
       f(i,j,3)=f(i,j,3)
     1   -dy1*(u(i,j,1)-u(i,j-1,1)+u(i-1,j,1)-u(i-1,j-1,1))
     2   +dx2*(v(i-1,j,3)-2.0*v(i,j,3)+v(i+1,j,3))
     3   +dy2*(v(i,j+1,3)-2.0*v(i,j,3)+v(i,j-1,3))
c
   60  continue
!xocl end spread
c
c      end of 1 time step advance
c
  100  continue
  200  continue
  300  continue
c
       call clock(zt2)
       zt1=zt1-zt0
       zt2=zt2-zt0
       zt=zt2-zt1
       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,ti1
       ti=1.0d0
       ti1=1.0d0
       call gettod(ti1)
       ti=1.0d-6*ti1
c      x=0.0
c      y=secnds(x)
c      ti=1.0d0*y
       return
       end 
       subroutine ainit1(ppin)
       parameter(npe=8)
!xocl processor pe(npe)
!xocl subprocessor pes(npe)=pe(1:npe)
       parameter(nx=100,ny=100)
       parameter(nb=3,iip0=8,iiq0=4,iir0=4,last=4)
       parameter(nx1=nx+1,nx2=nx+2,ny1=ny+1,ny2=ny+2)
       parameter(n1=nx2,n2=n1*ny2,n3=n2*nb,noinp=30)
!xocl index partition ind=(pes,index=1:ny2,part=band)
!xocl index partition indo=(pes,index=1:ny2,part=band,overlap=(1,1))
c
       real*8 f(nx2,ny2,nb)
       real*8 gf(nx2,ny2,nb)
       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)
       dxl=ppin(3)
       dyl=ppin(4)
       dn=ppin(5)
       dv=ppin(6)
c
       hx=xl/float(nx1)
       hy=yl/float(ny1)
c
!xocl spread do /ind
       do 10 j=1,ny2
       do 10 i=1,nx2
       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)
       if(ax1.le.dxl) dn1=dn
       f(i,j,1)=dn1
       f(i,j,2)=0.0
       f(i,j,3)=0.0
   10  continue
!xocl end spread
c
       return
       end