************************************************************************
*     Randall R. Correll
*     Mijan F. Huq
*     Scott A. Klasky
*     Richard A. Matzner
*
*     COPYRIGHT 1995-1996

* Front end routine to initial data:
*   - calls appropriate subroutine
************************************************************************
		subroutine initdata(nx,ny,nz,g11,g12,g13,g22,g23,g33,
     &                    gup11,gup12,gup13,gup22,gup23,gup33,
     &                    k11,k12,k13,k21,k22,k23,k31,k32,k33,
     &                    alpha, beta1,beta2,beta3,
     &                    t11,temp1,temp2,temp3,xx,yy,zz,bha,bhr)

      implicit none

      integer   nx,ny,nz

c     metric:  g_ij
      real*8   g11(nx,ny,nz), g12(nx,ny,nz), g13(nx,ny,nz),
     &                        g22(nx,ny,nz), g23(nx,ny,nz),
     &                                       g33(nx,ny,nz)

      real*8   gup11(nx,ny,nz), gup12(nx,ny,nz), gup13(nx,ny,nz),
     &                          gup22(nx,ny,nz), gup23(nx,ny,nz),
     &                                           gup33(nx,ny,nz)

c     extrinsic curvature:  K^i_j
      real*8   k11(nx,ny,nz), k12(nx,ny,nz), k13(nx,ny,nz),
     &         k21(nx,ny,nz), k22(nx,ny,nz), k23(nx,ny,nz),
     &         k31(nx,ny,nz), k32(nx,ny,nz), k33(nx,ny,nz)
      real*8    xx(nx), yy(ny), zz(nz), bha(2), bhr(2,3)
      real*8   temp1(nx,ny,nz), temp2(nx,ny,nz), temp3(nx,ny,nz)
      real*8   t11(nx,ny,nz)

     	real*8   beta1(nx,ny,nz), beta2(nx,ny,nz), beta3(nx,ny,nz),
     &         alpha(nx,ny,nz)

c		call init_iso(nx,ny,nz,g11,g12,g13,g22,g23,g33,
c    &           k11,k12,k13,k21,k22,k23,k31,k32,k33,xx,yy,zz,bha,bhr,
c    &           temp1)

c      call init_ef(nx,ny,nz,g11,g12,g13,g22,g23,g33,
c     &           gup11,gup12,gup13,gup22,gup23,gup33,
c     &           k11,k12,k13,k21,k22,k23,k31,k32,k33,
c     &           alpha, beta1, beta2, beta3, xx, yy, zz, bha,bhr,
c     *           t11,temp1, temp2, temp3)

c		call init_hdf(nx,ny,nz,g11,g12,g13,g22,g23,g33,
c     &           gup11,gup12,gup13,gup22,gup23,gup33,
c     &           k11,k12,k13,k21,k22,k23,k31,k32,k33,
c     &           alpha, beta1, beta2, beta3, xx, yy, zz, bha,bhr,
c     *           t11,temp1, temp2, temp3)

		call init_ascii(nx,ny,nz,g11,g12,g13,g22,g23,g33,
     &           gup11,gup12,gup13,gup22,gup23,gup33,
     &           k11,k12,k13,k21,k22,k23,k31,k32,k33,
     &           alpha, beta1, beta2, beta3, xx, yy, zz, bha,bhr,
     *           t11,temp1, temp2, temp3)

		return
		end
************************************************************************
*  Routine to generate exact Isotropic schwarzschild initial data
************************************************************************
      subroutine init_iso(nx,ny,nz,g11,g12,g13,g22,g23,g33,
     &           k11,k12,k13,k21,k22,k23,k31,k32,k33,x,y,z,bha,bhr,
     &           phi)

      implicit none

      integer   nx,ny,nz
c     metric:  g_ij
      real*8   g11(nx,ny,nz), g12(nx,ny,nz), g13(nx,ny,nz),
     &                        g22(nx,ny,nz), g23(nx,ny,nz),
     &                                       g33(nx,ny,nz)
c     extrinsic curvature:  K^i_j
      real*8   k11(nx,ny,nz), k12(nx,ny,nz), k13(nx,ny,nz),
     &         k21(nx,ny,nz), k22(nx,ny,nz), k23(nx,ny,nz),
     &         k31(nx,ny,nz), k32(nx,ny,nz), k33(nx,ny,nz)
      real*8    x(nx), y(ny), z(nz), bha(2), bhr(2,3)
      real*8   phi(nx,ny,nz)

c     local variables
      integer  i,j,k
      real*8   test, r

