!----------------------------------------------------------------------- 
! utility.f  
! Collection of utility routines for TCODE 
!----------------------------------------------------------------------- 
subroutine err(routine,reason) 
implicit none 
character*(*)                 :: routine,reason 
 
   write(0,*) 'Error in',routine 
   write(0,*) 'routine failed since',reason 
 
return 
end 
!-------------------------------------------------------------------------- 
real*8 function n2i(n,a,mask,m,nl) 
implicit none 
integer                       :: n,nl
real*8,dimension(:,:,:)       :: a,mask
integer                       :: i,j,k,l 
integer,dimension(:)       :: m

write(6,*) 'inside subroutine'
write(6,*) a(3,3,3:5)


! first change mask
   do l = 1 , nl
		mask(m(l),m(nl+1),m(2*nl+1)) = 0.0d0
	end do


 
   n2i = 0.0d0 
 
   l = 0
   do k = 3 , n - 2 
   do j = 3 , n - 2 
   do i = 3 , n - 2 
      n2i = n2i + mask(i,j,k)*abs(a(i,j,k))
      l = l + int(mask(i,j,k))
   end do 
   end do 
   end do 
   
   n2i = n2i / l 

   do l = 1 , nl
		mask(m(l),m(nl+1),m(2*nl+1)) = 1.0d0
	end do
 
 
return 
end 
!------------------------------------------------------------------------- 
      real*8 function n2ib(nx,ny,nz,a,m,nl) 
      implicit none 
      integer nx,ny,nz 
      real*8 a(nx,ny,nz) 
      integer nl 
      integer m(:) 
      integer l,i,j,k 
 
      n2ib = 0.0d0 
      do l = 1 , nl 
 
         i  = m(l) 
         j  = m(nl+l) 
         k  = m(2*nl+l) 
 
         n2ib = n2ib + a(i,j,k)*a(i,j,k)
 
      end do 
 
      n2ib = sqrt( n2ib / nl ) 
 
      return 
      end 
!------------------------------------------------------------------------- 
! routine to raise a rank (0 2)  symmetric tensor to a rank (2 0) tensor 
! This routine returns the new symmetric tensor in another set of arrays. 
! 
subroutine raise02to20s(gup11, gup12, gup13, gup22, gup23, gup33,        &
                        tin11, tin12, tin13, tin22, tin23, tin33,        &
                        tout11, tout12, tout13, tout22, tout23, tout33) 
implicit none 
real*8                        :: h, dt
real*8,dimension(:,:,:)       :: gup11,gup12,gup13,gup22,gup23,gup33
real*8,dimension(:,:,:)       :: tin11,tin12,tin13,tin22,tin23,tin33 
real*8,dimension(:,:,:)       :: tout11,tout12,tout13,tout22,tout23,tout33 


   tout11 = gup11*gup11*tin11+ gup11*gup12*tin12+                        &
            gup11*gup13*tin13+ gup12*gup11*tin12+                        &
            gup12*gup12*tin22+ gup12*gup13*tin23+                        &
            gup13*gup11*tin13+ gup13*gup12*tin23+                        &
            gup13*gup13*tin33      
 
   tout12 = gup11*gup12*tin11+ gup11*gup22*tin12+                        &
            gup11*gup23*tin13+ gup12*gup12*tin12+                        &
            gup12*gup22*tin22+ gup12*gup23*tin23+                        &
            gup13*gup12*tin13+ gup13*gup22*tin23+                        &
            gup13*gup23*tin33 
 
   tout13 = gup11*gup13*tin11+ gup11*gup23*tin12+                        &
            gup11*gup33*tin13+ gup12*gup13*tin12+                        &
            gup12*gup23*tin22+ gup12*gup33*tin23+                        &
            gup13*gup13*tin13+ gup13*gup23*tin23+                        &
            gup13*gup33*tin33 
 
! T22 
   tout22 = gup12*gup12*tin11+ gup12*gup22*tin12+                        &
            gup12*gup23*tin13+ gup22*gup12*tin12+                        &
            gup22*gup22*tin22+ gup22*gup23*tin23+                        &
            gup23*gup12*tin13+ gup23*gup22*tin23+                        &
            gup23*gup23*tin33 
 
! T23 
  tout23 = gup12*gup13*tin11+  gup12*gup23*tin12+                        &
           gup12*gup33*tin13+  gup22*gup13*tin12+                        &
           gup22*gup23*tin22+  gup22*gup33*tin23+                        &
           gup23*gup13*tin13+  gup23*gup23*tin23+                        &
           gup23*gup33*tin33 
 
! T33 
  tout33 = gup13*gup13*tin11+  gup13*gup23*tin12+                        &
           gup13*gup33*tin13+  gup23*gup13*tin12+                        &
           gup23*gup23*tin22+  gup23*gup33*tin23+                        &
           gup33*gup13*tin13+  gup33*gup23*tin23+                        &
           gup33*gup33*tin33 
 
return 
end 
 
!----------------------------------------------------------------------- 
! routine to raise index of a vector (0 1) -> (1 0)  
! 
subroutine raise01i1( gup11, gup12, gup13, gup22, gup23, gup33,          &
                      t1, t2, t3, temp1, temp2, temp3) 
