!-----------------------------------------------------------------------
! ricci.f
!
! overall routine checked to be numerically consistent
!  - 2nd convergent for ricci1-ricci6
!  - ricci7 correct to round-off
!  - routines checked with exact connections
!
!-----------------------------------------------------------------------
! This routine returns covariant form of ricci tensor. Raise index before
! using in evolution code.

subroutine ricci(g,gup,c,ca,h,t,rieman,m,d,nl,r)
implicit none

integer                              :: nl
real*8                               :: h
real*8,dimension(:,:,:,:)            :: g,gup,c,t,r
real*8,dimension(:,:,:,:)            :: rieman,ca 
real*8,dimension(:)               :: d
integer,dimension(:)              :: m

! interface---------------------------------------------------------
interface
   subroutine ricci1(c111, c112, c113, c122, c123, c133,        &
                     c211, c212, c213, c222, c223, c233,        &
                     c311, c312, c313, c322, c323, c333,        &
                     c111d1, c112d1, c113d1, c121d1, c122d1,    &
                     c123d1, c131d1, c132d1, c133d1,            &
                     r11, r12, r13, r22, r23, r33,              &
                     m,d,nl, h)
   implicit none
   integer                       :: nl
   real*8                        :: h
   real*8,dimension(:,:,:)    :: c111,c112,c113,c122,c123,c133
   real*8,dimension(:,:,:)    :: c211,c212,c213,c222,c223,c233
   real*8,dimension(:,:,:)    :: c311,c312,c313,c322,c323,c333
   real*8,dimension(:,:,:)    :: c111d1,c112d1,c113d1
   real*8,dimension(:,:,:)    :: c121d1,c122d1,c123d1
   real*8,dimension(:,:,:)    :: c131d1,c132d1,c133d1
   real*8,dimension(:,:,:)    :: r11,r12,r13,r22,r23,r33
   real*8,dimension(:)        :: d
   integer,dimension(:)       :: m
   integer                       :: i,j,k
   end subroutine
   subroutine ricci2(c111, c112, c113, c122, c123, c133,        &
                     c211, c212, c213, c222, c223, c233,        &
                     c311, c312, c313, c322, c323, c333,        &
                     c211d2, c212d2, c213d2, c221d2, c222d2,    &
                     c223d2, c231d2, c232d2, c233d2,            &
                     r11, r12, r13, r22, r23, r33,              &
                     m,d,nl,h)
      implicit none
      integer                       :: nl
      real*8                        :: h
      real*8,dimension(:,:,:)    :: c111,c112,c113,c122,c123,c133
      real*8,dimension(:,:,:)    :: c211,c212,c213,c222,c223,c233
      real*8,dimension(:,:,:)    :: c311,c312,c313,c322,c323,c333
      real*8,dimension(:,:,:)    :: c211d2,c212d2,c213d2
      real*8,dimension(:,:,:)    :: c221d2,c222d2,c223d2
      real*8,dimension(:,:,:)    :: c231d2,c232d2,c233d2
      real*8,dimension(:,:,:)    :: r11,r12,r13,r22,r23,r33
      real*8,dimension(:)        :: d
      integer,dimension(:)       :: m
      integer                       :: i,j,k
   end subroutine
   subroutine ricci3(c111, c112, c113, c122, c123, c133,         &
                     c211, c212, c213, c222, c223, c233,         &
                     c311, c312, c313, c322, c323, c333,         &
                     c311d3, c312d3, c313d3, c321d3, c322d3,     &
                     c323d3, c331d3, c332d3, c333d3,             &
                     r11, r12, r13, r22, r23, r33,               &
                     m,d,nl, h)
      implicit none
      integer                       :: nl
      real*8                        :: h
      real*8,dimension(:,:,:)    :: c111,c112,c113,c122,c123,c133
      real*8,dimension(:,:,:)    :: c211,c212,c213,c222,c223,c233
      real*8,dimension(:,:,:)    :: c311,c312,c313,c322,c323,c333
      real*8,dimension(:,:,:)    :: c311d3,c312d3,c313d3
      real*8,dimension(:,:,:)    :: c321d3,c322d3,c323d3
      real*8,dimension(:,:,:)    :: c331d3,c332d3,c333d3
      real*8,dimension(:,:,:)    :: r11,r12,r13,r22,r23,r33
      real*8,dimension(:)        :: d
      integer,dimension(:)       :: m
      integer                       :: i,j,k
   end subroutine
   subroutine ricci4(c111, c112, c113, c122, c123, c133,   &
                     c211, c212, c213, c222, c223, c233,   &
                     c311, c312, c313, c322, c323, c333,   &
                     c111d1, c111d2, c111d3, c121d1, c121d2,    &
                     c121d3, c131d1, c131d2, c131d3,   &
                     r11, r12, r13, r22, r23, r33,   &
                     m,d,nl, h)
      implicit none
      integer                       :: nl
      real*8                        :: h
      real*8,dimension(:,:,:)    :: c111,c112,c113,c122,c123,c133
      real*8,dimension(:,:,:)    :: c211,c212,c213,c222,c223,c233
      real*8,dimension(:,:,:)    :: c311,c312,c313,c322,c323,c333
      real*8,dimension(:,:,:)    :: c111d1, c111d2, c111d3
      real*8,dimension(:,:,:)    :: c121d1, c121d2, c121d3
      real*8,dimension(:,:,:)    :: c131d1, c131d2, c131d3
      real*8,dimension(:,:,:)    :: r11,r12,r13,r22,r23,r33
      real*8,dimension(:)        :: d
      integer,dimension(:)       :: m
      integer                       :: i,j,k
   end subroutine
   subroutine ricci5(c111, c112, c113, c122, c123, c133,         &
                     c211, c212, c213, c222, c223, c233,         &
                     c311, c312, c313, c322, c323, c333,         &
                     c212d1, c212d2, c212d3, c222d1, c222d2,     &
                     c222d3, c232d1, c232d2, c232d3,             &
                     r11, r12, r13, r22, r23, r33,               &
                     m,d,nl, h)
      implicit none
      integer                       :: nl
      real*8                        :: h
      real*8,dimension(:,:,:)    :: c111,c112,c113,c122,c123,c133
      real*8,dimension(:,:,:)    :: c211,c212,c213,c222,c223,c233
      real*8,dimension(:,:,:)    :: c311,c312,c313,c322,c323,c333
      real*8,dimension(:,:,:)    :: c212d1, c212d2, c212d3
      real*8,dimension(:,:,:)    :: c222d1, c222d2, c222d3
      real*8,dimension(:,:,:)    :: c232d1, c232d2, c232d3
      real*8,dimension(:,:,:)    :: r11,r12,r13,r22,r23,r33
      real*8,dimension(:)        :: d
      integer,dimension(:)       :: m
      integer                       :: i,j,k
   end subroutine
   subroutine ricci6(c111, c112, c113, c122, c123, c133,         &
                     c211, c212, c213, c222, c223, c233,         &
                     c311, c312, c313, c322, c323, c333,         &
                     c313d1, c313d2, c313d3, c323d1, c323d2,     &
                     c323d3, c333d1, c333d2, c333d3,             &
                     r11, r12, r13, r22, r23, r33,               &
                     m,d,nl, h)
      implicit none
      integer                       :: nl
      real*8                        :: h
      real*8,dimension(:,:,:)    :: c111,c112,c113,c122,c123,c133
      real*8,dimension(:,:,:)    :: c211,c212,c213,c222,c223,c233
      real*8,dimension(:,:,:)    :: c311,c312,c313,c322,c323,c333
      real*8,dimension(:,:,:)    :: c313d1, c313d2, c313d3
      real*8,dimension(:,:,:)    :: c323d1, c323d2, c323d3
      real*8,dimension(:,:,:)    :: c333d1, c333d2, c333d3
      real*8,dimension(:,:,:)    :: r11,r12,r13,r22,r23,r33
      real*8,dimension(:)        :: d
      integer,dimension(:)       :: m
      integer                       :: i,j,k
   end subroutine
   subroutine ricci7(c111, c112, c113, c122, c123, c133,    &
                     c211, c212, c213, c222, c223, c233,    &
                     c311, c312, c313, c322, c323, c333,    &
                     t1, t2, t3, t4, t5, t6, t7, t8, t9,    &
                     r11, r12, r13, r22, r23, r33,          &
                     m,d,nl, h)
      implicit none
      integer                       :: nl
      real*8                        :: h
      real*8,dimension(:,:,:)    :: c111,c112,c113,c122,c123,c133
      real*8,dimension(:,:,:)    :: c211,c212,c213,c222,c223,c233
      real*8,dimension(:,:,:)    :: c311,c312,c313,c322,c323,c333
      real*8,dimension(:,:,:)    :: t1,t2,t3,t4,t5,t6,t7,t8,t9
      real*8,dimension(:,:,:)    :: r11,r12,r13,r22,r23,r33
      real*8,dimension(:)        :: d
      integer,dimension(:)       :: m
      integer                       :: i,j,k
   end subroutine
