!----------------------------------------------------------------------- 
! connect.f 
!----------------------------------------------------------------------- 
! This routine returns covariant form of ricci tensor. Raise index before 
! using in evolution code. 
! This routine assumes that the associated gamma will not be defined!
 
subroutine gamma(g,gup,h,t,nl,m,d,c,ca)
implicit none 
 
real*8                                        :: h 
real*8,dimension(:,:,:,:),       intent(in)   :: g,gup
real*8,dimension(:,:,:,:),       intent(out)  :: c,t,ca
integer                                       :: nl 
real*8,dimension(:)                        :: d
integer,dimension(:)                       :: m

! interfaces-----------------------------------------------------------Y
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

subroutine connect1(g11, g12, g13, g22, g23, g33,      &
                    gup11, gup21, gup31, gup22, gup23, gup33,      &
                    c111, c112, c113, c122, c123, c133,      &
                    c211, c212, c213, c222, c223, c233,      &
                    c311, c312, c313, c322, c323, c333,      &
                    g111, g112, g113, g121, g122, g123,       &
                    g131, g132, g133, h)
implicit none

real*8                        :: h
real*8,dimension(:,:,:),intent(IN)    :: g11,g12,g13,g22,g23,g33
real*8,dimension(:,:,:),intent(IN)    :: gup11,gup21,gup31,gup22,gup23,gup33

real*8,dimension(:,:,:),intent(OUT)    :: c111,c112,c113,c122,c123,c133
real*8,dimension(:,:,:),intent(OUT)    :: c211,c212,c213,c222,c223,c233
real*8,dimension(:,:,:),intent(OUT)    :: c311,c312,c313,c322,c323,c333

real*8,dimension(:,:,:),intent(IN)     :: g111,g112,g113,g121, &
                                          g122,g123,g131,g132,g133
end subroutine

subroutine connect2(g11, g12, g13, g22, g23, g33,    &
                    gup11, gup12, gup31, gup22, gup32, gup33,    &
                    c111, c112, c113, c122, c123, c133,    &
                    c211, c212, c213, c222, c223, c233,    &
                    c311, c312, c313, c322, c323, c333,    &
                    g211, g212, g213, g222, g223, g233,     &
                    g111, g112, g333, h)
implicit none

real*8                        :: h
real*8,dimension(:,:,:),intent(IN)    :: g11,g12,g13,g22,g23,g33
real*8,dimension(:,:,:),intent(IN)    :: gup11,gup12,gup31,gup22,gup32,gup33

real*8,dimension(:,:,:),intent(OUT)    :: c111,c112,c113,c122,c123,c133
real*8,dimension(:,:,:),intent(OUT)    :: c211,c212,c213,c222,c223,c233
real*8,dimension(:,:,:),intent(OUT)    :: c311,c312,c313,c322,c323,c333

real*8,dimension(:,:,:),intent(IN)     :: g211,g212,g213,g222,g223,    &
                                          g233,g111,g112,g333

end subroutine

subroutine connect3(g11, g12, g13, g22, g23, g33,      &
                    gup11, gup12, gup13, gup22, gup23, gup33,      &
                    c111, c112, c113, c122, c123, c133,      &
                    c211, c212, c213, c222, c223, c233,      &
                    c311, c312, c313, c322, c323, c333,      &
                    g311, g312, g313, g322, g323, g333,       &
                    g332, g112, g233, h)

implicit none

real*8                        :: h
real*8,dimension(:,:,:),intent(IN)    :: g11,g12,g13,g22,g23,g33
real*8,dimension(:,:,:),intent(IN)    :: gup11,gup12,gup13,gup22,gup23,gup33

real*8,dimension(:,:,:),intent(out)    :: c111,c112,c113,c122,c123,c133
real*8,dimension(:,:,:),intent(out)    :: c211,c212,c213,c222,c223,c233
real*8,dimension(:,:,:),intent(out)    :: c311,c312,c313,c322,c323,c333

real*8,dimension(:,:,:),intent(in)     :: g311,g312,g313,g322,      &
                                          g323,g333,g332,g112,g233
end subroutine

subroutine connect4(g11, g12, g13, g22, g23, g33,      &
                    gup11, gup21, gup31, gup22, gup32, gup33,      &
                    c111, c112, c113, c122, c123, c133,      &
                    c211, c212, c213, c222, c223, c233,      &
                    c311, c312, c313, c322, c323, c333,      &
                    g111, g121, g131, g122, g132, g133,       &
                    g222, g333, g311, h)
implicit none
real*8                        :: h
real*8,dimension(:,:,:),intent(in)    :: g11,g12,g13,g22,g23,g33
real*8,dimension(:,:,:),intent(in)    :: gup11,gup21,gup31,gup22,gup32,gup33

real*8,dimension(:,:,:),intent(OUT)    :: c111,c112,c113,c122,c123,c133
real*8,dimension(:,:,:),intent(OUT)    :: c211,c212,c213,c222,c223,c233
real*8,dimension(:,:,:),intent(OUT)    :: c311,c312,c313,c322,c323,c333

real*8,dimension(:,:,:),intent(in)    :: g111,g121,g131,g122, &
                                         g132,g133,g222,g333,g311