implicit none 
real*8,dimension(:,:,:)       :: gup11,gup12,gup13,gup22,gup23,gup33
real*8,dimension(:,:,:)      :: t1,t2,t3,temp1,temp2,temp3 

   temp1 = gup11 * t1 + gup12 * t2 +  gup13 * t3  
   temp2 = gup12 * t1 + gup22 * t2 +  gup23 * t3  
   temp3 = gup13 * t1 + gup23 * t2 +  gup33 * t3  
 
return 
end 
!----------------------------------------------------------------------- 
!  routine to raise first index of (0 2) tensor. notation->(contra cov)index 
! 
subroutine rais02i1(gup11, gup12, gup13, gup22, gup23, gup33,         &
                    t11, t12, t13, t21, t22, t23, t31, t32, t33,      &
                    temp1, temp2, temp3) 
implicit none 
real*8,dimension(:,:,:)    :: gup11,gup12,gup13,gup22,gup23,gup33
real*8,dimension(:,:,:)    :: t11,t12,t13,t21,t22,t23,t31,t32,t33
real*8,dimension(:,:,:)    :: temp1,temp2,temp3 

! {i}1 
   temp1 = t11
   temp2 = t21
   temp3 = t31
   t11 = gup11 * temp1  + gup12 * temp2 + gup13 * temp3 
   t21 = gup12 * temp1  + gup22 * temp2 + gup23 * temp3 
   t31 = gup13 * temp1  + gup23 * temp2 + gup33 * temp3 
 
! {i}2 
   temp1 = t12 
   temp2 = t22 
   temp3 = t22 
 
   t12 = gup11 * temp1  + gup12 * temp2 + gup13 * temp3 
   t22 = gup12 * temp1  + gup22 * temp2 + gup23 * temp3 
   t32 = gup13 * temp1  + gup23 * temp2 + gup33 * temp3 
 
! {i}3 
   temp1 = t13 
   temp2 = t23 
   temp3 = t33 
 
   t13 = gup11 * temp1  + gup12 * temp2 + gup13 * temp3 
   t23 = gup12 * temp1  + gup22 * temp2 + gup23 * temp3 
   t33 = gup13 * temp1  + gup23 * temp2 + gup33 * temp3 
               
 
return 
end 
!----------------------------------------------------------------------- 
!  routine to take the trace of a rank 0 2 symmetric tensor  
subroutine trace02sym(gup11, gup12, gup13, gup22, gup23, gup33,       &
                      t11, t12, t13, t22, t23, t33, trT) 
implicit none
real*8,dimension(:,:,:)    :: gup11,gup12,gup13,gup22,gup23,gup33
real*8,dimension(:,:,:)    :: t11,t12,t13,t22,t23,t33
real*8,dimension(:,:,:)    :: trT 

   trT = gup11 * t11 + gup12 * t12 * 2.0d0 + gup13 * t13 * 2.0d0      &
       + gup22 * t22 + gup23 * t23 * 2.0d0 + gup33 * t33  
 
return 
end 
!----------------------------------------------------------------------- 
! routine to lower index of a vector (1 0) -> (0 1)  
! 
subroutine lower10i1(g11, g12, g13, g22, g23, g33,                    &
                     t1, t2, t3, temp1, temp2, temp3) 
implicit none
real*8                        :: h, dt
real*8,dimension(:,:,:)       :: g11,g12,g13,g22,g23,g33
real*8,dimension(:,:,:)       :: t1, t2, t3
real*8,dimension(:,:,:)       :: temp1, temp2, temp3  

   temp1 = g11 * t1 + g12 * t2 +  g13 * t3  
   temp2 = g12 * t1 + g22 * t2 +  g23 * t3  
   temp3 = g13 * t1 + g23 * t2 +  g33 * t3  
 
return 
end 
!----------------------------------------------------------------------- 
!  routine to lower first index of (1 1) tensor. notation->(contra cov)index 
! 
subroutine lower11i1(g11, g12, g13, g22, g23, g33,                          &
                     t11, t12, t13, t21, t22, t23, t31, t32, t33,           &
                     temp1, temp2, temp3) 
implicit none
real*8                        :: h, dt
real*8,dimension(:,:,:)       :: g11,g12,g13,g22,g23,g33
real*8,dimension(:,:,:)       :: t11,t12,t13,t21,t22,t23,t31,t32,t33 
real*8,dimension(:,:,:)       :: temp1, temp2, temp3 

   temp1 = t11 
   temp2 = t21 
   temp3 = t31 
 
   t11 = g11*temp1 + g12*temp2 + g13*temp3 
   t21 = g12*temp1 + g22*temp2 + g23*temp3 
   t31 = g13*temp1 + g23*temp2 + g33*temp3 
 
   temp1 = t12 
   temp2 = t22 
   temp3 = t32 
 
   t12 = g11*temp1 + g12*temp2 + g13*temp3 
   t22 = g12*temp1 + g22*temp2 + g23*temp3 
   t32 = g13*temp1 + g23*temp2 + g33*temp3 
 
   temp1 = t13 
   temp2 = t23 
   temp3 = t33 

   t13 = g11*temp1 + g12*temp2 + g13*temp3 
   t23 = g12*temp1 + g22*temp2 + g23*temp3 
   t33 = g13*temp1 + g23*temp2 + g33*temp3 

return 
end 
!----------------------------------------------------------------------- 
subroutine replace_hole_val(array, chkval, setval, mask) 
implicit none 
real*8                        :: chkval, setval 
real*8,dimension(:,:,:)       :: array, mask 
 
   array = (mask-chkval)*array + (1.0d0-(mask-chkval))*setval 

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