end interface
      ! end interface---------------------------------------------------------

!      call subsidiary routines to calculate ricci tensor piecewise.

   call ricci1(c(1,:,:,:),c(2,:,:,:),c(3,:,:,:),c(4,:,:,:),c(5,:,:,:),      &
               c(6,:,:,:),c(7,:,:,:),c(8,:,:,:),c(9,:,:,:),c(10,:,:,:),     &
               c(11,:,:,:),c(12,:,:,:),                                   &
               c(13,:,:,:),c(14,:,:,:),c(15,:,:,:),                      &
               c(16,:,:,:),c(17,:,:,:),c(18,:,:,:),t(1,:,:,:),t(2,:,:,:),&
               t(3,:,:,:),t(4,:,:,:),t(5,:,:,:),t(6,:,:,:),t(7,:,:,:),      &
               t(8,:,:,:),t(9,:,:,:),r(1,:,:,:),r(2,:,:,:),r(3,:,:,:),      &
               r(4,:,:,:), r(5,:,:,:), r(6,:,:,:),m,d,nl, h)

   call ricci2(c(1,:,:,:),c(2,:,:,:),c(3,:,:,:),c(4,:,:,:),c(5,:,:,:),      &
               c(6,:,:,:),c(7,:,:,:),c(8,:,:,:),c(9,:,:,:),c(10,:,:,:),     &
               c(11,:,:,:),c(12,:,:,:),                                   &
               c(13,:,:,:),c(14,:,:,:),c(15,:,:,:),                      &
               c(16,:,:,:),c(17,:,:,:),c(18,:,:,:),t(1,:,:,:),t(2,:,:,:),&
               t(3,:,:,:),t(4,:,:,:),t(5,:,:,:),t(6,:,:,:),t(7,:,:,:),      &
               t(8,:,:,:),t(9,:,:,:),r(1,:,:,:),r(2,:,:,:),r(3,:,:,:),      &
               r(4,:,:,:), r(5,:,:,:), r(6,:,:,:),m,d,nl, h)

   call ricci3(c(1,:,:,:),c(2,:,:,:),c(3,:,:,:),c(4,:,:,:),c(5,:,:,:),      &
               c(6,:,:,:),c(7,:,:,:),c(8,:,:,:),c(9,:,:,:),c(10,:,:,:),     &
               c(11,:,:,:),c(12,:,:,:),                                   &
               c(13,:,:,:),c(14,:,:,:),c(15,:,:,:),                      &
               c(16,:,:,:),c(17,:,:,:),c(18,:,:,:),t(1,:,:,:),t(2,:,:,:),&
               t(3,:,:,:),t(4,:,:,:),t(5,:,:,:),t(6,:,:,:),t(7,:,:,:),      &
               t(8,:,:,:),t(9,:,:,:),r(1,:,:,:),r(2,:,:,:),r(3,:,:,:),      &
               r(4,:,:,:), r(5,:,:,:), r(6,:,:,:),m,d,nl, h)

   call ricci4(c(1,:,:,:),c(2,:,:,:),c(3,:,:,:),c(4,:,:,:),c(5,:,:,:),      &
               c(6,:,:,:),c(7,:,:,:),c(8,:,:,:),c(9,:,:,:),c(10,:,:,:),     &
               c(11,:,:,:),c(12,:,:,:),                                   &
               c(13,:,:,:),c(14,:,:,:),c(15,:,:,:),                      &
               c(16,:,:,:),c(17,:,:,:),c(18,:,:,:),t(1,:,:,:),t(2,:,:,:),&
               t(3,:,:,:),t(4,:,:,:),t(5,:,:,:),t(6,:,:,:),t(7,:,:,:),      &
               t(8,:,:,:),t(9,:,:,:),r(1,:,:,:),r(2,:,:,:),r(3,:,:,:),      &
               r(4,:,:,:), r(5,:,:,:), r(6,:,:,:),m,d,nl, h)

   call ricci5(c(1,:,:,:),c(2,:,:,:),c(3,:,:,:),c(4,:,:,:),c(5,:,:,:),      &
               c(6,:,:,:),c(7,:,:,:),c(8,:,:,:),c(9,:,:,:),c(10,:,:,:),     &
               c(11,:,:,:),c(12,:,:,:),                                   &
               c(13,:,:,:),c(14,:,:,:),c(15,:,:,:),                      &
               c(16,:,:,:),c(17,:,:,:),c(18,:,:,:),t(1,:,:,:),t(2,:,:,:),&
               t(3,:,:,:),t(4,:,:,:),t(5,:,:,:),t(6,:,:,:),t(7,:,:,:),      &
               t(8,:,:,:),t(9,:,:,:),r(1,:,:,:),r(2,:,:,:),r(3,:,:,:),      &
               r(4,:,:,:), r(5,:,:,:), r(6,:,:,:),m,d,nl, h)

   call ricci6(c(1,:,:,:),c(2,:,:,:),c(3,:,:,:),c(4,:,:,:),c(5,:,:,:),      &
               c(6,:,:,:),c(7,:,:,:),c(8,:,:,:),c(9,:,:,:),c(10,:,:,:),     &
               c(11,:,:,:),c(12,:,:,:),                                   &
               c(13,:,:,:),c(14,:,:,:),c(15,:,:,:),                      &
               c(16,:,:,:),c(17,:,:,:),c(18,:,:,:),t(1,:,:,:),t(2,:,:,:),&
               t(3,:,:,:),t(4,:,:,:),t(5,:,:,:),t(6,:,:,:),t(7,:,:,:),      &
               t(8,:,:,:),t(9,:,:,:),r(1,:,:,:),r(2,:,:,:),r(3,:,:,:),      &
               r(4,:,:,:), r(5,:,:,:), r(6,:,:,:),m,d,nl, h)

   call ricci7(c(1,:,:,:),c(2,:,:,:),c(3,:,:,:),c(4,:,:,:),c(5,:,:,:),      &
               c(6,:,:,:),c(7,:,:,:),c(8,:,:,:),c(9,:,:,:),c(10,:,:,:),     &
               c(11,:,:,:),c(12,:,:,:),                                   &
               c(13,:,:,:),c(14,:,:,:),c(15,:,:,:),                      &
               c(16,:,:,:),c(17,:,:,:),c(18,:,:,:),t(1,:,:,:),t(2,:,:,:),&
               t(3,:,:,:),t(4,:,:,:),t(5,:,:,:),t(6,:,:,:),t(7,:,:,:),      &
               t(8,:,:,:),t(9,:,:,:),r(1,:,:,:),r(2,:,:,:),r(3,:,:,:),      &
               r(4,:,:,:), r(5,:,:,:), r(6,:,:,:),m,d,nl, h)