end subroutine

subroutine connect5(g11, g12, g13, g22, g23, g33,      &
                    gup11, gup12, gup31, gup22, gup32, gup33,      &
                    c111, c112, c113, c122, c123, c133,      &
                    c211, c212, c213, c222, c223, c233,      &
                    c311, c312, c313, c322, c323, c333,      &
                    g211, g221, g231, g222, g232, g233,       &
                    g131, g132, g133, h)

implicit none
real*8                        :: h
real*8,dimension(:,:,:),intent(in)    :: g11,g12,g13,g22,g23,g33
real*8,dimension(:,:,:),intent(in)    :: gup11,gup12,gup31,gup22,gup32,gup33

real*8,dimension(:,:,:),intent(out)    :: c111,c112,c113,c122,c123,c133
real*8,dimension(:,:,:),intent(out)    :: c211,c212,c213,c222,c223,c233
real*8,dimension(:,:,:),intent(out)    :: c311,c312,c313,c322,c323,c333
real*8,dimension(:,:,:),intent(in)     :: g211,g221,g231,g222,      &
                                          g232,g233,g131,g132,g133
end subroutine

subroutine connect6(g11, g12, g13, g22, g23, g33,      &
                    gup11, gup12, gup13, gup22, gup23, gup33,      &
                    c111, c112, c113, c122, c123, c133,      &
                    c211, c212, c213, c222, c223, c233,      &
                    c311, c312, c313, c322, c323, c333,      &
                    g311, g321, g331, g322, g332, g333,       &
                    g131, g132, g133, h)

implicit none
real*8                        :: h
real*8,dimension(:,:,:),intent(in)    :: g11,g12,g13,g22,g23,g33
real*8,dimension(:,:,:),intent(in)    :: gup11,gup12,gup13,gup22,gup23,gup33

real*8,dimension(:,:,:),intent(out)    :: c111,c112,c113,c122,c123,c133
real*8,dimension(:,:,:),intent(out)    :: c211,c212,c213,c222,c223,c233
real*8,dimension(:,:,:),intent(out)    :: c311,c312,c313,c322,c323,c333

real*8,dimension(:,:,:),intent(in)    :: g311,g321,g331,g322,        &
                                         g332,g333,g131,g132,g133
end subroutine

subroutine connect7(g11, g12, g13, g22, g23, g33,      &
                    gup11, gup21, gup31, gup22, gup32, gup33,      &
                    c111, c112, c113, c122, c123, c133,      &
                    c211, c212, c213, c222, c223, c233,      &
                    c311, c312, c313, c322, c323, c333,      &
                    g111, g121, g131, g221, g231, g331,       &
                    g333, g132, g133, h)

implicit none
real*8                        :: h
real*8,dimension(:,:,:),intent(in)    :: g11,g12,g13,g22,g23,g33
real*8,dimension(:,:,:),intent(in)    :: gup11,gup21,gup31,gup22,gup32,gup33

real*8,dimension(:,:,:),intent(out)    :: c111,c112,c113,c122,c123,c133
real*8,dimension(:,:,:),intent(out)    :: c211,c212,c213,c222,c223,c233
real*8,dimension(:,:,:),intent(out)    :: c311,c312,c313,c322,c323,c333

real*8,dimension(:,:,:),intent(in)    :: g111,g121,g131,g221,g231,  &
                                         g331,g333,g132,g133
end subroutine

subroutine connect8(g11, g12, g13, g22, g23, g33,      &
                    gup11, gup12, gup31, gup22, gup32, gup33,      &
                    c111, c112, c113, c122, c123, c133,      &
                    c211, c212, c213, c222, c223, c233,      &
                    c311, c312, c313, c322, c323, c333,      &
                    g112, g122, g132, g222, g232, g332,       &
                    g131, g333, g133, h)

implicit none
real*8                        :: h
real*8,dimension(:,:,:),intent(IN)    :: g11,g12,g13,g22,g23,g33
real*8,dimension(:,:,:),intent(IN)    :: gup11,gup12,gup31,gup22,gup32,gup33

real*8,dimension(:,:,:),intent(OUT)    :: c111,c112,c113,c122,c123,c133
real*8,dimension(:,:,:),intent(OUT)    :: c211,c212,c213,c222,c223,c233
real*8,dimension(:,:,:),intent(OUT)    :: c311,c312,c313,c322,c323,c333

real*8,dimension(:,:,:),intent(IN)  :: g112,g122,g132,g222,     &
                                       g232,g332,g131,g333,g133
end subroutine

subroutine connect9(g11, g12, g13, g22, g23, g33,      &
                    gup11, gup12, gup13, gup22, gup23, gup33,      &
                    c111, c112, c113, c122, c123, c133,      &
                    c211, c212, c213, c222, c223, c233,      &
                    c311, c312, c313, c322, c323, c333,      &
                    g113, g123, g133, g223, g233, g333,       &
                    g131, g132, g222, h)

implicit none
real*8                        :: h
real*8,dimension(:,:,:),intent(IN)    :: g11,g12,g13,g22,g23,g33
real*8,dimension(:,:,:),intent(IN)    :: gup11,gup12,gup13,gup22,gup23,gup33

