program t
real*8,dimension(:,:,:,:),allocatable :: temp
integer                               :: nx,ny,nz
real*8,dimension(:),allocatable       :: x,y,z
real*8                                :: h

interface 
   subroutine makegrid(n,h,x,y,z)
      implicit none
      real*8                       :: h
      real*8,dimension(:)          :: x,y,z
      integer                      :: i,n
   end subroutine

   subroutine makegrid3d(n,h,xx,yy,zz,x,y,z,r)
      implicit none
      real*8                       :: h
      real*8,dimension(:)          :: xx,yy,zz
      real*8,dimension(:,:,:)      :: x,y,z,r
      integer                      :: i,n
	end subroutine
end interface


	write(6,*) 'input n'
	read(5,*) nx
	ny = nx
	nz = nx
   h=  2.0d0 / (nx -1 )

	allocate(x(nx),y(ny),z(nz))
	allocate(temp(4,nx,ny,nz))

	call makegrid(nx,h,x,y,z)

	call makegrid3d(nx,h,x,y,z,temp(1,:,:,:),temp(2,:,:,:),    &
                   temp(3,:,:,:),temp(4,:,:,:) )
	
end



subroutine makegrid(n,h,x,y,z)
implicit none
real*8                       :: h
real*8,dimension(:)          :: x,y,z
integer                      :: i,n


   do i = 1, n 
      x(i) = -1.0d0 + (i-1.0d0)*h
   end do
   do i = 1, n 
      y(i) = -1.0d0 + (i-1.0d0)*h
   end do
   do i = 1, n 
      z(i) = -1.0d0 + (i-1.0d0)*h
   end do

return
end


subroutine makegrid3d(n,h,xx,yy,zz,x,y,z,r)
implicit none
real*8                       :: h
real*8,dimension(:)          :: xx,yy,zz
real*8,dimension(:,:,:)      :: x,y,z,r 
integer                      :: i,n

  do i=1,n 
     x(i,:,:)=xx(1)+(i-1)*h
  enddo
  do i=1,n 
     y(:,i,:)=yy(1)+(i-1)*h
  enddo
  do i=1,n 
     z(:,:,i)=zz(1)+(i-1)*h
  enddo

	r = sqrt( x**2 + y**2 + z**2 )


	do i = 1 , 17
		write(6,*) r(i,9,9)
	end do


return
end