c  if only one black hole, compute initial data as follows

c        compute conformal factor for a/2 < r

         test = 0.5e0 * bha(1)
         
         do k = 1, nz
         do j = 1, ny
         do i = 1, nx

            r = sqrt( (x(i)-bhr(1,1))**2
     &               +(y(j)-bhr(1,2))**2
     &               +(z(k)-bhr(1,3))**2 )
c           use masking function over singularity
            if ( r .le. test ) then
					phi(i,j,k) = 1.0e0
            else
               phi(i,j,k) = 1.0e0 + bha(1)/r
            end if

         end do
         end do
         end do


c  Compute metric from conformal factor and zero extrinsic curvature

      do k = 1, nz
      do j = 1, ny
      do i = 1, nx

         g11(i,j,k) = phi(i,j,k)**4
         g12(i,j,k) = 0.0e0
         g13(i,j,k) = 0.0e0
         g22(i,j,k) = phi(i,j,k)**4
         g23(i,j,k) = 0.0e0
         g33(i,j,k) = phi(i,j,k)**4

         k11(i,j,k) = 0.0e0
         k12(i,j,k) = 0.0e0
         k13(i,j,k) = 0.0e0
         k21(i,j,k) = 0.0e0
         k22(i,j,k) = 0.0e0
         k23(i,j,k) = 0.0e0
         k31(i,j,k) = 0.0e0
         k32(i,j,k) = 0.0e0
         k33(i,j,k) = 0.0e0

      end do
      end do
      end do

      return
      end
************************************************************************
* subroutine to generate initial data for Schwarzschild in 
* Eddington-Finkelstein coordinates
************************************************************************
      subroutine init_ef(nx,ny,nz,g11,g12,g13,g22,g23,g33,
     &           gup11,gup12,gup13,gup22,gup23,gup33,
     &           k11,k12,k13,k21,k22,k23,k31,k32,k33,
     &           alpha, beta1, beta2, beta3, xx, yy, zz, bha,bhr,
     *           mask,temp1, temp2, temp3)
		implicit none

		integer	nx, ny, nz

c     metric:  g_ij
      real*8   g11(nx,ny,nz), g12(nx,ny,nz), g13(nx,ny,nz),
     &                        g22(nx,ny,nz), g23(nx,ny,nz),
     &                                       g33(nx,ny,nz)

      real*8   gup11(nx,ny,nz), gup12(nx,ny,nz), gup13(nx,ny,nz),
     &                          gup22(nx,ny,nz), gup23(nx,ny,nz),
     &                                           gup33(nx,ny,nz)

c     extrinsic curvature:  K^i_j


      real*8   k11(nx,ny,nz), k12(nx,ny,nz), k13(nx,ny,nz),
     &         k21(nx,ny,nz), k22(nx,ny,nz), k23(nx,ny,nz),
     &         k31(nx,ny,nz), k32(nx,ny,nz), k33(nx,ny,nz)

c		lapse, shift and coordinates
     	real*8   beta1(nx,ny,nz), beta2(nx,ny,nz), beta3(nx,ny,nz),
     &         alpha(nx,ny,nz), xx(nx), yy(ny), zz(nz)

      real*8   temp1(nx,ny,nz), temp2(nx,ny,nz), temp3(nx,ny,nz)
      real*8   mask(nx,ny,nz)
      real*8   bha(2), bhr(2,3)

c	Local Variables

			integer	i,j,k,m,p
			real*8  rr, x,y,z, mass, fact

			mass = bha(1) * 0.5e0


			do i = 1, nx
				do j = 1, ny
					do k = 1, nz

                  rr = sqrt((xx(i)- bhr(1,1))**2 + (yy(j)-bhr(1,2))**2 
     &               + (zz(k)-bhr(1,3))**2)
					
					if(rr .eq. 0)then
						g11(i,j,k) = 1.0e0
						g12(i,j,k) = 0.0e0
						g13(i,j,k) = 0.0e0
						g22(i,j,k) = 1.0e0
						g23(i,j,k) = 0.0e0
						g33(i,j,k) = 1.0e0
						k11(i,j,k) = 0.0e0
						k12(i,j,k) = 0.0e0
						k13(i,j,k) = 0.0e0
						k21(i,j,k) = 0.0e0
						k22(i,j,k) = 0.0e0
						k23(i,j,k) = 0.0e0
						k31(i,j,k) = 0.0e0
						k32(i,j,k) = 0.0e0
						k33(i,j,k) = 0.0e0
					else

						x = xx(i) - bhr(1,1)
						y = yy(j) - bhr(1,2)
						z = zz(k) - bhr(1,3)

