!----------------------------------------------------------------------- 
! utility.f  
! Collection of utility routines for TCODE 
!----------------------------------------------------------------------- 
!  Routine to calculate the inverse 3-metric. Returns the determinant of 
!  the metric also.  
!  Routines require mask function for black hole case 
! 
!  routine verified  
! 
subroutine ginv(g,detg,mask,gup)
implicit none 
real*8,dimension(:,:,:,:)        :: g,gup
real*8,dimension(:,:,:)          :: detg,mask 
real*8,parameter                 :: zero = 0.0d0
real*8,parameter                 :: value  = 1.0d0

!interface  ----------------------------------------------------------------
interface
   subroutine calcdet(g11,g12,g13,g22,g23,g33,detg)
      implicit none
      real*8,dimension(:,:,:)    :: g11,g12,g13,g22,g23,g33
      real*8,dimension(:,:,:)    :: detg
	end subroutine calcdet

	subroutine replace_hole_val(array, chkval, setval, mask)
		implicit none
		real*8                        :: chkval, setval
		real*8,dimension(:,:,:)       :: array, mask
	end subroutine replace_hole_val 
		

end interface
!---------------------------------------------------------------------------


!   first calculate the determinant 
 
   call calcdet(g(1,:,:,:),g(2,:,:,:),g(3,:,:,:),g(4,:,:,:),g(5,:,:,:),   &
                g(6,:,:,:),detg) 
 
!    putting phony values for zero-determinant... 
         
   call replace_hole_val(detg,zero,value,mask) 
! 
 
!   construct the inverse component by component 
 
   gup(1,:,:,:) = (g(4,:,:,:)*g(6,:,:,:)-g(5,:,:,:)*g(5,:,:,:))/detg
   gup(2,:,:,:) = (g(3,:,:,:)*g(5,:,:,:)-g(2,:,:,:)*g(6,:,:,:))/detg
   gup(3,:,:,:) = (g(2,:,:,:)*g(5,:,:,:)-g(3,:,:,:)*g(4,:,:,:))/detg
   gup(4,:,:,:) = (g(1,:,:,:)*g(6,:,:,:)-g(3,:,:,:)*g(3,:,:,:))/detg
   gup(5,:,:,:) = (g(3,:,:,:)*g(2,:,:,:)-g(1,:,:,:)*g(5,:,:,:))/detg
   gup(6,:,:,:) = (g(1,:,:,:)*g(4,:,:,:)-g(2,:,:,:)*g(2,:,:,:))/detg
 
return 
end  
!----------------------------------------------------------------------- 
!  Routine to calculate the determinant of the 3-metric  
! 
!  routine verified  
 
subroutine calcdet(g11,g12,g13,g22,g23,g33,detg)
implicit none 
real*8,dimension(:,:,:)    :: g11,g12,g13,g22,g23,g33
real*8,dimension(:,:,:)    :: detg

   detg = g11*g22*g33 -  g11*g23*g23 -  g12*g12*g33  +  &
          g12*g13*g23 +  g13*g12*g23 -  g13*g13*g22 
 
return 
end 