real*8,dimension(:,:,:),intent(OUT)   :: c111,c112,c113,c122,c123,c133
real*8,dimension(:,:,:),intent(OUT)    :: c211,c212,c213,c222,c223,c233
real*8,dimension(:,:,:),intent(OUT)    :: c311,c312,c313,c322,c323,c333
real*8,dimension(:,:,:),intent(IN)    :: g113,g123,g133,g223,   &
                                         g233,g333,g131,g132,g222
end subroutine

end interface
! interfaces-----------------------------------------------------------Y

!--------------------------------------------------------------------
!  call subsidiary routines to calculate  Christoffell tensor piecewise. 
 
   c = 0.0d0
   ca= 0.0d0
 
   call diff(h,g(1,:,:,:),t(1,:,:,:),'x',m,d,nl) 
   call diff(h,g(1,:,:,:),t(2,:,:,:),'y',m,d,nl) 
   call diff(h,g(1,:,:,:),t(3,:,:,:),'z',m,d,nl) 
 
   call diff(h,g(2,:,:,:),t(4,:,:,:),'x',m,d,nl) 
   call diff(h,g(2,:,:,:),t(5,:,:,:),'y',m,d,nl) 
   call diff(h,g(2,:,:,:),t(6,:,:,:),'z',m,d,nl) 
 
   call diff(h,g(3,:,:,:),t(7,:,:,:),'x',m,d,nl) 
   call diff(h,g(3,:,:,:),t(8,:,:,:),'y',m,d,nl) 
   call diff(h,g(3,:,:,:),t(9,:,:,:),'z',m,d,nl) 
 
! note that the first 6 t's are really used, the last 3 are not.... 
 
   call connect1(g(1,:,:,:), g(2,:,:,:), g(3,:,:,:),                     &
                 g(4,:,:,:), g(5,:,:,:), g(6,:,:,:),                     &
                 gup(1,:,:,:), gup(2,:,:,:), gup(3,:,:,:),               &
                 gup(4,:,:,:), gup(5,:,:,:), gup(6,:,:,:),               &
                 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,:,:,:), h) 
 
   call connect4(g(1,:,:,:), g(2,:,:,:), g(3,:,:,:),                     &
                 g(4,:,:,:), g(5,:,:,:), g(6,:,:,:),                     &
                 gup(1,:,:,:), gup(2,:,:,:), gup(3,:,:,:),               &
                 gup(4,:,:,:), gup(5,:,:,:), gup(6,:,:,:),               &
                 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(4,:,:,:),t(7,:,:,:),                       &
                 t(5,:,:,:),t(8,:,:,:),t(9,:,:,:),                       &
                 t(2,:,:,:), t(3,:,:,:), t(6,:,:,:), h)

 
   call diff(h,g(2,:,:,:),t(1,:,:,:),'x',m,d,nl) 
   call diff(h,g(2,:,:,:),t(2,:,:,:),'y',m,d,nl) 
   call diff(h,g(2,:,:,:),t(3,:,:,:),'z',m,d,nl) 
 
   call diff(h,g(4,:,:,:),t(4,:,:,:),'x',m,d,nl) 
   call diff(h,g(4,:,:,:),t(5,:,:,:),'y',m,d,nl) 
   call diff(h,g(4,:,:,:),t(6,:,:,:),'z',m,d,nl) 
 
   call diff(h,g(5,:,:,:),t(7,:,:,:),'x',m,d,nl) 
   call diff(h,g(5,:,:,:),t(8,:,:,:),'y',m,d,nl) 
   call diff(h,g(5,:,:,:),t(9,:,:,:),'z',m,d,nl) 

   call connect2(g(1,:,:,:), g(2,:,:,:), g(3,:,:,:),                     &
                 g(4,:,:,:), g(5,:,:,:), g(6,:,:,:),                     &
                 gup(1,:,:,:), gup(2,:,:,:), gup(3,:,:,:),               &
                 gup(4,:,:,:), gup(5,:,:,:), gup(6,:,:,:),               &
                 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(5,:,:,:),t(6,:,:,:),t(9,:,:,:),                       &
                 t(7,:,:,:), t(8,:,:,:), t(6,:,:,:), h)

   call connect5(g(1,:,:,:), g(2,:,:,:), g(3,:,:,:),                     &
                 g(4,:,:,:), g(5,:,:,:), g(6,:,:,:),                     &
                 gup(1,:,:,:), gup(2,:,:,:), gup(3,:,:,:),               &
                 gup(4,:,:,:), gup(5,:,:,:), gup(6,:,:,:),               &
                 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(4,:,:,:),t(7,:,:,:),                       &
                 t(5,:,:,:),t(8,:,:,:),t(9,:,:,:),                       &
                 t(2,:,:,:), t(3,:,:,:), t(6,:,:,:), h)

 
 
   call diff(h,g(3,:,:,:),t(1,:,:,:),'x',m,d,nl) 
   call diff(h,g(3,:,:,:),t(2,:,:,:),'y',m,d,nl) 
   call diff(h,g(3,:,:,:),t(3,:,:,:),'z',m,d,nl) 
 
   call diff(h,g(5,:,:,:),t(4,:,:,:),'x',m,d,nl) 
   call diff(h,g(5,:,:,:),t(5,:,:,:),'y',m,d,nl) 
   call diff(h,g(5,:,:,:),t(6,:,:,:),'z',m,d,nl) 
 
   call diff(h,g(6,:,:,:),t(7,:,:,:),'x',m,d,nl) 
   call diff(h,g(6,:,:,:),t(8,:,:,:),'y',m,d,nl) 
   call diff(h,g(6,:,:,:),t(9,:,:,:),'z',m,d,nl) 
 

   call connect3(g(1,:,:,:), g(2,:,:,:), g(3,:,:,:),                     &
                 g(4,:,:,:), g(5,:,:,:), g(6,:,:,:),                     &
                 gup(1,:,:,:), gup(2,:,:,:), gup(3,:,:,:),               &
                 gup(4,:,:,:), gup(5,:,:,:), gup(6,:,:,:),               &
                 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(5,:,:,:),t(6,:,:,:),t(9,:,:,:),                       &
                 t(7,:,:,:), t(8,:,:,:), t(4,:,:,:), h)

 
   call connect6(g(1,:,:,:), g(2,:,:,:), g(3,:,:,:),                     &
                 g(4,:,:,:), g(5,:,:,:), g(6,:,:,:),                     &
                 gup(1,:,:,:), gup(2,:,:,:), gup(3,:,:,:),               &
                 gup(4,:,:,:), gup(5,:,:,:), gup(6,:,:,:),               &
                 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(4,:,:,:),t(7,:,:,:),                       &
                 t(5,:,:,:),t(8,:,:,:),t(9,:,:,:),                       &
                 t(2,:,:,:), t(3,:,:,:), t(6,:,:,:), h)
 