return 
end
!-----------------------------------------------------------------------
!  even though c{ijk}d{l} is symmetric in jk we keep other terms for
!  clean interfaces b/w cases where we have non-symmetric components later
!
!  routine checked to be second order in h

subroutine ricci1(c111, c112, c113, c122, c123, c133,        &
                  c211, c212, c213, c222, c223, c233,        &
                  c311, c312, c313, c322, c323, c333,        &
                  c111d1, c112d1, c113d1, c121d1, c122d1,    &
                  c123d1, c131d1, c132d1, c133d1,            &
                  r11, r12, r13, r22, r23, r33,              &
                  m,d,nl, h) 
implicit none
integer                       :: nl 
real*8                        :: h

real*8,dimension(:,:,:)    :: c111,c112,c113,c122,c123,c133
real*8,dimension(:,:,:)    :: c211,c212,c213,c222,c223,c233
real*8,dimension(:,:,:)    :: c311,c312,c313,c322,c323,c333
real*8,dimension(:,:,:)    :: c111d1,c112d1,c113d1
real*8,dimension(:,:,:)    :: c121d1,c122d1,c123d1
real*8,dimension(:,:,:)    :: c131d1,c132d1,c133d1

real*8,dimension(:,:,:)    :: r11,r12,r13,r22,r23,r33

