       Program main
*     nx, ny, nz:		grid sizes
*     x, y, z:			grid arrays
*     gn, go:			`new' and `old' field arrays
*     gy, gyy, gzz:		angular derivatives
*     l2n, l2o:			`new' and `old' laplacian on the sphere
 
       integer nx, ny, nz
       parameter (nx = 128, ny = 32, nz = 32)
       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 l2n(0:ny,0:nz), l2o(0:ny,0:nz)
       real timebegin, timeend, tarry(2)
 
*     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
 
       timebegin = dtime(tarry)

       open (unit = 11, file = 'plot.gnu', status = 'unknown')
 
       call grid(x, y, z, nx, ny, nz, dx, dy, dz, du)
 
       call getparam(u, uf, du, iterskip, islice, isliceno)
 
       call surfdata(gn, x, y, z, nx, ny, nz, u)
 
       call wslice(gn, x, y, z, nx, ny, nz, islice, isliceno, u)
 
       iter = 0
       do while (u .lt. uf)
          iter = iter + 1
          u = u + du
          call copy(gn, go, nx, ny, nz)
          call boundary(gn, go, nx, ny, nz, dx, dy, dz, du, u)
          call advance(gn, go, gy, gyy, gzz, l2o, l2n, x, y, z, 
     &                   nx, ny, nz, dx, dy, dz, du)
 
          if (mod(iter,iterskip) .eq. 0) 
     &      call wslice(gn, x, y, z, nx, ny, nz, islice, isliceno, u)
 
       end do
 
       close (unit = 11)
 
       timeend = dtime(tarry)

       print *, 'Time Taken = ', timeend- timebegin

       end

       subroutine grid(x, y, z, nx, ny, nz, dx, dy, dz, du)
       real x(0:nx), y(0:ny), z(0:nz)
       integer nx, ny, nz
       real dx, dy, dz, du
 
       real pi
       integer i, j, k
 
       pi = 3.14159265358979323846
 
       dx = 0.5 / float(nx)
       dy = 2. / float(ny)
       dz = 2. * pi / float(nz + 1)
 
       du = min(dx,dy,dz) / 2.
 
       do i = 0, nx
           x(i) = 0.5 + i * dx
       end do
 
       do j = 0, ny
          y(j) = -1. + j * dy
       end do
 
       do k = 0, nz
           z(k) = k * dz
       end do
  
       return
       end

       subroutine surfdata(g, x, y, z, nx, ny, nz, u)
       integer nx, ny, nz
       real x(0:nx), y(0:ny), z(0:nz)
       real g(0:nx,0:ny,0:nz)
       real u
 
       integer i, j, k
       real x1, x2
 
       x1 = 0.6
       x2 = 0.9
 
       do k = 0, nz
          do j = 0, ny
             do i = 0, nx
                if ((x(i) .gt. x1) .and. (x(i) .lt. x2)) then
                   g(i,j,k) = ((x(i) - x1) * (x2 - x(i))) ** 2
     &                     * (1. - y(j) ** 2) 
     &                     * cos(2. * z(k))
                else
                   g(i,j,k) = 0.
                end if
             end do
          end do
       end do
 
       return
       end

       subroutine boundary(gn, go, nx, ny, nz, dx, dy, dz, du, u)
       real gn(0:nx,0:ny,0:nz), go(0:nx,0:ny,0:nz)
       integer nx, ny, nz
       real dx, dy, dz, du, u
 
*     ... for the time being, set the boundary to zero,
*         i.e. equivalent to having a perfect reflector at x=0.5
 
       integer j, k
 
       do k = 0, nz
          do j = 0, ny
             gn(0,j,k) = 0.
CCCCCC       gn(nx/2,j,k) = 0. this line is modified with Raberto
          end do
       end do
 
       return
       end

       subroutine wslice(g, x, y, z, nx, ny, nz, islice, isliceno, u)
       integer nx, ny, nz, islice, isliceno
       real x(0:nx), y(0:ny), z(0:nz)
       real g(0:nx,0:ny,0:nz)
       real u
 
       integer i, j, k
       character filename*14
 
       filename(1:2) = 'u='
       write (filename(3:14),'(e12.6)') u
       open (unit = 10, file = filename, status = 'unknown')
 
       if (islice .eq. 1) then
          i = isliceno
          do k = 0, nz
             write (10,120) (y(j), z(k), g(i,j,k), j = 0, ny)
             write (10,130)
          end do
       else if (islice .eq. 2) then
          j = isliceno
          do k = 0, nz
             write (10,120) (x(i), z(k), g(i,j,k), i = 0, nx)
             write (10,130)
          end do
       else if (islice .eq. 3) then
          k = isliceno
          do j = 0, ny
             write (10,120) (x(i), y(j), g(i,j,k), i = 0, nx)
             write (10,130)
          end do
       end if
 
       close (unit = 10)
 
       write (unit = 11, fmt=150) 'splot "' // filename // '"'
 