! now we calculate the last part ... 
 
   call diff(h,g(1,:,:,:),t(1,:,:,:),'x',m,d,nl) 
   call diff(h,g(2,:,:,:),t(2,:,:,:),'x',m,d,nl) 
   call diff(h,g(3,:,:,:),t(3,:,:,:),'x',m,d,nl) 

   call diff(h,g(4,:,:,:),t(4,:,:,:),'x',m,d,nl) 
   call diff(h,g(5,:,:,:),t(5,:,:,:),'x',m,d,nl) 
   call diff(h,g(6,:,:,:),t(6,:,:,:),'x',m,d,nl) 
 
   call connect7(g(1,:,:,:), g(2,:,:,:), g(3,:,:,:),                     &
                 g(4,:,:,:), g(5,:,:,:), g(6,:,:,:),                     &
                 gup(1,:,:,:), gup(2,:,:,:), gup(3,:,:,:),               &
                 gup(4,:,:,:), gup(5,:,:,:), gup(6,:,:,:),               &
                 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,:,:,:), h)
 
   call diff(h,g(1,:,:,:),t(1,:,:,:),'y',m,d,nl) 
   call diff(h,g(2,:,:,:),t(2,:,:,:),'y',m,d,nl) 
   call diff(h,g(3,:,:,:),t(3,:,:,:),'y',m,d,nl) 

   call diff(h,g(4,:,:,:),t(4,:,:,:),'y',m,d,nl) 
   call diff(h,g(5,:,:,:),t(5,:,:,:),'y',m,d,nl) 
   call diff(h,g(6,:,:,:),t(6,:,:,:),'y',m,d,nl) 

   call connect8(g(1,:,:,:), g(2,:,:,:), g(3,:,:,:),                     &
                 g(4,:,:,:), g(5,:,:,:), g(6,:,:,:),                     &
                 gup(1,:,:,:), gup(2,:,:,:), gup(3,:,:,:),               &
                 gup(4,:,:,:), gup(5,:,:,:), gup(6,:,:,:),               &
                 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,:,:,:), h)
 
   call diff(h,g(1,:,:,:),t(1,:,:,:),'z',m,d,nl) 
   call diff(h,g(2,:,:,:),t(2,:,:,:),'z',m,d,nl) 
   call diff(h,g(3,:,:,:),t(3,:,:,:),'z',m,d,nl) 

   call diff(h,g(4,:,:,:),t(4,:,:,:),'z',m,d,nl) 
   call diff(h,g(5,:,:,:),t(5,:,:,:),'z',m,d,nl) 
   call diff(h,g(6,:,:,:),t(6,:,:,:),'z',m,d,nl) 
 
   call connect9(g(1,:,:,:), g(2,:,:,:), g(3,:,:,:),                     &
                 g(4,:,:,:), g(5,:,:,:), g(6,:,:,:),                     &
                 gup(1,:,:,:), gup(2,:,:,:), gup(3,:,:,:),               &
                 gup(4,:,:,:), gup(5,:,:,:), gup(6,:,:,:),               &
                 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,:,:,:), h)


return  
end 

!----------------------------------------------------------------------- 
! note: we accept in gup21 and gup31 = gup12 , gup13 
subroutine connect1(g11, g12, g13, g22, g23, g33,      &
                    gup11, gup21, gup31, gup22, gup23, gup33,      &
                    c111, c112, c113, c122, c123, c133,      &
                    c211, c212, c213, c222, c223, c233,      &
                    c311, c312, c313, c322, c323, c333,      &
                    g111, g112, g113, g121, g122, g123,       &
                    g131, g132, g133, h) 