real*8,dimension(:)        :: d
integer,dimension(:)       :: m
integer                       :: i,j,k
interface
   subroutine diff(h,f,df,dir,m,d,nl)
      implicit none
      integer                       :: nl
      real*8                        :: h
      real*8,dimension(:,:,:)    :: f, df
      real*8,dimension(:)        :: d
      integer,dimension(:)       :: m
      character                     :: dir
   end subroutine
end interface



!   Call difference routines for first term in Ricci
!   only calculate six derivatives 

   call diff(h, c111, c111d1, 'x',m,d,nl)
   call diff(h, c112, c112d1, 'x',m,d,nl)
   call diff(h, c113, c113d1, 'x',m,d,nl)
   call diff(h, c122, c122d1, 'x',m,d,nl)
   call diff(h, c123, c123d1, 'x',m,d,nl)
   call diff(h, c133, c133d1, 'x',m,d,nl)

!      Evaluate first term 
   r11 = c111d1
   r12 = c112d1
   r13 = c113d1
   r22 = c122d1
   r23 = c123d1
   r33 = c133d1

return
end
!-----------------------------------------------------------------------
!  routine checked to be second order in h
subroutine ricci2(c111, c112, c113, c122, c123, c133,        &
                  c211, c212, c213, c222, c223, c233,        &
                  c311, c312, c313, c322, c323, c333,        &
                  c211d2, c212d2, c213d2, c221d2, c222d2,    &
                  c223d2, c231d2, c232d2, c233d2,            &
                  r11, r12, r13, r22, r23, r33,              &
                  m,d,nl,h) 