c				metric functions first

               g11(i,j,k) = 1.0e0 + 2.0e0*Mass * x**2 / rr**3

               g12(i,j,k) = 2.0e0 * Mass * x * y/ rr**3

               g13(i,j,k) = 2.0e0 * Mass * z * x/rr**3

               g22(i,j,k) = 1.0e0 + 2.0e0*Mass * y**2 / rr**3

               g23(i,j,k) = 2.0e0 * Mass * y * z/rr**3

               g33(i,j,k) = 1.0e0 + 2.0e0*Mass * z**2 / rr**3

               fact = 2.0e0 * Mass / rr**4 /sqrt(1.0e0 + 2.0e0*Mass/rr)

               k11(i,j,k) = (rr**2 - Mass*x**2 / rr -2.0e0 * x**2)*fact
		
               k12(i,j,k) = (-x*y*(2.0e0 + Mass / rr))*fact

               k13(i,j,k) = (-x*z*(2.0e0 + Mass / rr))*fact

               k21(i,j,k) = k12(i,j,k)
	
               k22(i,j,k) = (rr**2 - Mass*y**2 / rr -2.0e0 * y**2)*fact

               k23(i,j,k) = (-y*z*(2.0e0 + Mass / rr))*fact

               k31(i,j,k) = k13(i,j,k)

               k32(i,j,k) = k23(i,j,k)

               k33(i,j,k) = (rr**2 - Mass*z**2 / rr -2.0e0 * z**2)*fact

					end if

c           set temporary mask

            mask(i,j,k) = 1.0e0

c		Define lapse

            alpha(i,j,k) = 1.0e0 / sqrt(1.0e0 + 2.0e0*mass / rr)

c		Define shift vector: beta^i

            beta1(i,j,k) = +2.0e0*mass*x/(rr * (2.0e0*mass + rr))
            beta2(i,j,k) = +2.0e0*mass*y/(rr * (2.0e0*mass + rr))
            beta3(i,j,k) = +2.0e0*mass*z/(rr * (2.0e0*mass + rr))

					end do
				end do
			end do

c     Raise K_ij -> K^i_j

               call inverse(g11,g12,g13,g22,g23,g33,
     *              gup11,gup12,gup13,gup22, gup23, gup33,temp1,
     *              mask,nx, ny, nz)

               call rais02i1(nx, ny, nz,
     *              gup11, gup12, gup13, gup22, gup23, gup33,
     *              k11, k12, k13, k21, k22, k23, k31, k32, k33,
     *              temp1, temp2, temp3)

			return
		end

c-----------------------------------------------------------------------
		subroutine init_hdf(nx,ny,nz,
     &           g11,g12,g13,g22,g23,g33,
     &           gup11,gup12,gup13,gup22,gup23,gup33,
     &           k11,k12,k13,k21,k22,k23,k31,k32,k33,
     &           alpha, beta1, beta2, beta3, xx, yy, zz, bha,bhr,
     *           mask,temp1, temp2, temp3)
		   implicit none
			integer	nx, ny, nz

c     metric:  g_ij
      real*8   g11(nx,ny,nz), g12(nx,ny,nz), g13(nx,ny,nz),
     &                        g22(nx,ny,nz), g23(nx,ny,nz),
     &                                       g33(nx,ny,nz)

      real*8   gup11(nx,ny,nz), gup12(nx,ny,nz), gup13(nx,ny,nz),
     &                          gup22(nx,ny,nz), gup23(nx,ny,nz),
     &                                           gup33(nx,ny,nz)

c     extrinsic curvature:  K^i_j


      real*8   k11(nx,ny,nz), k12(nx,ny,nz), k13(nx,ny,nz),
     &         k21(nx,ny,nz), k22(nx,ny,nz), k23(nx,ny,nz),
     &         k31(nx,ny,nz), k32(nx,ny,nz), k33(nx,ny,nz)

c		lapse, shift and coordinates
     	real*8   beta1(nx,ny,nz), beta2(nx,ny,nz), beta3(nx,ny,nz),
     &         alpha(nx,ny,nz), xx(nx), yy(ny), zz(nz)

      real*8   temp1(nx,ny,nz), temp2(nx,ny,nz), temp3(nx,ny,nz)
      real*8   mask(nx,ny,nz)
      real*8   bha(2), bhr(2,3)

