       module parameters
       integer, parameter :: PX=2,PY=2
       integer, parameter :: NP=PX*PY 
       end module parameters
         
       program wave_extract
       use parameters
       integer nx, ny, nz
       parameter (nx = 128, ny = 64, nz = 64)
c       parameter (nx = 8, ny = 8, nz = 8)
       real x(0:nx), y(0:ny), z(0:nz)
       real gn(0:nx,0:ny,0:nz), go(0:nx,0:ny,0:nz)
       real gy(0:ny,0:nz), gyy(0:ny,0:nz), gzz(0:ny,0:nz)
       real gx(0:ny,0:nz), gxx(0:ny,0:nz)
       real l2n(0:ny,0:nz), l2o(0:ny,0:nz)
 
*     dx, dy, dz:		grid spacings
*     u, uf, du:		initial time, final time, time step
*     iter, iterskip:		iteration counter, no. of iter between dumps
*     islice, isliceno:		selects x, y or z slicing, slice number
 
       real dx, dy, dz
       real u, uf, du
       integer iter, iterskip
       integer islice, isliceno
        

!HPF$ PROCESSORS PPP(PX,PY)
!HPF$ DISTRIBUTE GN(*,BLOCK,BLOCK) ONTO PPP
!HPF$ ALIGN (:,:,:) WITH GN(:,:,:) :: go
!HPF$ ALIGN (:,:)   WITH GN(*,:,:) :: gy, gyy, gzz, l2n, l2o, gx, gxx
!HPF$ ALIGN y(:) WITH GN(*,:,*)
!HPF$ ALIGN z(:) WITH GN(*,*,:)

          real t1(NP), t2(NP)
!HPF$   distribute t1(block)
!HPF$   distribute t2(block)


       INTERFACE
       subroutine mkl2(g, gy, gyy, gzz, l2, y, nx, ny, nz, dy, dz, i)
c      ----------------------------------------------------------
       use parameters
       integer nx, ny, nz
       real g(0:nx,0:ny,0:nz)
       real gy(0:ny,0:nz), gyy(0:ny,0:nz), gzz(0:ny,0:nz), l2(0:ny,0:nz)
       real y(0:ny)
       real dy, dz
       integer i
 
!HPF$ PROCESSORS PPP(PX,PY)
!HPF$ DISTRIBUTE G *(*,BLOCK,BLOCK) ONTO *PPP
!HPF$ ALIGN (:,:)   WITH *G(*,:,:) :: gy, gyy, gzz, l2
!HPF$ ALIGN y(:) WITH *G(*,:,*)
      end subroutine mkl2

       extrinsic(hpf_local) subroutine timer(return_time, initial_time)
       real, dimension(:), intent(in)  :: initial_time
       real, dimension(:), intent(out) :: return_time
!HPF$ distribute *(block) :: initial_time
!HPF$ align return_time with initial_time
       end subroutine timer


      END INTERFACE 

!          start_timer
       t1 = 0.0
       call timer(t2, t1)


 
ccc       open (unit = 11, file = 'plot.gnu', status = 'unknown')
c      --------------------------------------------------------------- 
c       call grid(x, y, z, nx, ny, nz, dx, dy, dz, du)
c      ----------------------------------------------------------------
       pi = 3.14159265358979323846
 
       dx = 0.5 / float(nx)
       dy = 2. / float(ny)
       dz = 2. * pi / float(nz + 1)
 
       du = min(dx,dy,dz) / 2.
 
       FORALL(i=0:nx) x(i) = 0.5 + i * dx
       FORALL(j=0:ny) y(j) = -1. + j * dy
       FORALL(k=0:nz) z(k) = k * dz

c      ----------------------------------------------------------------- 
c      call getparam(u, uf, du, iterskip, islice, isliceno)
c      -----------------------------------------------------------------
       write (6,*) 'u initial: 0' 
       u = 0
       write (6,*) 'u final: 1' 
       uf = 1
       write (6,*) 'estimated iters: ', int(uf / du)
 
c      -----------------------------------------------------------------
c       call surfdata(gn, x, y, z, nx, ny, nz, u)
c      -----------------------------------------------------------------
       x1 = 0.6
       x2 = 0.9

       gn = 0
       FORALL(i=0:nx, j=0:ny, k=0:nz, (x(i).gt.x1) .and. (x(i).lt.x2)) 
     &             gn(i,j,k) = ((x(i) - x1) * (x2 - x(i))) ** 2
     &                     * (1. - y(j) ** 2) 
     &                     * cos(2. * z(k))

c      -----------------------------------------------------------------
ccc       call wslice(gn, x, y, z, nx, ny, nz, islice, isliceno, u)
c      -----------------------------------------------------------------
 
       iter = 0
       do while (u .lt. uf)
          iter = iter + 1
          u = u + du