120     format(3(5x,e12.5))
130     format(' ')
150     format(a)
 
       return
       end

       subroutine copy(gn, go, nx, ny, nz)
       integer nx, ny, nz
       real gn(0:nx,0:ny,0:nz), go(0:nx,0:ny,0:nz)
 
       integer i, j, k
 
       do k = 0, nz
          do j = 0, ny
             do i = 0, nx
                go(i,j,k) = gn(i,j,k)
             end do
          end do
       end do
 
       return
       end

       subroutine advance(gn, go, gy, gyy, gzz, l2o, l2n, x, y, z, 
     &                   nx, ny, nz, dx, dy, dz, du)
       integer nx, ny, 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 l2o(0:ny,0:nz), l2n(0:ny,0:nz)
       real x(0:nx), y(0:ny), z(0:nz)
       real dx, dy, dz, du
 
       integer i, j, k
       real gx, gxx, xh
       real alpha, beta, gamma
*
*     ... 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
       do k = 0, nz
          do j = 0, ny
             gx = (go(1,j,k) - go(0,j,k)) / dx
             gxx = (3. * go(0,j,k) 
     &           - 7. * go(1,j,k) 
     &           + 5. * go(2,j,k)
     &               -  go(3,j,k)
     &            ) / (2. * dx * dx)
             gn(1,j,k) = gn(0,j,k)
     &                + go(1,j,k)
     &                - go(0,j,k)
     &                + ((1. - xh) ** 2 * gxx - 2. * (1. - xh) * gx
     &                   - 0.5 * (l2o(j,k) + l2n(j,k)) / (xh * xh)
     &                  ) * (0.5 * du * dx)
          end do
       end do
*
*     ... over the bulk of the grid ...
*
       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 ...
*
          do k = 0, nz
             do j = 0, ny
                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
          end do
 
 
       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 ...
*
       do k = 0, nz
          do j = 0, ny
                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.))
          end do
       end do
 
       return
       end

       subroutine mkl2(g, gy, gyy, gzz, l2, y, nx, ny, nz, dy, dz, i)
       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
 
       integer j, k
       real g_north,  g_south
       real gm_north, gm_south 
 
*
*        calculate the `phi' derivatives
*
 
          do j = 0, ny
             gzz(j,0) = (g(i,j,1)   
     &           - 2. * g(i,j,0)
     &                + g(i,j,nz)
     &                  ) / (dz * dz)
             do k = 1, nz - 1
                gzz(j,k) = (g(i,j+1,k)   
     &              - 2. * g(i,j,k)   
     &                   + g(i,j-1,k)
     &                    ) / (dz * dz)
             end do
             gzz(j,nz) = (g(i,j,0)   
     &            - 2. * g(i,j,nz)
     &                 + g(i,j,nz-1)
     &                  ) / (dz * dz)
          end do
 
*
*        ... estimate the average value of g on a small circle around
*            the poles
*
          g_north  = 0.
          g_south  = 0.
          gm_north = 0.
          gm_south = 0.
          do k = 0, nz
             g_north  = g_north  + g(i,ny,k)
             g_south  = g_south  + g(i,0,k)
             gm_north = gm_north + g(i,ny-1,k)
             gm_south = gm_south + g(i,1,k)
          end do
          g_north  = g_north / float(nz + 1)
          g_south  = g_south / float(nz + 1)
          gm_north = gm_north / float(nz + 1)
          gm_south = gm_south / float(nz + 1)
 
*
*        ... now the `theta' derivatives
*
 
          do k = 0, nz
             gy(0,k) = (- 3. * g(i,0,k) 
     &                 + 4. * g(i,1,k)
     &                      - g(i,2,k)
     &                ) / (2. * dy)
             gyy(0,k) = (2. * g(i,0,k) 
     &                - 5. * g(i,1,k)
     &                + 4. * g(i,2,k)
     &                     - g(i,3,k)
     &                 ) / (dy * dy)
             do j = 1, ny - 1
                gy(j,k) = (g(i,j+1,k)
     &                  - g(i,j-1,k)
     &                   ) / (2. * dy)
                gyy(j,k) = (g(i,j+1,k)
     &              - 2. * g(i,j,k)   
     &                   + g(i,j-1,k)
     &                    ) / (dy * dy)
             end do
             gy(ny,k) = (3. * g(i,ny,k) 
     &                - 4. * g(i,ny-1,k)
     &                     + g(i,ny-2,k)
     &                 ) / (2. * dy)
             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)
          end do
 
*
*        ... assemble the L^2 operator ...
*
 
          do k = 0, nz
             l2(0,k) = - 4. * (gm_south - g_south) / (dz * dz)
             do j = 1, ny - 1
                l2(j,k) = - ((1. - y(j) ** 2) * gyy(j,k) 
     &                      - 2. * y(j) * gy(j,k) 
     &                      + gzz(j,k) / (1. - y(j) * y(j))) 
             end do
             l2(ny,k) = - 4. * (gm_north - g_north) / (dz * dz)
          end do
 
       return
       end

       subroutine getparam(ubegin, ufinal, du, iterskip, islice, 
     &                    isliceno)
       real ubegin, ufinal, du
       integer iterskip
       integer islice, isliceno
 
       write (*,'(a$)') 'u initial: '
       read (*,*) ubegin
       write (*,'(a$)') 'u final: '
       read (*,*) ufinal
       write (*,'(a,i4)') 'estimated iters: ', int(ufinal / du)
       write (*,'(a$)') 'iterskip: '
       read (*,*) iterskip
       write (*,'(a$)') 'slicing choice [1=x, 2=y, 3=z]: '
       read (*,*) islice
       write (*,'(a$)') 'slice #: '
       read (*,*) isliceno
 
       return
       end