c		local variables

		integer	i, j, k, level, indlnb, rank

		real*8	time

		character*32 cname1, cname2, cname3, name

		write(0,*)'In readhdf3d: ',nx,ny,nz

		write(name,'(a11)')'psiinit.hdf'
		write(cname1,'(a1)')'x'
		write(cname2,'(a1)')'y'
		write(cname3,'(a1)')'z'
		level = 1
		time = 0.0d0

c			Set up to read in psi 
			call readhdf3d(nx, ny, nz, temp1, xx, yy, zz, time,
     .   cname1, cname2, cname3, level,name(1:indlnb(name)),temp2)

c			Set up g11
			do i = 1, nx
			  do j = 1, ny
			    do k = 1, nz
					g11(i,j,k) = temp1(i,j,k)**4
			    end do
			  end do
			end do

c			Set up g22
			do i = 1, nx
			  do j = 1, ny
			    do k = 1, nz
					g22(i,j,k) = temp1(i,j,k)**4
			    end do
			  end do
			end do

c			Set up g33
			do i = 1, nx
			  do j = 1, ny
			    do k = 1, nz
					g33(i,j,k) = temp1(i,j,k)**4
			    end do
			  end do
			end do

c			Zero out other components of the metric

			do i = 1, nx
			  do j = 1, ny
			    do k = 1, nz
					g12(i,j,k) = 0.0
			    end do
			  end do
			end do

			do i = 1, nx
			  do j = 1, ny
			    do k = 1, nz
					g13(i,j,k) = 0.0
			    end do
			  end do
			end do

			do i = 1, nx
			  do j = 1, ny
			    do k = 1, nz
					g23(i,j,k) = 0.0
			    end do
			  end do
			end do

c		Read in covariant form of extrinsic curvature tensor
		write(name,'(a11)')'k11init.hdf'
		write(cname1,'(a1)')'x'
		write(cname2,'(a1)')'y'
		write(cname3,'(a1)')'z'
			call readhdf3d(nx, ny, nz, k11, xx, yy, zz, time,
     .    cname1, cname2, cname3, level,name(1:indlnb(name)),temp2)

		write(name,'(a11)')'k12init.hdf'
		write(cname1,'(a1)')'x'
		write(cname2,'(a1)')'y'
		write(cname3,'(a1)')'z'
			call readhdf3d(nx, ny, nz, k12, xx, yy, zz, time,
     .    cname1, cname2, cname3, level,name(1:indlnb(name)),temp2)

		write(name,'(a11)')'k13init.hdf'
		write(cname1,'(a1)')'x'
		write(cname2,'(a1)')'y'
		write(cname3,'(a1)')'z'
			call readhdf3d(nx, ny, nz, k13, xx, yy, zz, time,
     .    cname1, cname2, cname3, level,name(1:indlnb(name)),temp2)

		write(name,'(a11)')'k22init.hdf'
		write(cname1,'(a1)')'x'
		write(cname2,'(a1)')'y'
		write(cname3,'(a1)')'z'
			call readhdf3d(nx, ny, nz, k22, xx, yy, zz, time,
     .    cname1, cname2, cname3, level,name(1:indlnb(name)),temp2)

		write(name,'(a11)')'k23init.hdf'
		write(cname1,'(a1)')'x'
		write(cname2,'(a1)')'y'
		write(cname3,'(a1)')'z'
			call readhdf3d(nx, ny, nz, k23, xx, yy, zz, time,
     .   cname1, cname2, cname3, level,name(1:indlnb(name)),temp2)

		write(name,'(a11)')'k33init.hdf'
		write(cname1,'(a1)')'x'
		write(cname2,'(a1)')'y'
		write(cname3,'(a1)')'z'
			call readhdf3d(nx, ny, nz, k33, xx, yy, zz, time,
     .   cname1, cname2, cname3, level,name(1:indlnb(name)),temp2)

c		Set mask to 1 everywhere or could read in the mask this way too.

			do i = 1, nx
			  do j = 1, ny
			    do k = 1, nz
			      mask(i,j,k) = 1.0d0
			    end do
			  end do
			end do

c     Raise K_ij -> K^i_j

               call inverse(g11,g12,g13,g22,g23,g33,
     *              gup11,gup12,gup13,gup22, gup23, gup33,temp1,
     *              mask,nx, ny, nz)

               call rais02i1(nx, ny, nz,
     *              gup11, gup12, gup13, gup22, gup23, gup33,
     *              k11, k12, k13, k21, k22, k23, k31, k32, k33,
     *              temp1, temp2, temp3)



		  return
		end