implicit none
integer                       :: nl
real*8                        :: h
real*8,dimension(:,:,:)    :: c111,c112,c113,c122,c123,c133
real*8,dimension(:,:,:)    :: c211,c212,c213,c222,c223,c233
real*8,dimension(:,:,:)    :: c311,c312,c313,c322,c323,c333
real*8,dimension(:,:,:)    :: c211d2,c212d2,c213d2
real*8,dimension(:,:,:)    :: c221d2,c222d2,c223d2 
real*8,dimension(:,:,:)    :: c231d2,c232d2,c233d2

real*8,dimension(:,:,:)    :: r11,r12,r13,r22,r23,r33
real*8,dimension(:)        :: d
integer,dimension(:)       :: m
integer                       :: i,j,k

interface
   subroutine diff(h,f,df,dir,m,d,nl)
      implicit none
      integer                       :: nl
      real*8                        :: h
      real*8,dimension(:,:,:)    :: f, df
      real*8,dimension(:)        :: d
      integer,dimension(:)       :: m
      character                     :: dir
   end subroutine
end interface


!   Call difference routines for first term in Ricci
!   only calculate six derivatives 

   call diff(h, c211, c211d2, 'y',m,d,nl)
   call diff(h, c212, c212d2, 'y',m,d,nl)
   call diff(h, c213, c213d2, 'y',m,d,nl)
   call diff(h, c222, c222d2, 'y',m,d,nl)
   call diff(h, c223, c223d2, 'y',m,d,nl)
   call diff(h, c233, c233d2, 'y',m,d,nl)

!      Evaluate second term 
   r11 = r11 + c211d2
   r12 = r12 + c212d2
   r13 = r13 + c213d2
   r22 = r22 + c222d2
   r23 = r23 + c223d2
   r33 = r33 + c233d2

return
end
!-----------------------------------------------------------------------
!  routine checked to be second order in h
subroutine ricci3(c111, c112, c113, c122, c123, c133,         &
                  c211, c212, c213, c222, c223, c233,         &
                  c311, c312, c313, c322, c323, c333,         &
                  c311d3, c312d3, c313d3, c321d3, c322d3,     &
                  c323d3, c331d3, c332d3, c333d3,             &
                  r11, r12, r13, r22, r23, r33,               &
                  m,d,nl, h)
implicit none
integer                       :: nl
real*8                        :: h
real*8,dimension(:,:,:)    :: c111,c112,c113,c122,c123,c133
real*8,dimension(:,:,:)    :: c211,c212,c213,c222,c223,c233
real*8,dimension(:,:,:)    :: c311,c312,c313,c322,c323,c333
real*8,dimension(:,:,:)    :: c311d3,c312d3,c313d3 
real*8,dimension(:,:,:)    :: c321d3,c322d3,c323d3 
real*8,dimension(:,:,:)    :: c331d3,c332d3,c333d3 