implicit none

real*8                        :: h
real*8,dimension(:,:,:),intent(IN)    :: g11,g12,g13,g22,g23,g33
real*8,dimension(:,:,:),intent(IN)    :: gup11,gup21,gup31,gup22,gup23,gup33

real*8,dimension(:,:,:),intent(OUT)    :: c111,c112,c113,c122,c123,c133
real*8,dimension(:,:,:),intent(OUT)    :: c211,c212,c213,c222,c223,c233
real*8,dimension(:,:,:),intent(OUT)    :: c311,c312,c313,c322,c323,c333

real*8,dimension(:,:,:),intent(IN)     :: g111,g112,g113,g121, &
                                          g122,g123,g131,g132,g133 

   c111 = c111  + 0.5d0*gup11 *g111  
   c112 = c112  + 0.5d0*gup11 *g112  
   c113 = c113  + 0.5d0*gup11 *g113  
   c122 = c122  + 0.5d0*gup11 *g122  
   c123 = c123  + 0.5d0*gup11 *g123  
   c133 = c133  + 0.5d0*gup11 *g133  
 
   c211 = c211  + 0.5d0*gup21 *g111  
   c212 = c212  + 0.5d0*gup21 *g112  
   c213 = c213  + 0.5d0*gup21 *g113  
   c222 = c222  + 0.5d0*gup21 *g122  
   c223 = c223  + 0.5d0*gup21 *g123  
   c233 = c233  + 0.5d0*gup21 *g133  

   c311 = c311  + 0.5d0*gup31 *g111  
   c312 = c312  + 0.5d0*gup31 *g112  
   c313 = c313  + 0.5d0*gup31 *g113  
   c322 = c322  + 0.5d0*gup31 *g122  
   c323 = c323  + 0.5d0*gup31 *g123  
   c333 = c333  + 0.5d0*gup31 *g133  
 
 
return 
end 
!----------------------------------------------------------------------- 
! note gup23 -> gup32 
subroutine connect2(g11, g12, g13, g22, g23, g33,    &
                    gup11, gup12, gup31, gup22, gup32, gup33,    &
                    c111, c112, c113, c122, c123, c133,    &
                    c211, c212, c213, c222, c223, c233,    &
                    c311, c312, c313, c322, c323, c333,    &
                    g211, g212, g213, g222, g223, g233,     &
                    g111, g112, g333, h) 
implicit none

real*8                        :: h
real*8,dimension(:,:,:),intent(IN)    :: g11,g12,g13,g22,g23,g33
real*8,dimension(:,:,:),intent(IN)    :: gup11,gup12,gup31,gup22,gup32,gup33

real*8,dimension(:,:,:),intent(OUT)    :: c111,c112,c113,c122,c123,c133
real*8,dimension(:,:,:),intent(OUT)    :: c211,c212,c213,c222,c223,c233
real*8,dimension(:,:,:),intent(OUT)    :: c311,c312,c313,c322,c323,c333

real*8,dimension(:,:,:),intent(IN)     :: g211,g212,g213,g222,g223,    &
                                          g233,g111,g112,g333

 
   c111 = c111  + 0.5d0*gup12 *g211  
   c112 = c112  + 0.5d0*gup12 *g212  
   c113 = c113  + 0.5d0*gup12 *g213  
   c122 = c122  + 0.5d0*gup12 *g222  
   c123 = c123  + 0.5d0*gup12 *g223  
   c133 = c133  + 0.5d0*gup12 *g233  
 
   c211 = c211  + 0.5d0*gup22 *g211  
   c212 = c212  + 0.5d0*gup22 *g212  
   c213 = c213  + 0.5d0*gup22 *g213  
   c222 = c222  + 0.5d0*gup22 *g222  
   c223 = c223  + 0.5d0*gup22 *g223  
   c233 = c233  + 0.5d0*gup22 *g233  
 
   c311 = c311  + 0.5d0*gup32 *g211  
   c312 = c312  + 0.5d0*gup32 *g212  
   c313 = c313  + 0.5d0*gup32 *g213  
   c322 = c322  + 0.5d0*gup32 *g222  
   c323 = c323  + 0.5d0*gup32 *g223  
   c333 = c333  + 0.5d0*gup32 *g233  
 
 
return 
end 
 
!----------------------------------------------------------------------- 
subroutine connect3(g11, g12, g13, g22, g23, g33,      &
                    gup11, gup12, gup13, gup22, gup23, gup33,      &
                    c111, c112, c113, c122, c123, c133,      &
                    c211, c212, c213, c222, c223, c233,      &
                    c311, c312, c313, c322, c323, c333,      &
                    g311, g312, g313, g322, g323, g333,       &
                    g332, g112, g233, h)

implicit none

real*8                        :: h
real*8,dimension(:,:,:),intent(IN)    :: g11,g12,g13,g22,g23,g33
real*8,dimension(:,:,:),intent(IN)    :: gup11,gup12,gup13,gup22,gup23,gup33

