!!*********************************************************************** 
! Front end routine to initial data: 
!   - calls appropriate subroutine 
!*********************************************************************** 
subroutine initial(x,y,z,params,t,g,k)
implicit none
real*8,dimension(:,:,:,:)    :: g,k,t
real*8,dimension(:)          :: x,y,z
real*8,dimension(:)          :: params

interface

   subroutine init_iso(xx,yy,zz,params,t,g,k)
      implicit none
      real*8,dimension(:,:,:,:)           :: g,k,t
      real*8,dimension(:)                 :: params
      real*8,dimension(:)                 :: xx,yy,zz
      integer                             :: i
   end subroutine

end interface

!params ..........
! (4) type of initial condition: 1 = isotropic Schwarzschild

	if (int(params(4)).eq.1) then
      call init_iso(x,y,z,params,t,g,k)
	else
		write(6,*) '*** error in initial data type ***'
		stop
	end if

return 
end 
!*********************************************************************** 
!  Routine to generate exact Isotropic schwarzschild initial data 
!  params(1) = nx 
!  params(2) = ny 
!  params(3) = nz 
!  params(4) = type of initial condition
!  params(5) = amplitude of linear wave
!  params(6) = width of linear wave
!  params(7) = radius of hole 1
!  params(8) = radius of hole 2
!  params(9) = center of hole 1 in x direction 
!  params(10) = center of hole 1 in y direction 
!  params(11) = center of hole 1 in z direction 
!  params(12) = center of hole 2 in x direction 
!  params(13) = center of hole 2 in y direction 
!  params(14) = center of hole 2 in z direction 
!   params(15) = xmin
!   params(16) = xmax
!   params(17) = ymin
!   params(18) = ymax
!   params(19) = zmin
!   params(20) = zmax
!   params(21) = gauge type
!   params(22) = tolerance
!   params(23) = nt

!*********************************************************************** 
subroutine init_iso(xx,yy,zz,params,t,g,k)
implicit none
real*8,dimension(:,:,:,:)           :: g,k,t
real*8,dimension(:)                 :: params
real*8,dimension(:)                 :: xx,yy,zz
integer                             :: i

real*8,dimension(2)                :: bha
real*8,dimension(3,2)              :: bhr
real*8                             :: test,h

	bha(1)       = params(7)
	bha(2)       = params(8)
	bhr(1,1)     = params(9)
	bhr(2,1)     = params(10)
	bhr(3,1)     = params(11)
	bhr(1,2)     = params(12)
	bhr(2,2)     = params(13)
	bhr(3,2)     = params(14)


!  if only one black hole, compute initial data as follows 
!        compute conformal factor for a/2 < r 
 
         test = bha(1)  
	h = xx(2) - xx(1)

  ! Set x, y, z  to coordinates of grid points
  do i=1,int(params(1))
     t(1,i,:,:)=xx(1)+(i-1)*h
  enddo
  do i=1,int(params(2))
     t(2,:,i,:)=yy(1)+(i-1)*h
  enddo
  do i=1,int(params(3))
     t(3,:,:,i)=zz(1)+(i-1)*h
  enddo

   t(4,:,:,:) = sqrt( t(1,:,:,:)**2 + t(2,:,:,:)**2 + t(3,:,:,:)**2 ) 

   where ( t(4,:,:,:) <= test )
      t(5,:,:,:) = 1.0d0 
   elsewhere
      t(5,:,:,:) = 1.0d0 + bha(1)/t(4,:,:,:) 
   end where

!  Compute metric from conformal factor and zero extrinsic curvature 
 
   g(1,:,:,:)  = t(5,:,:,:) **4 
   g(2,:,:,:)  = 0.0d0 
   g(3,:,:,:)  = 0.0d0 
   g(4,:,:,:)  = t(5,:,:,:) **4 
   g(5,:,:,:)  = 0.0d0 
   g(6,:,:,:)  = t(5,:,:,:) **4 

   k(1,:,:,:)  = 0.0d0 
   k(2,:,:,:)  = 0.0d0 
   k(3,:,:,:)  = 0.0d0 
   k(4,:,:,:)  = 0.0d0
   k(5,:,:,:)  = 0.0d0 
   k(6,:,:,:)  =  0.0d0

return 
end 