real*8,dimension(:,:,:)    :: r11,r12,r13,r22,r23,r33
real*8,dimension(:)        :: d
integer,dimension(:)       :: m
integer                       :: i,j,k

interface
   subroutine diff(h,f,df,dir,m,d,nl)
      implicit none
      integer                       :: nl
      real*8                        :: h
      real*8,dimension(:,:,:)    :: f, df
      real*8,dimension(:)        :: d
      integer,dimension(:)       :: m
      character                     :: dir
   end subroutine
end interface



!   Call difference routines for first term in Ricci
!   only calculate six derivatives 

   call diff(h, c311, c311d3, 'z',m,d,nl)
   call diff(h, c312, c312d3, 'z',m,d,nl)
   call diff(h, c313, c313d3, 'z',m,d,nl)
   call diff(h, c322, c322d3, 'z',m,d,nl)
   call diff(h, c323, c323d3, 'z',m,d,nl)
   call diff(h, c333, c333d3, 'z',m,d,nl)

!      Evaluate third term 
   r11 = r11 + c311d3
   r12 = r12 + c312d3
   r13 = r13 + c313d3
   r22 = r22 + c322d3
   r23 = r23 + c323d3
   r33 = r33 + c333d3

return
end
!-----------------------------------------------------------------------
!  routine checked to be second order in h
subroutine ricci4(c111, c112, c113, c122, c123, c133,   &
                  c211, c212, c213, c222, c223, c233,   &
                  c311, c312, c313, c322, c323, c333,   &
                  c111d1, c111d2, c111d3, c121d1, c121d2,    &
                  c121d3, c131d1, c131d2, c131d3,   &
                  r11, r12, r13, r22, r23, r33,   &
                  m,d,nl, h)
implicit none
integer                       :: nl
real*8                        :: h
real*8,dimension(:,:,:)    :: c111,c112,c113,c122,c123,c133
real*8,dimension(:,:,:)    :: c211,c212,c213,c222,c223,c233
real*8,dimension(:,:,:)    :: c311,c312,c313,c322,c323,c333
real*8,dimension(:,:,:)    :: c111d1, c111d2, c111d3 
real*8,dimension(:,:,:)    :: c121d1, c121d2, c121d3 
real*8,dimension(:,:,:)    :: c131d1, c131d2, c131d3 

real*8,dimension(:,:,:)    :: r11,r12,r13,r22,r23,r33
real*8,dimension(:)        :: d
integer,dimension(:)       :: m
integer                       :: i,j,k

interface
   subroutine diff(h,f,df,dir,m,d,nl)
      implicit none
      integer                       :: nl
      real*8                        :: h
      real*8,dimension(:,:,:)    :: f, df
      real*8,dimension(:)        :: d
      integer,dimension(:)       :: m
      character                     :: dir
   end subroutine
end interface


!   Call difference routines for first term in Ricci
!   only calculate six derivatives 

   call diff(h, c111, c111d1, 'x',m,d,nl)
   call diff(h, c111, c111d2, 'y',m,d,nl)
   call diff(h, c111, c111d3, 'z',m,d,nl)
   call diff(h, c112, c121d2, 'y',m,d,nl)
   call diff(h, c112, c121d3, 'z',m,d,nl)
   call diff(h, c113, c131d3, 'z',m,d,nl)

!      Evaluate first term 
   r11 = r11 - c111d1
   r12 = r12 - c111d2
   r13 = r13 - c111d3
   r22 = r22 - c121d2
   r23 = r23 - c121d3
   r33 = r33 - c131d3

return
end
!-----------------------------------------------------------------------
!  routine checked to be second order in h
subroutine ricci5(c111, c112, c113, c122, c123, c133,         &
                  c211, c212, c213, c222, c223, c233,         &
                  c311, c312, c313, c322, c323, c333,         &
                  c212d1, c212d2, c212d3, c222d1, c222d2,     &
                  c222d3, c232d1, c232d2, c232d3,             &
                  r11, r12, r13, r22, r23, r33,               &
                  m,d,nl, h) 