c-----------------------------------------------------------------------
		subroutine init_ascii(nx,ny,nz,
     &           g11,g12,g13,g22,g23,g33,
     &           gup11,gup12,gup13,gup22,gup23,gup33,
     &           k11,k12,k13,k21,k22,k23,k31,k32,k33,
     &           alpha, beta1, beta2, beta3, xx, yy, zz, bha,bhr,
     *           mask,temp1, temp2, temp3)
		   implicit none
			integer	nx, ny, nz

c     metric:  g_ij
      real*8   g11(nx,ny,nz), g12(nx,ny,nz), g13(nx,ny,nz),
     &                        g22(nx,ny,nz), g23(nx,ny,nz),
     &                                       g33(nx,ny,nz)

      real*8   gup11(nx,ny,nz), gup12(nx,ny,nz), gup13(nx,ny,nz),
     &                          gup22(nx,ny,nz), gup23(nx,ny,nz),
     &                                           gup33(nx,ny,nz)

c     extrinsic curvature:  K^i_j


      real*8   k11(nx,ny,nz), k12(nx,ny,nz), k13(nx,ny,nz),
     &         k21(nx,ny,nz), k22(nx,ny,nz), k23(nx,ny,nz),
     &         k31(nx,ny,nz), k32(nx,ny,nz), k33(nx,ny,nz)

c		lapse, shift and coordinates
     	real*8   beta1(nx,ny,nz), beta2(nx,ny,nz), beta3(nx,ny,nz),
     &         alpha(nx,ny,nz), xx(nx), yy(ny), zz(nz)

      real*8   temp1(nx,ny,nz), temp2(nx,ny,nz), temp3(nx,ny,nz)
      real*8   mask(nx,ny,nz)
      real*8   bha(2), bhr(2,3)

c		local variables

		integer	i, j, k, level, indlnb, rank

		real*8	time

		character*32 cname1, cname2, cname3, name

		write(0,*)'In readhdf3d: ',nx,ny,nz

		write(name,'(a11)')'psiinit.dat'

c			Set up to read in psi 

			call read_ascii(nx,ny,nz,name,temp1)

c			Set up g11
			do i = 1, nx
			  do j = 1, ny
			    do k = 1, nz
					g11(i,j,k) = temp1(i,j,k)**4
			    end do
			  end do
			end do

c			Set up g22
			do i = 1, nx
			  do j = 1, ny
			    do k = 1, nz
					g22(i,j,k) = temp1(i,j,k)**4
			    end do
			  end do
			end do

c			Set up g33
			do i = 1, nx
			  do j = 1, ny
			    do k = 1, nz
					g33(i,j,k) = temp1(i,j,k)**4
			    end do
			  end do
			end do

c			Zero out other components of the metric

			do i = 1, nx
			  do j = 1, ny
			    do k = 1, nz
					g12(i,j,k) = 0.0
			    end do
			  end do
			end do

			do i = 1, nx
			  do j = 1, ny
			    do k = 1, nz
					g13(i,j,k) = 0.0
			    end do
			  end do
			end do

			do i = 1, nx
			  do j = 1, ny
			    do k = 1, nz
					g23(i,j,k) = 0.0
			    end do
			  end do
			end do

c		Read in covariant form of extrinsic curvature tensor

		call read_ascii(nx,ny,nz,name,temp1)

		write(name,'(a11)')'k11init.dat'
		call read_ascii(nx,ny,nz,name,k11)

		write(name,'(a11)')'k12init.dat'
		call read_ascii(nx,ny,nz,name,k12)

		write(name,'(a11)')'k13init.dat'
		call read_ascii(nx,ny,nz,name,k13)

		write(name,'(a11)')'k22init.dat'
		call read_ascii(nx,ny,nz,name,k22)

		write(name,'(a11)')'k23init.dat'
		call read_ascii(nx,ny,nz,name,k23)

		write(name,'(a11)')'k33init.dat'
		call read_ascii(nx,ny,nz,name,k33)

c		Set mask to 1 everywhere or could read in the mask this way too.

			do i = 1, nx
			  do j = 1, ny
			    do k = 1, nz
			      mask(i,j,k) = 1.0d0
			    end do
			  end do
			end do

c     Raise K_ij -> K^i_j

               call inverse(g11,g12,g13,g22,g23,g33,
     *              gup11,gup12,gup13,gup22, gup23, gup33,temp1,
     *              mask,nx, ny, nz)

               call rais02i1(nx, ny, nz,
     *              gup11, gup12, gup13, gup22, gup23, gup33,
     *              k11, k12, k13, k21, k22, k23, k31, k32, k33,
     *              temp1, temp2, temp3)



		  return
		end
c-----------------------------------------------------------------------