real*8,dimension(:,:,:),intent(out)    :: c111,c112,c113,c122,c123,c133
real*8,dimension(:,:,:),intent(out)    :: c211,c212,c213,c222,c223,c233
real*8,dimension(:,:,:),intent(out)    :: c311,c312,c313,c322,c323,c333

real*8,dimension(:,:,:),intent(in)     :: g311,g312,g313,g322,      &
                                          g323,g333,g332,g112,g233 


   c111 = c111  + 0.5d0*gup13 *g311  
   c112 = c112  + 0.5d0*gup13 *g312  
   c113 = c113  + 0.5d0*gup13 *g313  
   c122 = c122  + 0.5d0*gup13 *g322  
   c123 = c123  + 0.5d0*gup13 *g323  
   c133 = c133  + 0.5d0*gup13 *g333  
 
   c211 = c211  + 0.5d0*gup23 *g311  
   c212 = c212  + 0.5d0*gup23 *g312  
   c213 = c213  + 0.5d0*gup23 *g313  
   c222 = c222  + 0.5d0*gup23 *g322  
   c223 = c223  + 0.5d0*gup23 *g323  
   c233 = c233  + 0.5d0*gup23 *g333  
 
   c311 = c311  + 0.5d0*gup33 *g311  
   c312 = c312  + 0.5d0*gup33 *g312  
   c313 = c313  + 0.5d0*gup33 *g313  
   c322 = c322  + 0.5d0*gup33 *g322  
   c323 = c323  + 0.5d0*gup33 *g323  
   c333 = c333  + 0.5d0*gup33 *g333  
 
 
return 
end
!----------------------------------------------------------------------- 
      subroutine connect4(g11, g12, g13, g22, g23, g33,      &
                          gup11, gup21, gup31, gup22, gup32, gup33,      &
                          c111, c112, c113, c122, c123, c133,      &
                          c211, c212, c213, c222, c223, c233,      &
                          c311, c312, c313, c322, c323, c333,      &
                          g111, g121, g131, g122, g132, g133,       &
                          g222, g333, g311, h)
implicit none
real*8                        :: h
real*8,dimension(:,:,:),intent(in)    :: g11,g12,g13,g22,g23,g33
real*8,dimension(:,:,:),intent(in)    :: gup11,gup21,gup31,gup22,gup32,gup33

real*8,dimension(:,:,:),intent(OUT)    :: c111,c112,c113,c122,c123,c133
real*8,dimension(:,:,:),intent(OUT)    :: c211,c212,c213,c222,c223,c233
real*8,dimension(:,:,:),intent(OUT)    :: c311,c312,c313,c322,c323,c333

real*8,dimension(:,:,:),intent(in)    :: g111,g121,g131,g122, &
                                         g132,g133,g222,g333,g311
 
   c111 = c111  + 0.5d0*gup11 *g111  
   c112 = c112  + 0.5d0*gup11 *g121  
   c113 = c113  + 0.5d0*gup11 *g131  
   c122 = c122  + 0.5d0*gup11 *g122  
   c123 = c123  + 0.5d0*gup11 *g132  
   c133 = c133  + 0.5d0*gup11 *g133  
 
   c211 = c211  + 0.5d0*gup21 *g111  
   c212 = c212  + 0.5d0*gup21 *g121  
   c213 = c213  + 0.5d0*gup21 *g131  
   c222 = c222  + 0.5d0*gup21 *g122  
   c223 = c223  + 0.5d0*gup21 *g132  
   c233 = c233  + 0.5d0*gup21 *g133  
 
   c311 = c311  + 0.5d0*gup31 *g111  
   c312 = c312  + 0.5d0*gup31 *g121  
   c313 = c313  + 0.5d0*gup31 *g131  
   c322 = c322  + 0.5d0*gup31 *g122  
   c323 = c323  + 0.5d0*gup31 *g132  
   c333 = c333  + 0.5d0*gup31 *g133  
 
 
return 
end 
!----------------------------------------------------------------------- 
subroutine connect5(g11, g12, g13, g22, g23, g33,      &
                    gup11, gup12, gup31, gup22, gup32, gup33,      &
                    c111, c112, c113, c122, c123, c133,      &
                    c211, c212, c213, c222, c223, c233,      &
                    c311, c312, c313, c322, c323, c333,      &
                    g211, g221, g231, g222, g232, g233,       &
                    g131, g132, g133, h)

implicit none
real*8                        :: h
real*8,dimension(:,:,:),intent(in)    :: g11,g12,g13,g22,g23,g33
real*8,dimension(:,:,:),intent(in)    :: gup11,gup12,gup31,gup22,gup32,gup33

