     1  c
     2         program pwave3
     3  c
     4  cc     pwave3.f  main program of 3-dimensional wave equtaion
     5  c      modified leap-frog method  1991.12.14
     6         parameter(npe=2)
     7  !xocl processor pe(npe)
     8         parameter(nx=100,ny=100,nz=100)
     9         parameter(nb=4,iip0=8,iiq0=1,iir0=1,last=4)
    10         parameter(nx1=nx+1,nx2=nx+2,ny1=ny+1,ny2=ny+2)
    11         parameter(nz1=nz+1,nz2=nz+2)
    12         parameter(n1=nx2,n2=n1*ny2,n3=n2*nz2,noinp=30)
    13         parameter(n4=n3*nb)
    14         parameter(thx=0.5,tam=0.5,vis0=0.10)
    15  !xocl index partition ind=(pe,index=1:nz2,part=band)
    16  !xocl index partition indo=(pe,index=1:nz2,part=band,overlap=(1,1))
    17  c
    18         real*8 f(nx2,ny2,nz2,nb),u(nx2,ny2,nz2,nb)
    19         real*8 gf(nx2,ny2,nz2,nb),gu(nx2,ny2,nz2,nb)
    20         real*8 v(nx2,ny2,nz2,nb),p(nb)
    21         real*8 gv(nx2,ny2,nz2,nb)
    22         dimension ppin(10)
    23         real*8 zt0,zt1,zt2,zt
    24  !xocl local f(:,:,/indo,:),u(:,:,/indo,:),v(:,:,/indo,:)
    25  !xocl global gf,gu,gv
    26         equivalence (gf,f),(gu,u),(gv,v)
    27         common /blk/gf
    28  c
    29  c           ppin xl   yl   zl   dxl  dyl  dzl  dn  dv   9  10
    30         data ppin/62.0,62.0,62.0,10.0,10.0,10.0,1.0,0.0, 9, 10/
    31  c
    32  !xocl parallel region
    33  c

    55  c      initial parameters
    56  c
    57         write(6,122) nx,ny,nz,last,nx2,ny2,nz2
    58         write(6,124) t1,hx,hy,hz,dx1,dy1,dz1
    59    122  format(1h ,10i8)
    60    124  format(1h ,8(1pe10.3))
    61  c
    62  c      initial conditions
    63         call ainit1(ppin)
    64  c
    65  c      start of calculation
    66         call clock(zt0)
    67  c
    68         do 500 ii=1,last
    69         call clock(zt1)
    70         do 300 iir=1,iir0
    71         do 200 iiq=1,iiq0
    72         do 100 iip=1,iip0
    73  !xocl overlapfix(f) (id)
    74  !xocl movewait (id)
    75  c
    76  c      boundary conditions
    77  c
    78         do 20 m=1,nb
    79  !xocl spread do /ind
    80         do 24 k=2,nz1
    81         do 24 j=2,ny1
    82         f(1,j,k,m)=f(nx1,j,k,m)
    83         f(nx2,j,k,m)=f(2,j,k,m)
    84     24 continue
    85  !xocl end spread
    86  !xocl spread do /ind
    87         do 22 k=2,nz1
    88         do 22 i=1,nx2
    89         f(i,1,k,m)=f(i,ny1,k,m)
    90         f(i,ny2,k,m)=f(i,2,k,m)
    91     22  continue
    92  !xocl end spread
    93  !xocl spread move /indo,:,:
    94         do 26 k=1,1
    95         do 26 j=1,ny2
    96         do 26 i=1,nx2
    97  c      f(i,j,1,m)=gf(i,j,nz1,m)
    98  c      f(i,j,nz2,m)=gf(i,j,2,m)
    99         f(i,j,k,m)=gf(i,j,k+nz,m)
   100     26  continue
   101  !xocl end spread(id)
   102  !xocl movewait(id)
   103  !xocl spread move /indo,:,:
   104         do 28 k=nz2,nz2
   105         do 28 j=1,ny2
   106         do 28 i=1,nx2
   107  c      f(i,j,1,m)=gf(i,j,nz1,m)
   108  c      f(i,j,nz2,m)=gf(i,j,2,m)
   109         f(i,j,k,m)=gf(i,j,k-nz,m)
   110     28  continue
   111  !xocl end spread(id)
   112  !xocl movewait(id)
   113     20  continue
   114  c
   115  !xocl overlapfix(f) (id)
   116  !xocl movewait (id)
   117  c
   118  c      case of iip=1  2-step Lax-Wendroff method
   119         if(iip.eq.1) then
   120         t=0.5*t1
   121         dx=0.5*dx1
   122         dy=0.5*dy1
   123         dz=0.5*dz1
   124  !xocl spread do /ind
   125         do 32 k=1,nz1
   126         do 32 m=1,nb
   127         do 32 j=1,ny1
   128         do 32 i=1,nx1
   129         u(i,j,k,m)=0.125*(f(i,j,k,m)+f(i+1,j,k,m)+

   133     32  continue
   134  !xocl end spread
   135  c
   136         else
   137         t=t1
   138         dx=dx1
   139         dy=dy1
   140         dz=dz1
   141         endif
   142  c
   143  !xocl overlapfix(f,u) (id)
   144  !xocl movewait (id)
   145  c
   146  c      first step
   147  c
   148  !xocl spread do /ind
   149         do 40 k=1,nz1
   150         do 40 j=1,ny1
   151         do 40 i=1,nx1
   152  c
   153         p(1)=0.125*(f(i,j,k,1)+f(i+1,j,k,1)+

   199     40  continue
   200  !xocl end spread
   201  c
   202  !xocl overlapfix(u) (id)
   203  !xocl movewait (id)
   204  c
   205  c      preparation of second step
   206  c      second step
   207  c
   208  !xocl spread do /ind
   209         do 62 k=1,nz2
   210         do 62 m=1,nb
   211         do 62 j=1,ny2
   212         do 62 i=1,nx2
   213         v(i,j,k,m)=f(i,j,k,m)
   214     62  continue
   215  !xocl end spread
   216  c
   217  !xocl overlapfix(v) (id)
   218  !xocl movewait (id)
   219  c
   220  !xocl spread do /ind
   221         do 60 k=2,nz1
   222         do 60 j=2,ny1
   223         do 60 i=2,nx1
   224  c
   225         p(1)=0.125*(u(i,j,k,1)+u(i-1,j,k,1)+

   283     60  continue
   284  !xocl end spread
   285  c
   286  c      end of 1 time step advance
   287  c
   288    100  continue
   289    200  continue
   290    300  continue
   291  c
   292         call clock(zt2)
   293         zt1=zt1-zt0
   294         zt2=zt2-zt0
   295         zt=zt2-zt1
   296         write(6,402) ii,zt0,zt1,zt2,zt
   297    402  format(1h , i6,1pe12.3,3(0pf12.5))
   298  c
   299  c      write the output data
   300  c
   301    500  continue
   302      9  continue
   303  !xocl end parallel
   304  c
   305         stop
   306         end
   307         subroutine clock(ti)
   308         real*8 ti,ti1
   309         ti=1.0d0
   310         ti1=1.0d0
   311         call gettod(ti1)
   312         ti=1.0d-6*ti1
   313  c      x=0.0
   314  c      y=secnds(x)
   315  c      ti=1.0d0*y
   316         return
   317         end
   318         subroutine ainit1(ppin)
   319         parameter(npe=2)
   320  !xocl processor pe(npe)
   321  !xocl subprocessor pes(npe)=pe(1:npe)
   322         parameter(nx=100,ny=100,nz=100)
   323         parameter(nb=4,iip0=8,iiq0=1,iir0=1,last=4)
   324         parameter(nx1=nx+1,nx2=nx+2,ny1=ny+1,ny2=ny+2)
   325         parameter(nz1=nz+1,nz2=nz+2)
   326         parameter(n1=nx2,n2=n1*ny2,n3=n2*nz2,noinp=30)
   327         parameter(n4=n3*nb)
   328  !xocl index partition ind=(pes,index=1:nz2,part=band)
   329  !xocl index partition indo=(pes,index=1:nz2,part=band,overlap=(1,1))
   330  c
   331         real*8 f(nx2,ny2,nz2,nb)
   332         real*8 gf(nx2,ny2,nz2,nb)
   333         dimension ppin(10)
   334  !xocl local f(:,:,/indo,:)
   335  !xocl global gf
   336         equivalence (gf,f)
   337         common /blk/gf
   338  c
   339  !xocl overlapfix(f) (id)
   340  !xocl movewait (id)
   341  c

   354  !xocl spread do /ind
   355         do 10 k=1,nz2
   356         z=0.5*hz*(2*k-nz2-1)
   357         do 10 j=1,ny2
   358         y=0.5*hy*(2*j-ny2-1)
   359         do 10 i=1,nx2

   369     10  continue
   370  !xocl end spread
   371         return
   372         end