c      -----------------------------------------------------------------
c          call copy(gn, go, nx, ny, nz)
c      -----------------------------------------------------------------
           FORALL(i=0:nx, j=0:ny, k=0:nz) go(i,j,k) = gn(i,j,k)
c      -----------------------------------------------------------------
c          call boundary(gn, go, nx, ny, nz, dx, dy, dz, du, u)
c      -----------------------------------------------------------------
           FORALL(j=0:ny,k=0:nz) gn(0,j,k) = 0.

c      -----------------------------------------------------------------
c          call advance(gn, go, gy, gyy, gzz, l2o, l2n, x, y, z, 
c     &                   nx, ny, nz, dx, dy, dz, du)
c      -----------------------------------------------------------------
c      ..........................................................
*
*     ... at the inner boundary
*         evaluate the L^2 operator acting on g ...
*         on the new and old levels ...
*
          call mkl2(go, gy, gyy, gzz, l2o, y, nx, ny, nz, dy, dz, 1)
          call mkl2(gn, gy, gyy, gzz, l2n, y, nx, ny, nz, dy, dz, 0)
*
*     ... and insert into a modified form of the evolution algorithm ...
*
       xh = x(0) + 0.5 * dx

       FORALL(j=0:ny, k=0:nz) gx(j,k) = (go(1,j,k) - go(0,j,k)) / dx
       FORALL(j=0:ny, k=0:nz) gxx(j,k) = (3. * go(0,j,k) 
     &           - 7. * go(1,j,k) 
     &           + 5. * go(2,j,k)
     &               -  go(3,j,k)
     &            ) / (2. * dx * dx)

       FORALL(j=0:ny, k=0:nz) gn(1,j,k) = gn(0,j,k)
     &       + go(1,j,k)
     &         - go(0,j,k)
     &        + ((1. - xh) ** 2 * gxx(j,k) - 2. * (1. - xh) * gx(j,k)
     &                   - 0.5 * (l2o(j,k) + l2n(j,k)) / (xh * xh)
     &                  ) * (0.5 * du * dx)
*
*     ... over the bulk of the grid ...
*
c       stop
       do i = 2, nx - 1
 
          alpha = (1. - (x(i) - 0.5 * dx)) / dx
          beta  = 4. / (du * dx)
          gamma = - 2. / (x(i) - 0.5 * dx) ** 2

*
*        ... evaluate the L^2 operator acting on g
*            on the new and old levels ...
*
          call mkl2(go, gy, gyy, gzz, l2o, y, nx, ny, nz, dy, dz, i)
          call mkl2(gn, gy, gyy, gzz, l2n, y, nx, ny, nz, dy, dz, i-1)

*
*        ... and insert into the evolution algorithm ...
*
      
         FORALL(j=0:ny, k=0:nz) 
     &          gn(i,j,k) = (gamma * 0.5 * (l2o(j,k) + l2n(j,k))
     &         + (beta - 2 * alpha * alpha)      * gn(i-1,j,k)
     &         + alpha * (alpha + 1.)            * gn(i-2,j,k)
     &         + alpha * (alpha - 1.)            * go(i+1,j,k)
     &         + (beta - 2 * alpha * alpha)      * go(i,j,k)
     &         + (- beta + alpha * (alpha + 1.)) * go(i-1,j,k)
     &         ) / (beta - alpha * (alpha - 1.))
 
       end do
*
*     ... at the outer boundary
*         evaluate the L^2 operator acting on g
*         on the old level ...
*
       call mkl2(go, gy, gyy, gzz, l2o, y, nx, ny, nz, dy, dz, nx)
       call mkl2(gn, gy, gyy, gzz, l2n, y, nx, ny, nz, dy, dz, nx-1)
 
*
*     ... and insert into a modified form of the evolution algorithm ...
*
         FORALL(j=0:ny, k=0:nz) 
     &          gn(nx,j,k) = (gamma * 0.5 * (l2o(j,k) + l2n(j,k))
     &         + (beta - 2 * alpha * alpha)        * gn(nx-1,j,k)
     &         + alpha * (alpha + 1.)              * gn(nx-2,j,k)
     &         + (beta + alpha * (2. * alpha - 3.)) * go(nx,j,k)
     &         + (-beta - alpha * (5. * alpha - 4.)) * go(nx-1,j,k)
     &         + alpha * (4. * alpha - 1.) * go(nx-2,j,k)
     &         + (-alpha * alpha) * go(nx-3,j,k)
     &         ) / (beta - alpha * (alpha - 1.))
 
 
ccc          if (mod(iter,iterskip) .eq. 0) 
ccc     &      call wslice(gn, x, y, z, nx, ny, nz, islice, isliceno, u)
 
       end do

!         stop_timer
          call timer(t1, t2)
          print *,'time',t1,' seconds'
  

 