real*8,dimension(:,:,:),intent(out)    :: c111,c112,c113,c122,c123,c133
real*8,dimension(:,:,:),intent(out)    :: c211,c212,c213,c222,c223,c233
real*8,dimension(:,:,:),intent(out)    :: c311,c312,c313,c322,c323,c333
real*8,dimension(:,:,:),intent(in)     :: g211,g221,g231,g222,      &
                                          g232,g233,g131,g132,g133 
 
   c111 = c111  + 0.5d0*gup12 *g211  
   c112 = c112  + 0.5d0*gup12 *g221  
   c113 = c113  + 0.5d0*gup12 *g231  
   c122 = c122  + 0.5d0*gup12 *g222  
   c123 = c123  + 0.5d0*gup12 *g232  
   c133 = c133  + 0.5d0*gup12 *g233  
 
   c211 = c211  + 0.5d0*gup22 *g211  
   c212 = c212  + 0.5d0*gup22 *g221  
   c213 = c213  + 0.5d0*gup22 *g231  
   c222 = c222  + 0.5d0*gup22 *g222  
   c223 = c223  + 0.5d0*gup22 *g232  
   c233 = c233  + 0.5d0*gup22 *g233  
 
   c311 = c311  + 0.5d0*gup32 *g211  
   c312 = c312  + 0.5d0*gup32 *g221  
   c313 = c313  + 0.5d0*gup32 *g231  
   c322 = c322  + 0.5d0*gup32 *g222  
   c323 = c323  + 0.5d0*gup32 *g232  
   c333 = c333  + 0.5d0*gup32 *g233  
 
 
return 
end 
!----------------------------------------------------------------------- 
subroutine connect6(g11, g12, g13, g22, g23, g33,      &
                    gup11, gup12, gup13, gup22, gup23, gup33,      &
                    c111, c112, c113, c122, c123, c133,      &
                    c211, c212, c213, c222, c223, c233,      &
                    c311, c312, c313, c322, c323, c333,      &
                    g311, g321, g331, g322, g332, g333,       &
                    g131, g132, g133, h) 

implicit none
real*8                        :: h
real*8,dimension(:,:,:),intent(in)    :: g11,g12,g13,g22,g23,g33
real*8,dimension(:,:,:),intent(in)    :: gup11,gup12,gup13,gup22,gup23,gup33

real*8,dimension(:,:,:),intent(out)    :: c111,c112,c113,c122,c123,c133
real*8,dimension(:,:,:),intent(out)    :: c211,c212,c213,c222,c223,c233
real*8,dimension(:,:,:),intent(out)    :: c311,c312,c313,c322,c323,c333

real*8,dimension(:,:,:),intent(in)    :: g311,g321,g331,g322,        &
                                         g332,g333,g131,g132,g133  
 
   c111 = c111  + 0.5d0*gup13 *g311  
   c112 = c112  + 0.5d0*gup13 *g321  
   c113 = c113  + 0.5d0*gup13 *g331  
   c122 = c122  + 0.5d0*gup13 *g322  
   c123 = c123  + 0.5d0*gup13 *g332  
   c133 = c133  + 0.5d0*gup13 *g333  

   c211 = c211  + 0.5d0*gup23 *g311  
   c212 = c212  + 0.5d0*gup23 *g321  
   c213 = c213  + 0.5d0*gup23 *g331  
   c222 = c222  + 0.5d0*gup23 *g322  
   c223 = c223  + 0.5d0*gup23 *g332  
   c233 = c233  + 0.5d0*gup23 *g333  
 
   c311 = c311  + 0.5d0*gup33 *g311  
   c312 = c312  + 0.5d0*gup33 *g321  
   c313 = c313  + 0.5d0*gup33 *g331  
   c322 = c322  + 0.5d0*gup33 *g322  
   c323 = c323  + 0.5d0*gup33 *g332  
   c333 = c333  + 0.5d0*gup33 *g333  

 
return 
end 
!----------------------------------------------------------------------- 
subroutine connect7(g11, g12, g13, g22, g23, g33,      &
                    gup11, gup21, gup31, gup22, gup32, gup33,      &
                    c111, c112, c113, c122, c123, c133,      &
                    c211, c212, c213, c222, c223, c233,      &
                    c311, c312, c313, c322, c323, c333,      &
                    g111, g121, g131, g221, g231, g331,       &
                    g333, g132, g133, h) 
 
implicit none
real*8                        :: h
real*8,dimension(:,:,:),intent(in)    :: g11,g12,g13,g22,g23,g33
real*8,dimension(:,:,:),intent(in)    :: gup11,gup21,gup31,gup22,gup32,gup33

real*8,dimension(:,:,:),intent(out)    :: c111,c112,c113,c122,c123,c133
real*8,dimension(:,:,:),intent(out)    :: c211,c212,c213,c222,c223,c233
real*8,dimension(:,:,:),intent(out)    :: c311,c312,c313,c322,c323,c333

real*8,dimension(:,:,:),intent(in)    :: g111,g121,g131,g221,g231,  &
                                         g331,g333,g132,g133 
 
   c111 = c111  - 0.5d0*gup11 *g111  
   c112 = c112  - 0.5d0*gup11 *g121  
   c113 = c113  - 0.5d0*gup11 *g131  
   c122 = c122  - 0.5d0*gup11 *g221  
   c123 = c123  - 0.5d0*gup11 *g231  
   c133 = c133  - 0.5d0*gup11 *g331  

   c211 = c211  - 0.5d0*gup21 *g111  
   c212 = c212  - 0.5d0*gup21 *g121  
   c213 = c213  - 0.5d0*gup21 *g131  
   c222 = c222  - 0.5d0*gup21 *g221  
   c223 = c223  - 0.5d0*gup21 *g231  
   c233 = c233  - 0.5d0*gup21 *g331  
 
   c311 = c311  - 0.5d0*gup31 *g111  
   c312 = c312  - 0.5d0*gup31 *g121  
   c313 = c313  - 0.5d0*gup31 *g131  
   c322 = c322  - 0.5d0*gup31 *g221  
   c323 = c323  - 0.5d0*gup31 *g231  
   c333 = c333  - 0.5d0*gup31 *g331  
 
 