implicit none
integer                       :: nl
real*8                        :: h
real*8,dimension(:,:,:)    :: c111,c112,c113,c122,c123,c133
real*8,dimension(:,:,:)    :: c211,c212,c213,c222,c223,c233
real*8,dimension(:,:,:)    :: c311,c312,c313,c322,c323,c333
real*8,dimension(:,:,:)    :: c212d1, c212d2, c212d3 
real*8,dimension(:,:,:)    :: c222d1, c222d2, c222d3 
real*8,dimension(:,:,:)    :: c232d1, c232d2, c232d3 

real*8,dimension(:,:,:)    :: r11,r12,r13,r22,r23,r33
real*8,dimension(:)        :: d
integer,dimension(:)       :: m
integer                       :: i,j,k
interface
   subroutine diff(h,f,df,dir,m,d,nl)
      implicit none
      integer                       :: nl
      real*8                        :: h
      real*8,dimension(:,:,:)    :: f, df
      real*8,dimension(:)        :: d
      integer,dimension(:)       :: m
      character                     :: dir
   end subroutine
end interface


!   Call difference routines for first term in Ricci
!   only calculate six derivatives 

   call diff(h, c212, c212d1, 'x',m,d,nl)
   call diff(h, c212, c212d2, 'y',m,d,nl)
   call diff(h, c212, c212d3, 'z',m,d,nl)
   call diff(h, c222, c222d2, 'y',m,d,nl)
   call diff(h, c222, c222d3, 'z',m,d,nl)
   call diff(h, c223, c232d3, 'z',m,d,nl)

!      Evaluate first term 
   r11 = r11 - c212d1
   r12 = r12 - c212d2
   r13 = r13 - c212d3
   r22 = r22 - c222d2
   r23 = r23 - c222d3
   r33 = r33 - c232d3

return
end
!-----------------------------------------------------------------------
subroutine ricci6(c111, c112, c113, c122, c123, c133,         &
                  c211, c212, c213, c222, c223, c233,         &
                  c311, c312, c313, c322, c323, c333,         &
                  c313d1, c313d2, c313d3, c323d1, c323d2,     &
                  c323d3, c333d1, c333d2, c333d3,             &
                  r11, r12, r13, r22, r23, r33,               &
                  m,d,nl, h) 
implicit none
integer                       :: nl 
real*8                        :: h
real*8,dimension(:,:,:)    :: c111,c112,c113,c122,c123,c133
real*8,dimension(:,:,:)    :: c211,c212,c213,c222,c223,c233
real*8,dimension(:,:,:)    :: c311,c312,c313,c322,c323,c333
real*8,dimension(:,:,:)    :: c313d1, c313d2, c313d3 
real*8,dimension(:,:,:)    :: c323d1, c323d2, c323d3 
real*8,dimension(:,:,:)    :: c333d1, c333d2, c333d3 

real*8,dimension(:,:,:)    :: r11,r12,r13,r22,r23,r33
real*8,dimension(:)        :: d
integer,dimension(:)       :: m
integer                       :: i,j,k

interface
   subroutine diff(h,f,df,dir,m,d,nl)
      implicit none
      integer                       :: nl
      real*8                        :: h
      real*8,dimension(:,:,:)    :: f, df
      real*8,dimension(:)        :: d
      integer,dimension(:)       :: m
      character                     :: dir
   end subroutine
end interface


!   Call difference routines for first term in Ricci
!   only calculate six derivatives 

   call diff(h, c313, c313d1, 'x',m,d,nl)
   call diff(h, c313, c313d2, 'y',m,d,nl)
   call diff(h, c313, c313d3, 'z',m,d,nl)
   call diff(h, c323, c323d2, 'y',m,d,nl)
   call diff(h, c323, c323d3, 'z',m,d,nl)
   call diff(h, c333, c333d3, 'z',m,d,nl)

!      Evaluate first term 
   r11 = r11 - c313d1
   r12 = r12 - c313d2
   r13 = r13 - c313d3
   r22 = r22 - c323d2
   r23 = r23 - c323d3
   r33 = r33 - c333d3

return
end
!-----------------------------------------------------------------------
!
! routine found to be correct to round off with exact connections
subroutine ricci7(c111, c112, c113, c122, c123, c133,    &
                  c211, c212, c213, c222, c223, c233,    &
                  c311, c312, c313, c322, c323, c333,    &
                  t1, t2, t3, t4, t5, t6, t7, t8, t9,    &
                  r11, r12, r13, r22, r23, r33,          &
                  m,d,nl, h) 