ccc       close (unit = 11)

       write(6,*) "gn(0,      0,   0)=", gn(0,0,0)
       write(6,*) "gn(1,      1,   1)=", gn(1,1,1)
       write(6,*) "gn(2,      2,   2)=", gn(2,2,2)
       write(6,*) "gn(1,   ny-5,   2)=", gn(1,ny-5,2)
       write(6,*) "gn(nx-2,   0,   0)=", gn(nx-2,0,0)
       write(6,*) "gn(nx-2,ny-5,nz-5)=", gn(nx-2,ny-5,nz-5)
       write(6,*) "gn(nx,  ny-1,nz-1)=", gn(nx,ny-1,nz-1)

 
       end




c      ----------------------------------------------------------
       subroutine mkl2(g, gy, gyy, gzz, l2, y, nx, ny, nz, dy, dz, i)
c      ----------------------------------------------------------
       use parameters
       integer nx, ny, nz
       real g(0:nx,0:ny,0:nz)
       real gy(0:ny,0:nz), gyy(0:ny,0:nz), gzz(0:ny,0:nz), l2(0:ny,0:nz)
       real y(0:ny)
       real dy, dz
       integer i
 

!HPF$ PROCESSORS PPP(PX,PY)
!HPF$ DISTRIBUTE G *(*,BLOCK,BLOCK) ONTO *PPP
!HPF$ ALIGN (:,:)   WITH *G(*,:,:) :: gy, gyy, gzz, l2
!HPF$ ALIGN y(:) WITH *G(*,:,*)
  
       integer j, k
       real g_north,  g_south
       real gm_north, gm_south 

c      ..........................................................
*
*        calculate the `phi' derivatives
*
       FORALL(j=0:ny) gzz(j,0) = (g(i,j,1)   
     &           - 2. * g(i,j,0)
     &                + g(i,j,nz)
     &                  ) / (dz * dz)

       FORALL(j=0:ny, k=1:nz-1) gzz(j,k) = (g(i,j,k+1)   
     &              - 2. * g(i,j,k)   
     &                   + g(i,j,k-1)
     &                    ) / (dz * dz)
       FORALL(j=0:ny) gzz(j,nz) = (g(i,j,0)   
     &            - 2. * g(i,j,nz)
     &                 + g(i,j,nz-1)
     &                  ) / (dz * dz)
 
*
*        ... estimate the average value of g on a small circle around
*            the poles
*
          g_north  = SUM(g(i, ny,   0:nz)) / float(nz + 1)
          g_south  = SUM(g(i, 0,    0:nz)) / float(nz + 1)
          gm_north = SUM(g(i, ny-1, 0:nz)) / float(nz + 1)
          gm_south = SUM(g(i, 1,    0:nz)) / float(nz + 1)

*
*        ... now the `theta' derivatives
*
 

          FORALL(k=0:nz) gy(0,k) = (- 3. * g(i,0,k) 
     &                 + 4. * g(i,1,k)
     &                      - g(i,2,k)
     &                ) / (2. * dy)

          FORALL(k=0:nz) gyy(0,k) = (2. * g(i,0,k) 
     &                - 5. * g(i,1,k)
     &                + 4. * g(i,2,k)
     &                     - g(i,3,k)
     &                 ) / (dy * dy)

          FORALL(j=1:ny-1, k=0:nz)    gy(j,k) = (g(i,j+1,k)
     &                  - g(i,j-1,k)
     &                   ) / (2. * dy)
          FORALL(j=1:ny-1, k=0:nz)    gyy(j,k) = (g(i,j+1,k)
     &              - 2. * g(i,j,k)   
     &                   + g(i,j-1,k)
     &                    ) / (dy * dy)

          FORALL(k=0:nz) gy(ny,k) = (3. * g(i,ny,k) 
     &                - 4. * g(i,ny-1,k)
     &                     + g(i,ny-2,k)
     &                 ) / (2. * dy)
          FORALL(k=0:nz) gyy(ny,k) = (2. * g(i,ny,k) 
     &                 - 5. * g(i,ny-1,k)
     &                 + 4. * g(i,ny-2,k)
     &                      - g(i,ny-3,k)
     &                 ) / (dy * dy)

 
*
*        ... assemble the L^2 operator ...
*
 
          FORALL(k=0:nz)
     &        l2(0,k) = - 4. * (gm_south - g_south) / (dz * dz)

          FORALL(j=1: ny-1, k=0:nz)
     &          l2(j,k) = - ((1. - y(j) ** 2) * gyy(j,k) 
     &                      - 2. * y(j) * gy(j,k) 
     &                      + gzz(j,k) / (1. - y(j) * y(j))) 

          FORALL(k=0:nz)
     &       l2(ny,k) = - 4. * (gm_north - g_north) / (dz * dz)


       return
       end

       extrinsic(hpf_local) subroutine timer(return_time, initial_time)
       real, dimension(:), intent(in)  :: initial_time
       real, dimension(:), intent(out) :: return_time

       return_time(lbound(return_time, 1)) = 
     &       secnds(initial_time(lbound(initial_time, 1)))

       end