return 
end 
!----------------------------------------------------------------------- 
      subroutine connect8(g11, g12, g13, g22, g23, g33,      &
                          gup11, gup12, gup31, gup22, gup32, gup33,      &
                          c111, c112, c113, c122, c123, c133,      &
                          c211, c212, c213, c222, c223, c233,      &
                          c311, c312, c313, c322, c323, c333,      &
                          g112, g122, g132, g222, g232, g332,       &
                          g131, g333, g133, h) 
 
implicit none
real*8                        :: h
real*8,dimension(:,:,:),intent(IN)    :: g11,g12,g13,g22,g23,g33
real*8,dimension(:,:,:),intent(IN)    :: gup11,gup12,gup31,gup22,gup32,gup33

real*8,dimension(:,:,:),intent(OUT)    :: c111,c112,c113,c122,c123,c133
real*8,dimension(:,:,:),intent(OUT)    :: c211,c212,c213,c222,c223,c233
real*8,dimension(:,:,:),intent(OUT)    :: c311,c312,c313,c322,c323,c333

real*8,dimension(:,:,:),intent(IN)  :: g112,g122,g132,g222,     &
                                       g232,g332,g131,g333,g133 
 
   c111 = c111  - 0.5d0*gup12 *g112  
   c112 = c112  - 0.5d0*gup12 *g122  
   c113 = c113  - 0.5d0*gup12 *g132  
   c122 = c122  - 0.5d0*gup12 *g222  
   c123 = c123  - 0.5d0*gup12 *g232  
   c133 = c133  - 0.5d0*gup12 *g332  
 
   c211 = c211  - 0.5d0*gup22 *g112  
   c212 = c212  - 0.5d0*gup22 *g122  
   c213 = c213  - 0.5d0*gup22 *g132  
   c222 = c222  - 0.5d0*gup22 *g222  
   c223 = c223  - 0.5d0*gup22 *g232  
   c233 = c233  - 0.5d0*gup22 *g332  

   c311 = c311  - 0.5d0*gup32 *g112  
   c312 = c312  - 0.5d0*gup32 *g122  
   c313 = c313  - 0.5d0*gup32 *g132  
   c322 = c322  - 0.5d0*gup32 *g222  
   c323 = c323  - 0.5d0*gup32 *g232  
   c333 = c333  - 0.5d0*gup32 *g332  
 
 
return 
end 
!----------------------------------------------------------------------- 
subroutine connect9(g11, g12, g13, g22, g23, g33,      &
                    gup11, gup12, gup13, gup22, gup23, gup33,      &
                    c111, c112, c113, c122, c123, c133,      &
                    c211, c212, c213, c222, c223, c233,      &
                    c311, c312, c313, c322, c323, c333,      &
                    g113, g123, g133, g223, g233, g333,       &
                    g131, g132, g222, h) 
 
implicit none
real*8                        :: h
real*8,dimension(:,:,:),intent(IN)    :: g11,g12,g13,g22,g23,g33
real*8,dimension(:,:,:),intent(IN)    :: gup11,gup12,gup13,gup22,gup23,gup33

real*8,dimension(:,:,:),intent(OUT)   :: c111,c112,c113,c122,c123,c133
real*8,dimension(:,:,:),intent(OUT)    :: c211,c212,c213,c222,c223,c233
real*8,dimension(:,:,:),intent(OUT)    :: c311,c312,c313,c322,c323,c333
real*8,dimension(:,:,:),intent(IN)    :: g113,g123,g133,g223,   &
                                         g233,g333,g131,g132,g222 
 
   c111 = c111  - 0.5d0*gup13 *g113  
   c112 = c112  - 0.5d0*gup13 *g123  
   c113 = c113  - 0.5d0*gup13 *g133  
   c122 = c122  - 0.5d0*gup13 *g223  
   c123 = c123  - 0.5d0*gup13 *g233  
   c133 = c133  - 0.5d0*gup13 *g333  
 
   c211 = c211  - 0.5d0*gup23 *g113  
   c212 = c212  - 0.5d0*gup23 *g123  
   c213 = c213  - 0.5d0*gup23 *g133  
   c222 = c222  - 0.5d0*gup23 *g223  
   c223 = c223  - 0.5d0*gup23 *g233  
   c233 = c233  - 0.5d0*gup23 *g333  

   c311 = c311  - 0.5d0*gup33 *g113  
   c312 = c312  - 0.5d0*gup33 *g123  
   c313 = c313  - 0.5d0*gup33 *g133  
   c322 = c322  - 0.5d0*gup33 *g223  
   c323 = c323  - 0.5d0*gup33 *g233  
   c333 = c333  - 0.5d0*gup33 *g333  
 
 
return 
end 
