!/*@@
!
! $Id: ricci_l2.f90,v 1.4 1997/03/07 14:19:32 scott Exp $
!
!@@*/

!/*@@
!   @routine ricci
!   @author M.Huq S.Klasky 
!   @desc
!   This routine returns covariant form of ricci tensor. Raise index before
!   using in evolution code.
!   Overall routine checked to be numerically consistent
!     - 2nd convergent for ricci1-ricci6
!     - ricci7 correct to roune-off
!     - routines checked with exact connections
!   @enddesc
!   @calledby tstep
!   @calls diff
!   @comment
!   uses module tensors and diffinterface
!   @endcomment
!@@*/
!-----------------------------------------------------------------------

subroutine ricci(params,g,gup,c,ca,h,ri,m,d,nl,r)
use input_params
use tensors
use diffinterface
implicit none
real*8                                        :: h
integer                              :: nl
real*8,dimension(:)                  :: d,params
integer,dimension(:)                 :: m
type(GridFunction),dimension(3,3)             :: gup,g 
type(GridFunction),dimension(3,3)             :: r 
type(GridFunction),dimension(3,3,3)           :: c,ca
type(GridFunction),dimension(3,3,3,3)         :: ri 

integer  :: i, j, k, l,n,p,q
integer,dimension(3,3)                        :: wd
integer, dimension(3,3,3,3),parameter                   ::symmetry= &
  reshape( (/0,0,0,0,0,0,0,0,0, 0,1,1,-1,0,1,-1,-1,0, &
       & 0,1,1,-1,0,1,-1,-1,0, &
       & 0,-1,-1,1,0,-1,1,1,0,0,0,0,0,0,0,0,0,0,0,1,1,-1,0,1,-1,-1,0,&
       & 0,-1,-1,1,0,-1,1,1,0,0,-1,-1,1,0,-1,1,1,0,0,0,0,0,0,0,0,0,0 /),&
       & (/3,3,3,3/))



  wd(1,1) = 4
  wd(1,2) = 7
  wd(1,3) = 8
  wd(2,1) = 7
  wd(2,2) = 5
  wd(2,3) = 9
  wd(3,1) = 8
  wd(3,2) = 9
  wd(3,3) = 6

   call diff(params,g(1,2)%d,g(1,2)%d,1.0d0,wd(1,2),h, &
             int(params(P_BNDRY_DERIV_ORDER)), &
             2,m,d,nl,r(1,1)%d)
   call diff(params,g(1,1)%d,g(1,1)%d,-0.5d0,wd(2,2),h, &
             int(params(P_BNDRY_DERIV_ORDER)), &
             2,m,d,nl,r(1,2)%d)
   call diff(params,g(2,2)%d,g(2,2)%d,-0.5d0,wd(1,1),h, &
             int(params(P_BNDRY_DERIV_ORDER)), &
             2,m,d,nl,r(1,3)%d)
	ri(1,2,1,2)%d = r(1,1)%d + r(1,2)%d + r(1,3)%d 

   call diff(params,g(1,3)%d,g(1,3)%d, 0.5d0,wd(1,2),h, &
             int(params(P_BNDRY_DERIV_ORDER)), &
             2,m,d,nl,r(1,1)%d)
   call diff(params,g(2,1)%d,g(2,1)%d, 0.5d0,wd(1,3),h, &
             int(params(P_BNDRY_DERIV_ORDER)), &
              2,m,d,nl,r(1,2)%d)
   call diff(params,g(1,1)%d,g(1,1)%d,-0.5d0,wd(2,3),h, &
             int(params(P_BNDRY_DERIV_ORDER)), &
              2,m,d,nl,r(1,3)%d)
   call diff(params,g(2,3)%d,g(2,3)%d,-0.5d0,wd(1,1),h, &
             int(params(P_BNDRY_DERIV_ORDER)), &
             2,m,d,nl,r(2,2)%d)
	ri(1,2,1,3)%d = r(1,1)%d + r(1,2)%d + r(1,3)%d  + r(2,2)%d

   call diff(params,g(1,3)%d,g(1,3)%d, 0.5d0,wd(2,2),h, &
             int(params(P_BNDRY_DERIV_ORDER)), &
             2,m,d,nl,r(1,1)%d)
   call diff(params,g(2,2)%d,g(2,2)%d, 0.5d0,wd(1,3),h, &
             int(params(P_BNDRY_DERIV_ORDER)),2 &
             ,m,d,nl,r(1,2)%d)
   call diff(params,g(1,2)%d,g(1,2)%d,-0.5d0,wd(2,3),h, &
             int(params(P_BNDRY_DERIV_ORDER)),2, &
             m,d,nl,r(1,3)%d)
   call diff(params,g(2,3)%d,g(2,3)%d,-0.5d0,wd(1,2),h, &
             int(params(P_BNDRY_DERIV_ORDER)),2, &
              m,d,nl,r(2,2)%d)
	ri(1,2,2,3)%d = r(1,1)%d + r(1,2)%d + r(1,3)%d  + r(2,2)%d

   call diff(params,g(1,3)%d,g(1,3)%d, 1.0d0,wd(1,3),h, &
             int(params(P_BNDRY_DERIV_ORDER)), &
             2,m,d,nl,r(1,1)%d)
   call diff(params,g(1,1)%d,g(1,1)%d,-0.5d0,wd(3,3),h, &
             int(params(P_BNDRY_DERIV_ORDER)), &
             2,m,d,nl,r(1,2)%d)
   call diff(params,g(3,3)%d,g(3,3)%d,-0.5d0,wd(1,1),h, &
             int(params(P_BNDRY_DERIV_ORDER)), &
             2,m,d,nl,r(1,3)%d)
	ri(1,3,1,3)%d = r(1,1)%d + r(1,2)%d + r(1,3)%d 

   call diff(params,g(1,3)%d,g(1,3)%d, 0.5d0,wd(2,3),h, &
             int(params(P_BNDRY_DERIV_ORDER)), &
             2,m,d,nl,r(1,1)%d)
   call diff(params,g(2,3)%d,g(2,3)%d, 0.5d0,wd(1,3),h, &
             int(params(P_BNDRY_DERIV_ORDER)), &
             2,m,d,nl,r(1,2)%d)
   call diff(params,g(1,2)%d,g(1,2)%d,-0.5d0,wd(3,3),h, &
             int(params(P_BNDRY_DERIV_ORDER)), &
             2,m,d,nl,r(1,3)%d)
   call diff(params,g(3,3)%d,g(3,3)%d,-0.5d0,wd(1,2),h, &
             int(params(P_BNDRY_DERIV_ORDER)), &
             2,m,d,nl,r(2,2)%d)
	ri(1,3,2,3)%d = r(1,1)%d + r(1,2)%d + r(1,3)%d  + r(2,2)%d

   call diff(params,g(2,3)%d,g(2,3)%d, 1.0d0,wd(2,3),h, &
             int(params(P_BNDRY_DERIV_ORDER)), &
              2,m,d,nl,r(1,1)%d)
   call diff(params,g(2,2)%d,g(2,2)%d,-0.5d0,wd(3,3),h, &
             int(params(P_BNDRY_DERIV_ORDER)), &
              2,m,d,nl,r(1,2)%d)
   call diff(params,g(3,3)%d,g(3,3)%d,-0.5d0,wd(2,2),h, &
             int(params(P_BNDRY_DERIV_ORDER)), &
             2,m,d,nl,r(1,3)%d)
	ri(2,3,2,3)%d = r(1,1)%d + r(1,2)%d + r(1,3)%d 

!
! Add the cfel terms to ri.  The funny limits on the do loops
! are to insure that we only calculate the non-zero, independent parts of
! ri.
!



  do i=1,2
     do k=i+1,3
        do n=k,3
           do l=i,n-1
              do p=1,3
                 ri(i,k,l,n)%d=ri(i,k,l,n)%d &
                      & + ca(p,k,l)%d*c(p,i,n)%d &
                      & - ca(p,k,n)%d*c(p,i,l)%d
              enddo
           enddo
        enddo
     enddo
  enddo

  do k=1,3
     do n=1,k
        r(k,n)%d=0.d0
        do i=1,3
           do l=1,3
              if(i/=k .and. l/=n) then
                 r(k,n)%d=r(k,n)%d+gup(i,l)%d*symmetry(i,k,l,n) &
                      & * ri(i,k,l,n)%d
              endif
           enddo
        enddo
     enddo
  enddo

return 
end