implicit none
integer                       :: nl
real*8                        :: h
real*8,dimension(:,:,:)    :: c111,c112,c113,c122,c123,c133
real*8,dimension(:,:,:)    :: c211,c212,c213,c222,c223,c233
real*8,dimension(:,:,:)    :: c311,c312,c313,c322,c323,c333
real*8,dimension(:,:,:)    :: t1,t2,t3,t4,t5,t6,t7,t8,t9 

real*8,dimension(:,:,:)    :: r11,r12,r13,r22,r23,r33
real*8,dimension(:)        :: d
integer,dimension(:)       :: m
integer                       :: i,j,k


!   Calculate third term
!     rij = rij + (c111 + c221 + c331)*c1ij +
!         (c112 + c222 + c332)*c2ij +
!         (c113 + c223 + c333)*c3ij

!           rij = rij + 
!            (c111 + c212 + c313)*c1ij +
!            (c112 + c222 + c323)*c2ij +
!            (c113 + c223 + c333)*c3ij

   r11 = r11 + (c111 + c212 + c313)*c111 +      &
               (c112 + c222 + c323)*c211 +      &
               (c113 + c223 + c333)*c311

   r12 = r12 + (c111 + c212 + c313)*c112 +      &
               (c112 + c222 + c323)*c212 +      &
               (c113 + c223 + c333)*c312

   r13 = r13 + (c111 + c212 + c313)*c113 +      &
               (c112 + c222 + c323)*c213 +      &
               (c113 + c223 + c333)*c313

   r22 = r22 + (c111 + c212 + c313)*c122 +      &
               (c112 + c222 + c323)*c222 +      &
               (c113 + c223 + c333)*c322

   r23 = r23 + (c111 + c212 + c313)*c123 +      &
               (c112 + c222 + c323)*c223 +      &
               (c113 + c223 + c333)*c323

   r33 = r33 + (c111 + c212 + c313)*c133 +      &
               (c112 + c222 + c323)*c233 +      &
               (c113 + c223 + c333)*c333

!   Calculate fourth term
!  
!   rij  = rij - ( c1i1  c11j + c2i1  c12j + c3i1*c13j)
!           - ( c1i2  c21j + c2i2  c22j + c3i2*c23j)
!           - ( c1i3  c31j + c2i3  c32j + c3i3*c33j)

!           rij  = rij - ( c1i1  c11j + 
!            c2i1  c12j + c3i1*c13j) -
!           ( c1i2  c21j + c2i2  c22j + 
!           c3i2*c23j) - ( c1i3  c31j + 
!           c2i3  c32j + c3i3*c33j)

   r11  = r11 - ( c111 * c111 +                    &
                  c211 * c112 + c311*c113) -       &
                ( c112 * c211 + c212 * c212 +      &
                  c312*c213) - ( c113 * c311 +     &
                  c213 * c312 + c313*c313)

   r12  = r12 - ( c111 * c112 +                    &
                  c211 * c122 + c311*c123) -       &
                ( c112 * c212 + c212 * c222 +      &
                  c312*c223) - ( c113 * c312 +     &
                  c213 * c322 + c313*c323)

   r13  = r13 - ( c111 * c113 +                    &
                  c211 * c123 + c311*c133) -       &
                ( c112 * c213 + c212 * c223 +      &
                  c312*c233) - ( c113 * c313 +     &
                  c213 * c323 + c313*c333)

   r22  = r22 - ( c112 * c112 +                    &
                  c212 * c122 + c312*c123) -       &
                ( c122 * c212 + c222 * c222 +      &
                  c322*c223) - ( c123 * c312 +     &
                  c223 * c322 + c323*c323)

   r23  = r23 - ( c112 * c113 +                    &
                  c212 * c123 + c312*c133) -       &
                ( c122 * c213 + c222 * c223 +      &
                  c322*c233) - ( c123 * c313 +     &
                  c223 * c323 + c323*c333)

   r33  = r33 - ( c113 * c113 +                    &
                  c213 * c123 + c313*c133) -       &
                ( c123 * c213 + c223 * c223 +      &
                  c323*c233) - ( c133 * c313 +     &
                  c233 * c323 + c333*c333)

return
end
!-----------------------------------------------------------------------
