!*********************************************************************** 
!  Computes the lapse and shift vector 
!*********************************************************************** 
subroutine lapse(x,y,z,g,k,gup,t,mask,params,alpha)
implicit none
real*8,dimension(:,:,:,:)            :: g,k,gup,t
real*8,dimension(:,:,:)              :: mask
real*8,dimension(:,:,:)              :: alpha 
real*8,dimension(:)                  :: x,y,z
real*8,dimension(:)                  :: params
real*8                               :: bha(2),bhr(3,2)

interface

   subroutine lapse_iso(nx,ny,nz,mask,xx,yy,zz,x,y,z,r,bha,bhr,alpha)
      implicit none
      integer                       :: nx, ny, nz
      real*8,dimension(:,:,:)       :: mask, alpha
      real*8,dimension(:)           :: xx,yy,zz
      real*8,dimension(:,:,:)       :: x,y,z,r
      real*8                        :: bha(2),bhr(2,3)
      integer                       :: i,j,k
      real*8                        :: r1,r2,alpha1,alpha2
      real*8                        :: h
	end subroutine
end interface


   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 (params(21).eq.1) then !isometric Schwarschild lapse and gauge

		call lapse_iso(int(params(1)),int(params(2)),int(params(3)),      &
                        mask,x,y,z,t(1,:,:,:),t(2,:,:,:),               &
                        t(3,:,:,:),t(4,:,:,:),bha,bhr,alpha)
	else
		write(6,*) '*** error gauge condition unknown ***'
		stop

	end if

return
end
!*********************************************************************** 
subroutine shift(x,y,z,g,k,gup,t,mask,params,beta)
implicit none
real*8,dimension(:,:,:,:),target     :: g,k,gup,t
real*8,dimension(:,:,:)              :: mask
real*8,dimension(:,:,:,:)            :: beta
real*8,dimension(:)                  :: x,y,z
real*8,dimension(:)                  :: params
real*8                               :: bha(2),bhr(3,2)

   if (params(21).eq.1) then !isometric Schwarschild lapse and gauge

	beta  = 0.0d0

   else
      write(6,*) '*** error gauge condition unknown ***'
      stop

   end if

return
end

!*********************************************************************** 
subroutine lapse_iso(nx,ny,nz,mask,xx,yy,zz,x,y,z,r,bha,bhr,alpha) 
 
implicit none 
integer                       :: nx, ny, nz
real*8,dimension(:,:,:)       :: mask, alpha
real*8,dimension(:)           :: xx,yy,zz
real*8,dimension(:,:,:)       :: x,y,z,r
real*8                        :: bha(2),bhr(2,3)
integer                       :: i,j,k
real*8                        :: r1,r2,alpha1,alpha2
real*8                        :: h


	h = xx(2) - xx(1)

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

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

	alpha = mask * (r - bha(1)) / (r + bha(1))

 
return 
end 
!*********************************************************************** 
