c-----------------------------------------------------------------------
c connect.f
*     Randall R. Correll
*     Mijan F. Huq
*     Scott A. Klasky
*     Richard A. Matzner
*
*     COPYRIGHT 1995-1996
c-----------------------------------------------------------------------
c This routine returns covariant form of ricci tensor. Raise index before
c using in evolution code.

      subroutine connect(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,
     *           t1, t2, t3, t4, t5, t6, t7, t8, t9,
     *           r11, r12, r13, r22, r23, r33, h, nx, ny, nz,nl,m,d)
      implicit none

      integer  nx, ny, nz
      real*8 h

      real*8   g11(nx,ny,nz), g12(nx,ny,nz), g13(nx,ny,nz), 
     *           g22(nx,ny,nz), g23(nx,ny,nz), g33(nx,ny,nz),
     *           gup11(nx,ny,nz), gup12(nx,ny,nz), gup13(nx,ny,nz), 
     *           gup22(nx,ny,nz), gup23(nx,ny,nz), gup33(nx,ny,nz)

      real*8   c111(nx,ny,nz), c112(nx,ny,nz), c113(nx,ny,nz), 
     *           c122(nx,ny,nz), c123(nx,ny,nz), c133(nx,ny,nz),
     *           c211(nx,ny,nz), c212(nx,ny,nz), c213(nx,ny,nz), 
     *           c222(nx,ny,nz), c223(nx,ny,nz), c233(nx,ny,nz),
     *           c311(nx,ny,nz), c312(nx,ny,nz), c313(nx,ny,nz), 
     *           c322(nx,ny,nz), c323(nx,ny,nz), c333(nx,ny,nz)

      real*8   t1(nx,ny,nz), t2(nx,ny,nz), t3(nx,ny,nz), t4(nx,ny,nz), 
     *           t5(nx,ny,nz), t6(nx,ny,nz), t7(nx,ny,nz), t8(nx,ny,nz), 
     *           t9(nx,ny,nz)

      real*8   r11(nx,ny,nz), r12(nx,ny,nz), r13(nx,ny,nz), 
     *           r22(nx,ny,nz), r23(nx,ny,nz), r33(nx,ny,nz)

		integer nl
		logical ltrace
		parameter ( ltrace = .false. )

		real*8  d(3*nl)
		integer m(3*nl)

c      call subsidiary routines to calculate  Christoffell tensor piecewise.

      call setval(nx,ny,nz,c111,0.0e0)
      call setval(nx,ny,nz,c112,0.0e0)
      call setval(nx,ny,nz,c113,0.0e0)
      call setval(nx,ny,nz,c122,0.0e0)
      call setval(nx,ny,nz,c123,0.0e0)
      call setval(nx,ny,nz,c133,0.0e0)

      call setval(nx,ny,nz,c211,0.0e0)
      call setval(nx,ny,nz,c212,0.0e0)
      call setval(nx,ny,nz,c213,0.0e0)
      call setval(nx,ny,nz,c222,0.0e0)
      call setval(nx,ny,nz,c223,0.0e0)
      call setval(nx,ny,nz,c233,0.0e0)

      call setval(nx,ny,nz,c311,0.0e0)
      call setval(nx,ny,nz,c312,0.0e0)
      call setval(nx,ny,nz,c313,0.0e0)
      call setval(nx,ny,nz,c322,0.0e0)
      call setval(nx,ny,nz,c323,0.0e0)
      call setval(nx,ny,nz,c333,0.0e0)

      call diff(nx,ny,nz,h,g11,t1,'x',m,d,nl)
      call diff(nx,ny,nz,h,g11,t2,'y',m,d,nl)
      call diff(nx,ny,nz,h,g11,t3,'z',m,d,nl)

      call diff(nx,ny,nz,h,g12,t4,'x',m,d,nl)
      call diff(nx,ny,nz,h,g12,t5,'y',m,d,nl)
      call diff(nx,ny,nz,h,g12,t6,'z',m,d,nl)

      call diff(nx,ny,nz,h,g13,t7,'x',m,d,nl)
      call diff(nx,ny,nz,h,g13,t8,'y',m,d,nl)
      call diff(nx,ny,nz,h,g13,t9,'z',m,d,nl)

C note that the first 6 t's are really used, the last 3 are not....

		if (ltrace) write(0,*) 'calling connect1'
      call connect1(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,
     *                t1,   t2  ,   t3,   t4,   t5,   t6, 
     *                t7,   t8  ,   t9, h, nx , ny , nz)

		if (ltrace) write(0,*) 'calling connect4'
      call connect4(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,
     *                t1,   t4  ,   t7,   t5,   t8,   t9,
     *                t2,   t3  ,   t6, h, nx , ny , nz)

      call diff(nx,ny,nz,h,g12,t1,'x',m,d,nl)
      call diff(nx,ny,nz,h,g12,t2,'y',m,d,nl)
      call diff(nx,ny,nz,h,g12,t3,'z',m,d,nl)

      call diff(nx,ny,nz,h,g22,t4,'x',m,d,nl)
      call diff(nx,ny,nz,h,g22,t5,'y',m,d,nl)
      call diff(nx,ny,nz,h,g22,t6,'z',m,d,nl)

      call diff(nx,ny,nz,h,g23,t7,'x',m,d,nl)
      call diff(nx,ny,nz,h,g23,t8,'y',m,d,nl)
      call diff(nx,ny,nz,h,g23,t9,'z',m,d,nl)

		if (ltrace) write(0,*) 'calling connect2'
      call connect2(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,
     *                t1,   t2  ,   t3,   t5,   t6,   t9,
     *                t7,   t8  ,   t6, h, nx , ny , nz)


		if (ltrace) write(0,*) 'calling connect5'
      call connect5(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,
     *                t1,   t4  ,   t7,   t5,   t8,   t9,
     *                t2,   t3  ,   t6, h, nx , ny , nz)

      call diff(nx,ny,nz,h,g13,t1,'x',m,d,nl)
      call diff(nx,ny,nz,h,g13,t2,'y',m,d,nl)
      call diff(nx,ny,nz,h,g13,t3,'z',m,d,nl)

      call diff(nx,ny,nz,h,g23,t4,'x',m,d,nl)
      call diff(nx,ny,nz,h,g23,t5,'y',m,d,nl)
      call diff(nx,ny,nz,h,g23,t6,'z',m,d,nl)

      call diff(nx,ny,nz,h,g33,t7,'x',m,d,nl)
      call diff(nx,ny,nz,h,g33,t8,'y',m,d,nl)
      call diff(nx,ny,nz,h,g33,t9,'z',m,d,nl)

		if (ltrace) write(0,*) 'calling connect3'
      call 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,
     *                t1,   t2  ,   t3,   t5,   t6,   t9,
     *                t7,   t8  ,   t4, h, nx , ny , nz)

		if (ltrace) write(0,*) 'calling connect6'
      call 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,
     *                t1,   t4  ,   t7,   t5,   t8,   t9,
     *                t2,   t3  ,   t6, h, nx , ny , nz)


C now we calculate the last part ...

      call diff(nx,ny,nz,h,g11,t1,'x',m,d,nl)
      call diff(nx,ny,nz,h,g12,t2,'x',m,d,nl)
      call diff(nx,ny,nz,h,g13,t3,'x',m,d,nl)

      call diff(nx,ny,nz,h,g22,t4,'x',m,d,nl)
      call diff(nx,ny,nz,h,g23,t5,'x',m,d,nl)
      call diff(nx,ny,nz,h,g33,t6,'x',m,d,nl)

		if (ltrace) write(0,*) 'calling connect7'
      call connect7(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,
     *                t1,   t2  ,   t3,   t4,   t5,   t6,
     *                t7,   t8  ,   t9, h, nx , ny , nz)

      call diff(nx,ny,nz,h,g11,t1,'y',m,d,nl)
      call diff(nx,ny,nz,h,g12,t2,'y',m,d,nl)
      call diff(nx,ny,nz,h,g13,t3,'y',m,d,nl)

      call diff(nx,ny,nz,h,g22,t4,'y',m,d,nl)
      call diff(nx,ny,nz,h,g23,t5,'y',m,d,nl)
      call diff(nx,ny,nz,h,g33,t6,'y',m,d,nl)

		if (ltrace) write(0,*) 'calling connect8'
      call connect8(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,
     *                t1,   t2  ,   t3,   t4,   t5,   t6,
     *                t7,   t8  ,   t9, h, nx , ny , nz)

      call diff(nx,ny,nz,h,g11,t1,'z',m,d,nl)
      call diff(nx,ny,nz,h,g12,t2,'z',m,d,nl)
      call diff(nx,ny,nz,h,g13,t3,'z',m,d,nl)

      call diff(nx,ny,nz,h,g22,t4,'z',m,d,nl)
      call diff(nx,ny,nz,h,g23,t5,'z',m,d,nl)
      call diff(nx,ny,nz,h,g33,t6,'z',m,d,nl)

		if (ltrace) write(0,*) 'calling connect9'
      call 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,
     *                t1,   t2  ,   t3,   t4,   t5,   t6,
     *                t7,   t8  ,   t9, h, nx , ny , nz)

      return 
      end
c-----------------------------------------------------------------------
C 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, nx , ny , nz)

      implicit none

      integer  nx, ny, nz
      real*8   g11(nx,ny,nz), g12(nx,ny,nz), g13(nx,ny,nz),
     *           g22(nx,ny,nz), g23(nx,ny,nz), g33(nx,ny,nz),
     *           gup11(nx,ny,nz), gup21(nx,ny,nz), gup31(nx,ny,nz),
     *           gup22(nx,ny,nz), gup23(nx,ny,nz), gup33(nx,ny,nz)

      real*8   c111(nx,ny,nz), c112(nx,ny,nz), c113(nx,ny,nz),
     *           c122(nx,ny,nz), c123(nx,ny,nz), c133(nx,ny,nz),
     *           c211(nx,ny,nz), c212(nx,ny,nz), c213(nx,ny,nz),
     *           c222(nx,ny,nz), c223(nx,ny,nz), c233(nx,ny,nz),
     *           c311(nx,ny,nz), c312(nx,ny,nz), c313(nx,ny,nz),
     *           c322(nx,ny,nz), c323(nx,ny,nz), c333(nx,ny,nz)

      real*8 g111(nx,ny,nz),g112(nx,ny,nz),g113(nx,ny,nz)
      real*8 g121(nx,ny,nz),g122(nx,ny,nz),g123(nx,ny,nz)
      real*8 g131(nx,ny,nz),g132(nx,ny,nz),g133(nx,ny,nz)

      real*8 h

      integer i,j,k

      do k = 1 , nz
         do j = 1 ,ny
            do i = 1 , nx

               c111(i,j,k)= c111(i,j,k) + 0.5e0*gup11(i,j,k)*g111(i,j,k)
               c112(i,j,k)= c112(i,j,k) + 0.5e0*gup11(i,j,k)*g112(i,j,k)
               c113(i,j,k)= c113(i,j,k) + 0.5e0*gup11(i,j,k)*g113(i,j,k)
               c122(i,j,k)= c122(i,j,k) + 0.5e0*gup11(i,j,k)*g122(i,j,k)
               c123(i,j,k)= c123(i,j,k) + 0.5e0*gup11(i,j,k)*g123(i,j,k)
               c133(i,j,k)= c133(i,j,k) + 0.5e0*gup11(i,j,k)*g133(i,j,k)

               c211(i,j,k)= c211(i,j,k) + 0.5e0*gup21(i,j,k)*g111(i,j,k)
               c212(i,j,k)= c212(i,j,k) + 0.5e0*gup21(i,j,k)*g112(i,j,k)
               c213(i,j,k)= c213(i,j,k) + 0.5e0*gup21(i,j,k)*g113(i,j,k)
               c222(i,j,k)= c222(i,j,k) + 0.5e0*gup21(i,j,k)*g122(i,j,k)
               c223(i,j,k)= c223(i,j,k) + 0.5e0*gup21(i,j,k)*g123(i,j,k)
               c233(i,j,k)= c233(i,j,k) + 0.5e0*gup21(i,j,k)*g133(i,j,k)

               c311(i,j,k)= c311(i,j,k) + 0.5e0*gup31(i,j,k)*g111(i,j,k)
               c312(i,j,k)= c312(i,j,k) + 0.5e0*gup31(i,j,k)*g112(i,j,k)
               c313(i,j,k)= c313(i,j,k) + 0.5e0*gup31(i,j,k)*g113(i,j,k)
               c322(i,j,k)= c322(i,j,k) + 0.5e0*gup31(i,j,k)*g122(i,j,k)
               c323(i,j,k)= c323(i,j,k) + 0.5e0*gup31(i,j,k)*g123(i,j,k)
               c333(i,j,k)= c333(i,j,k) + 0.5e0*gup31(i,j,k)*g133(i,j,k)

            end do
         end do
      end do

      return
      end
c-----------------------------------------------------------------------
C 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, nx , ny , nz)

      implicit none

      integer  nx, ny, nz
      real*8   g11(nx,ny,nz), g12(nx,ny,nz), g13(nx,ny,nz),
     *           g22(nx,ny,nz), g23(nx,ny,nz), g33(nx,ny,nz),
     *           gup11(nx,ny,nz), gup12(nx,ny,nz), gup31(nx,ny,nz),
     *           gup22(nx,ny,nz), gup32(nx,ny,nz), gup33(nx,ny,nz)

      real*8   c111(nx,ny,nz), c112(nx,ny,nz), c113(nx,ny,nz),
     *           c122(nx,ny,nz), c123(nx,ny,nz), c133(nx,ny,nz),
     *           c211(nx,ny,nz), c212(nx,ny,nz), c213(nx,ny,nz),
     *           c222(nx,ny,nz), c223(nx,ny,nz), c233(nx,ny,nz),
     *           c311(nx,ny,nz), c312(nx,ny,nz), c313(nx,ny,nz),
     *           c322(nx,ny,nz), c323(nx,ny,nz), c333(nx,ny,nz)

      real*8 g211(nx,ny,nz),g212(nx,ny,nz),g213(nx,ny,nz)
      real*8 g222(nx,ny,nz),g223(nx,ny,nz),g233(nx,ny,nz)
      real*8 g111(nx,ny,nz),g112(nx,ny,nz),g333(nx,ny,nz)

      real*8 h

      integer i,j,k

      do k = 1 , nz
         do j = 1 ,ny
            do i = 1 , nx

               c111(i,j,k)= c111(i,j,k) + 0.5e0*gup12(i,j,k)*g211(i,j,k)
               c112(i,j,k)= c112(i,j,k) + 0.5e0*gup12(i,j,k)*g212(i,j,k)
               c113(i,j,k)= c113(i,j,k) + 0.5e0*gup12(i,j,k)*g213(i,j,k)
               c122(i,j,k)= c122(i,j,k) + 0.5e0*gup12(i,j,k)*g222(i,j,k)
               c123(i,j,k)= c123(i,j,k) + 0.5e0*gup12(i,j,k)*g223(i,j,k)
               c133(i,j,k)= c133(i,j,k) + 0.5e0*gup12(i,j,k)*g233(i,j,k)

               c211(i,j,k)= c211(i,j,k) + 0.5e0*gup22(i,j,k)*g211(i,j,k)
               c212(i,j,k)= c212(i,j,k) + 0.5e0*gup22(i,j,k)*g212(i,j,k)
               c213(i,j,k)= c213(i,j,k) + 0.5e0*gup22(i,j,k)*g213(i,j,k)
               c222(i,j,k)= c222(i,j,k) + 0.5e0*gup22(i,j,k)*g222(i,j,k)
               c223(i,j,k)= c223(i,j,k) + 0.5e0*gup22(i,j,k)*g223(i,j,k)
               c233(i,j,k)= c233(i,j,k) + 0.5e0*gup22(i,j,k)*g233(i,j,k)

               c311(i,j,k)= c311(i,j,k) + 0.5e0*gup32(i,j,k)*g211(i,j,k)
               c312(i,j,k)= c312(i,j,k) + 0.5e0*gup32(i,j,k)*g212(i,j,k)
               c313(i,j,k)= c313(i,j,k) + 0.5e0*gup32(i,j,k)*g213(i,j,k)
               c322(i,j,k)= c322(i,j,k) + 0.5e0*gup32(i,j,k)*g222(i,j,k)
               c323(i,j,k)= c323(i,j,k) + 0.5e0*gup32(i,j,k)*g223(i,j,k)
               c333(i,j,k)= c333(i,j,k) + 0.5e0*gup32(i,j,k)*g233(i,j,k)

            end do
         end do
      end do

      return
      end

c-----------------------------------------------------------------------
      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, nx , ny , nz)

      implicit none

      integer  nx, ny, nz
      real*8   g11(nx,ny,nz), g12(nx,ny,nz), g13(nx,ny,nz),
     *           g22(nx,ny,nz), g23(nx,ny,nz), g33(nx,ny,nz),
     *           gup11(nx,ny,nz), gup12(nx,ny,nz), gup13(nx,ny,nz),
     *           gup22(nx,ny,nz), gup23(nx,ny,nz), gup33(nx,ny,nz)

      real*8   c111(nx,ny,nz), c112(nx,ny,nz), c113(nx,ny,nz),
     *           c122(nx,ny,nz), c123(nx,ny,nz), c133(nx,ny,nz),
     *           c211(nx,ny,nz), c212(nx,ny,nz), c213(nx,ny,nz),
     *           c222(nx,ny,nz), c223(nx,ny,nz), c233(nx,ny,nz),
     *           c311(nx,ny,nz), c312(nx,ny,nz), c313(nx,ny,nz),
     *           c322(nx,ny,nz), c323(nx,ny,nz), c333(nx,ny,nz)

      real*8 g311(nx,ny,nz),g312(nx,ny,nz),g313(nx,ny,nz)
      real*8 g322(nx,ny,nz),g323(nx,ny,nz),g333(nx,ny,nz)
      real*8 g332(nx,ny,nz),g112(nx,ny,nz),g233(nx,ny,nz)

      real*8 h

      integer i,j,k

      do k = 1 , nz
         do j = 1 ,ny
            do i = 1 , nx

               c111(i,j,k)= c111(i,j,k) + 0.5e0*gup13(i,j,k)*g311(i,j,k)
               c112(i,j,k)= c112(i,j,k) + 0.5e0*gup13(i,j,k)*g312(i,j,k)
               c113(i,j,k)= c113(i,j,k) + 0.5e0*gup13(i,j,k)*g313(i,j,k)
               c122(i,j,k)= c122(i,j,k) + 0.5e0*gup13(i,j,k)*g322(i,j,k)
               c123(i,j,k)= c123(i,j,k) + 0.5e0*gup13(i,j,k)*g323(i,j,k)
               c133(i,j,k)= c133(i,j,k) + 0.5e0*gup13(i,j,k)*g333(i,j,k)

               c211(i,j,k)= c211(i,j,k) + 0.5e0*gup23(i,j,k)*g311(i,j,k)
               c212(i,j,k)= c212(i,j,k) + 0.5e0*gup23(i,j,k)*g312(i,j,k)
               c213(i,j,k)= c213(i,j,k) + 0.5e0*gup23(i,j,k)*g313(i,j,k)
               c222(i,j,k)= c222(i,j,k) + 0.5e0*gup23(i,j,k)*g322(i,j,k)
               c223(i,j,k)= c223(i,j,k) + 0.5e0*gup23(i,j,k)*g323(i,j,k)
               c233(i,j,k)= c233(i,j,k) + 0.5e0*gup23(i,j,k)*g333(i,j,k)

               c311(i,j,k)= c311(i,j,k) + 0.5e0*gup33(i,j,k)*g311(i,j,k)
               c312(i,j,k)= c312(i,j,k) + 0.5e0*gup33(i,j,k)*g312(i,j,k)
               c313(i,j,k)= c313(i,j,k) + 0.5e0*gup33(i,j,k)*g313(i,j,k)
               c322(i,j,k)= c322(i,j,k) + 0.5e0*gup33(i,j,k)*g322(i,j,k)
               c323(i,j,k)= c323(i,j,k) + 0.5e0*gup33(i,j,k)*g323(i,j,k)
               c333(i,j,k)= c333(i,j,k) + 0.5e0*gup33(i,j,k)*g333(i,j,k)

            end do
         end do
      end do

      return
      end

c-----------------------------------------------------------------------
      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, nx , ny , nz)

      implicit none

      integer  nx, ny, nz
      real*8   g11(nx,ny,nz), g12(nx,ny,nz), g13(nx,ny,nz),
     *           g22(nx,ny,nz), g23(nx,ny,nz), g33(nx,ny,nz),
     *           gup11(nx,ny,nz), gup21(nx,ny,nz), gup31(nx,ny,nz),
     *           gup22(nx,ny,nz), gup32(nx,ny,nz), gup33(nx,ny,nz)

      real*8   c111(nx,ny,nz), c112(nx,ny,nz), c113(nx,ny,nz),
     *           c122(nx,ny,nz), c123(nx,ny,nz), c133(nx,ny,nz),
     *           c211(nx,ny,nz), c212(nx,ny,nz), c213(nx,ny,nz),
     *           c222(nx,ny,nz), c223(nx,ny,nz), c233(nx,ny,nz),
     *           c311(nx,ny,nz), c312(nx,ny,nz), c313(nx,ny,nz),
     *           c322(nx,ny,nz), c323(nx,ny,nz), c333(nx,ny,nz)

      real*8 g111(nx,ny,nz),g121(nx,ny,nz),g131(nx,ny,nz)
      real*8 g122(nx,ny,nz),g132(nx,ny,nz),g133(nx,ny,nz)
      real*8 g222(nx,ny,nz),g333(nx,ny,nz),g311(nx,ny,nz)

      real*8 h

      integer i,j,k

      do k = 1 , nz
         do j = 1 ,ny
            do i = 1 , nx

               c111(i,j,k)= c111(i,j,k) + 0.5e0*gup11(i,j,k)*g111(i,j,k)
               c112(i,j,k)= c112(i,j,k) + 0.5e0*gup11(i,j,k)*g121(i,j,k)
               c113(i,j,k)= c113(i,j,k) + 0.5e0*gup11(i,j,k)*g131(i,j,k)
               c122(i,j,k)= c122(i,j,k) + 0.5e0*gup11(i,j,k)*g122(i,j,k)
               c123(i,j,k)= c123(i,j,k) + 0.5e0*gup11(i,j,k)*g132(i,j,k)
               c133(i,j,k)= c133(i,j,k) + 0.5e0*gup11(i,j,k)*g133(i,j,k)

               c211(i,j,k)= c211(i,j,k) + 0.5e0*gup21(i,j,k)*g111(i,j,k)
               c212(i,j,k)= c212(i,j,k) + 0.5e0*gup21(i,j,k)*g121(i,j,k)
               c213(i,j,k)= c213(i,j,k) + 0.5e0*gup21(i,j,k)*g131(i,j,k)
               c222(i,j,k)= c222(i,j,k) + 0.5e0*gup21(i,j,k)*g122(i,j,k)
               c223(i,j,k)= c223(i,j,k) + 0.5e0*gup21(i,j,k)*g132(i,j,k)
               c233(i,j,k)= c233(i,j,k) + 0.5e0*gup21(i,j,k)*g133(i,j,k)

               c311(i,j,k)= c311(i,j,k) + 0.5e0*gup31(i,j,k)*g111(i,j,k)
               c312(i,j,k)= c312(i,j,k) + 0.5e0*gup31(i,j,k)*g121(i,j,k)
               c313(i,j,k)= c313(i,j,k) + 0.5e0*gup31(i,j,k)*g131(i,j,k)
               c322(i,j,k)= c322(i,j,k) + 0.5e0*gup31(i,j,k)*g122(i,j,k)
               c323(i,j,k)= c323(i,j,k) + 0.5e0*gup31(i,j,k)*g132(i,j,k)
               c333(i,j,k)= c333(i,j,k) + 0.5e0*gup31(i,j,k)*g133(i,j,k)

            end do
         end do
      end do

      return
      end
c-----------------------------------------------------------------------
      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, nx , ny , nz)

      implicit none

      integer  nx, ny, nz
      real*8   g11(nx,ny,nz), g12(nx,ny,nz), g13(nx,ny,nz),
     *           g22(nx,ny,nz), g23(nx,ny,nz), g33(nx,ny,nz),
     *           gup11(nx,ny,nz), gup12(nx,ny,nz), gup31(nx,ny,nz),
     *           gup22(nx,ny,nz), gup32(nx,ny,nz), gup33(nx,ny,nz)

      real*8   c111(nx,ny,nz), c112(nx,ny,nz), c113(nx,ny,nz),
     *           c122(nx,ny,nz), c123(nx,ny,nz), c133(nx,ny,nz),
     *           c211(nx,ny,nz), c212(nx,ny,nz), c213(nx,ny,nz),
     *           c222(nx,ny,nz), c223(nx,ny,nz), c233(nx,ny,nz),
     *           c311(nx,ny,nz), c312(nx,ny,nz), c313(nx,ny,nz),
     *           c322(nx,ny,nz), c323(nx,ny,nz), c333(nx,ny,nz)

      real*8 g211(nx,ny,nz),g221(nx,ny,nz),g231(nx,ny,nz)
      real*8 g222(nx,ny,nz),g232(nx,ny,nz),g233(nx,ny,nz)
      real*8 g131(nx,ny,nz),g132(nx,ny,nz),g133(nx,ny,nz)

      real*8 h

      integer i,j,k

      do k = 1 , nz
         do j = 1 ,ny
            do i = 1 , nx

               c111(i,j,k)= c111(i,j,k) + 0.5e0*gup12(i,j,k)*g211(i,j,k)
               c112(i,j,k)= c112(i,j,k) + 0.5e0*gup12(i,j,k)*g221(i,j,k)
               c113(i,j,k)= c113(i,j,k) + 0.5e0*gup12(i,j,k)*g231(i,j,k)
               c122(i,j,k)= c122(i,j,k) + 0.5e0*gup12(i,j,k)*g222(i,j,k)
               c123(i,j,k)= c123(i,j,k) + 0.5e0*gup12(i,j,k)*g232(i,j,k)
               c133(i,j,k)= c133(i,j,k) + 0.5e0*gup12(i,j,k)*g233(i,j,k)

               c211(i,j,k)= c211(i,j,k) + 0.5e0*gup22(i,j,k)*g211(i,j,k)
               c212(i,j,k)= c212(i,j,k) + 0.5e0*gup22(i,j,k)*g221(i,j,k)
               c213(i,j,k)= c213(i,j,k) + 0.5e0*gup22(i,j,k)*g231(i,j,k)
               c222(i,j,k)= c222(i,j,k) + 0.5e0*gup22(i,j,k)*g222(i,j,k)
               c223(i,j,k)= c223(i,j,k) + 0.5e0*gup22(i,j,k)*g232(i,j,k)
               c233(i,j,k)= c233(i,j,k) + 0.5e0*gup22(i,j,k)*g233(i,j,k)

               c311(i,j,k)= c311(i,j,k) + 0.5e0*gup32(i,j,k)*g211(i,j,k)
               c312(i,j,k)= c312(i,j,k) + 0.5e0*gup32(i,j,k)*g221(i,j,k)
               c313(i,j,k)= c313(i,j,k) + 0.5e0*gup32(i,j,k)*g231(i,j,k)
               c322(i,j,k)= c322(i,j,k) + 0.5e0*gup32(i,j,k)*g222(i,j,k)
               c323(i,j,k)= c323(i,j,k) + 0.5e0*gup32(i,j,k)*g232(i,j,k)
               c333(i,j,k)= c333(i,j,k) + 0.5e0*gup32(i,j,k)*g233(i,j,k)

            end do
         end do
      end do

      return
      end
c-----------------------------------------------------------------------
      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, nx , ny , nz)

      implicit none

      integer  nx, ny, nz
      real*8   g11(nx,ny,nz), g12(nx,ny,nz), g13(nx,ny,nz),
     *           g22(nx,ny,nz), g23(nx,ny,nz), g33(nx,ny,nz),
     *           gup11(nx,ny,nz), gup12(nx,ny,nz), gup13(nx,ny,nz),
     *           gup22(nx,ny,nz), gup23(nx,ny,nz), gup33(nx,ny,nz)

      real*8   c111(nx,ny,nz), c112(nx,ny,nz), c113(nx,ny,nz),
     *           c122(nx,ny,nz), c123(nx,ny,nz), c133(nx,ny,nz),
     *           c211(nx,ny,nz), c212(nx,ny,nz), c213(nx,ny,nz),
     *           c222(nx,ny,nz), c223(nx,ny,nz), c233(nx,ny,nz),
     *           c311(nx,ny,nz), c312(nx,ny,nz), c313(nx,ny,nz),
     *           c322(nx,ny,nz), c323(nx,ny,nz), c333(nx,ny,nz)

      real*8 g311(nx,ny,nz),g321(nx,ny,nz),g331(nx,ny,nz)
      real*8 g322(nx,ny,nz),g332(nx,ny,nz),g333(nx,ny,nz)
      real*8 g131(nx,ny,nz),g132(nx,ny,nz),g133(nx,ny,nz)

      real*8 h

      integer i,j,k

      do k = 1 , nz
         do j = 1 ,ny
            do i = 1 , nx

               c111(i,j,k)= c111(i,j,k) + 0.5e0*gup13(i,j,k)*g311(i,j,k)
               c112(i,j,k)= c112(i,j,k) + 0.5e0*gup13(i,j,k)*g321(i,j,k)
               c113(i,j,k)= c113(i,j,k) + 0.5e0*gup13(i,j,k)*g331(i,j,k)
               c122(i,j,k)= c122(i,j,k) + 0.5e0*gup13(i,j,k)*g322(i,j,k)
               c123(i,j,k)= c123(i,j,k) + 0.5e0*gup13(i,j,k)*g332(i,j,k)
               c133(i,j,k)= c133(i,j,k) + 0.5e0*gup13(i,j,k)*g333(i,j,k)

               c211(i,j,k)= c211(i,j,k) + 0.5e0*gup23(i,j,k)*g311(i,j,k)
               c212(i,j,k)= c212(i,j,k) + 0.5e0*gup23(i,j,k)*g321(i,j,k)
               c213(i,j,k)= c213(i,j,k) + 0.5e0*gup23(i,j,k)*g331(i,j,k)
               c222(i,j,k)= c222(i,j,k) + 0.5e0*gup23(i,j,k)*g322(i,j,k)
               c223(i,j,k)= c223(i,j,k) + 0.5e0*gup23(i,j,k)*g332(i,j,k)
               c233(i,j,k)= c233(i,j,k) + 0.5e0*gup23(i,j,k)*g333(i,j,k)

               c311(i,j,k)= c311(i,j,k) + 0.5e0*gup33(i,j,k)*g311(i,j,k)
               c312(i,j,k)= c312(i,j,k) + 0.5e0*gup33(i,j,k)*g321(i,j,k)
               c313(i,j,k)= c313(i,j,k) + 0.5e0*gup33(i,j,k)*g331(i,j,k)
               c322(i,j,k)= c322(i,j,k) + 0.5e0*gup33(i,j,k)*g322(i,j,k)
               c323(i,j,k)= c323(i,j,k) + 0.5e0*gup33(i,j,k)*g332(i,j,k)
               c333(i,j,k)= c333(i,j,k) + 0.5e0*gup33(i,j,k)*g333(i,j,k)

            end do
         end do
      end do

      return
      end
c-----------------------------------------------------------------------
      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, nx , ny , nz)

      implicit none

      integer  nx, ny, nz
      real*8   g11(nx,ny,nz), g12(nx,ny,nz), g13(nx,ny,nz),
     *           g22(nx,ny,nz), g23(nx,ny,nz), g33(nx,ny,nz),
     *           gup11(nx,ny,nz), gup21(nx,ny,nz), gup31(nx,ny,nz),
     *           gup22(nx,ny,nz), gup32(nx,ny,nz), gup33(nx,ny,nz)

      real*8   c111(nx,ny,nz), c112(nx,ny,nz), c113(nx,ny,nz),
     *           c122(nx,ny,nz), c123(nx,ny,nz), c133(nx,ny,nz),
     *           c211(nx,ny,nz), c212(nx,ny,nz), c213(nx,ny,nz),
     *           c222(nx,ny,nz), c223(nx,ny,nz), c233(nx,ny,nz),
     *           c311(nx,ny,nz), c312(nx,ny,nz), c313(nx,ny,nz),
     *           c322(nx,ny,nz), c323(nx,ny,nz), c333(nx,ny,nz)

      real*8 g111(nx,ny,nz),g121(nx,ny,nz),g131(nx,ny,nz)
      real*8 g221(nx,ny,nz),g231(nx,ny,nz),g331(nx,ny,nz)
      real*8 g333(nx,ny,nz),g132(nx,ny,nz),g133(nx,ny,nz)

      real*8 h

      integer i,j,k

      do k = 1 , nz
         do j = 1 ,ny
            do i = 1 , nx

               c111(i,j,k)= c111(i,j,k) - 0.5e0*gup11(i,j,k)*g111(i,j,k)
               c112(i,j,k)= c112(i,j,k) - 0.5e0*gup11(i,j,k)*g121(i,j,k)
               c113(i,j,k)= c113(i,j,k) - 0.5e0*gup11(i,j,k)*g131(i,j,k)
               c122(i,j,k)= c122(i,j,k) - 0.5e0*gup11(i,j,k)*g221(i,j,k)
               c123(i,j,k)= c123(i,j,k) - 0.5e0*gup11(i,j,k)*g231(i,j,k)
               c133(i,j,k)= c133(i,j,k) - 0.5e0*gup11(i,j,k)*g331(i,j,k)

               c211(i,j,k)= c211(i,j,k) - 0.5e0*gup21(i,j,k)*g111(i,j,k)
               c212(i,j,k)= c212(i,j,k) - 0.5e0*gup21(i,j,k)*g121(i,j,k)
               c213(i,j,k)= c213(i,j,k) - 0.5e0*gup21(i,j,k)*g131(i,j,k)
               c222(i,j,k)= c222(i,j,k) - 0.5e0*gup21(i,j,k)*g221(i,j,k)
               c223(i,j,k)= c223(i,j,k) - 0.5e0*gup21(i,j,k)*g231(i,j,k)
               c233(i,j,k)= c233(i,j,k) - 0.5e0*gup21(i,j,k)*g331(i,j,k)

               c311(i,j,k)= c311(i,j,k) - 0.5e0*gup31(i,j,k)*g111(i,j,k)
               c312(i,j,k)= c312(i,j,k) - 0.5e0*gup31(i,j,k)*g121(i,j,k)
               c313(i,j,k)= c313(i,j,k) - 0.5e0*gup31(i,j,k)*g131(i,j,k)
               c322(i,j,k)= c322(i,j,k) - 0.5e0*gup31(i,j,k)*g221(i,j,k)
               c323(i,j,k)= c323(i,j,k) - 0.5e0*gup31(i,j,k)*g231(i,j,k)
               c333(i,j,k)= c333(i,j,k) - 0.5e0*gup31(i,j,k)*g331(i,j,k)

            end do
         end do
      end do

      return
      end
c-----------------------------------------------------------------------
      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, nx , ny , nz)

      implicit none

      integer  nx, ny, nz
      real*8   g11(nx,ny,nz), g12(nx,ny,nz), g13(nx,ny,nz),
     *           g22(nx,ny,nz), g23(nx,ny,nz), g33(nx,ny,nz),
     *           gup11(nx,ny,nz), gup12(nx,ny,nz), gup31(nx,ny,nz),
     *           gup22(nx,ny,nz), gup32(nx,ny,nz), gup33(nx,ny,nz)

      real*8   c111(nx,ny,nz), c112(nx,ny,nz), c113(nx,ny,nz),
     *           c122(nx,ny,nz), c123(nx,ny,nz), c133(nx,ny,nz),
     *           c211(nx,ny,nz), c212(nx,ny,nz), c213(nx,ny,nz),
     *           c222(nx,ny,nz), c223(nx,ny,nz), c233(nx,ny,nz),
     *           c311(nx,ny,nz), c312(nx,ny,nz), c313(nx,ny,nz),
     *           c322(nx,ny,nz), c323(nx,ny,nz), c333(nx,ny,nz)

      real*8 g112(nx,ny,nz),g122(nx,ny,nz),g132(nx,ny,nz)
      real*8 g222(nx,ny,nz),g232(nx,ny,nz),g332(nx,ny,nz)
      real*8 g131(nx,ny,nz),g333(nx,ny,nz),g133(nx,ny,nz)

      real*8 h

      integer i,j,k

      do k = 1 , nz
         do j = 1 ,ny
            do i = 1 , nx

               c111(i,j,k)= c111(i,j,k) - 0.5e0*gup12(i,j,k)*g112(i,j,k)
               c112(i,j,k)= c112(i,j,k) - 0.5e0*gup12(i,j,k)*g122(i,j,k)
               c113(i,j,k)= c113(i,j,k) - 0.5e0*gup12(i,j,k)*g132(i,j,k)
               c122(i,j,k)= c122(i,j,k) - 0.5e0*gup12(i,j,k)*g222(i,j,k)
               c123(i,j,k)= c123(i,j,k) - 0.5e0*gup12(i,j,k)*g232(i,j,k)
               c133(i,j,k)= c133(i,j,k) - 0.5e0*gup12(i,j,k)*g332(i,j,k)

               c211(i,j,k)= c211(i,j,k) - 0.5e0*gup22(i,j,k)*g112(i,j,k)
               c212(i,j,k)= c212(i,j,k) - 0.5e0*gup22(i,j,k)*g122(i,j,k)
               c213(i,j,k)= c213(i,j,k) - 0.5e0*gup22(i,j,k)*g132(i,j,k)
               c222(i,j,k)= c222(i,j,k) - 0.5e0*gup22(i,j,k)*g222(i,j,k)
               c223(i,j,k)= c223(i,j,k) - 0.5e0*gup22(i,j,k)*g232(i,j,k)
               c233(i,j,k)= c233(i,j,k) - 0.5e0*gup22(i,j,k)*g332(i,j,k)

               c311(i,j,k)= c311(i,j,k) - 0.5e0*gup32(i,j,k)*g112(i,j,k)
               c312(i,j,k)= c312(i,j,k) - 0.5e0*gup32(i,j,k)*g122(i,j,k)
               c313(i,j,k)= c313(i,j,k) - 0.5e0*gup32(i,j,k)*g132(i,j,k)
               c322(i,j,k)= c322(i,j,k) - 0.5e0*gup32(i,j,k)*g222(i,j,k)
               c323(i,j,k)= c323(i,j,k) - 0.5e0*gup32(i,j,k)*g232(i,j,k)
               c333(i,j,k)= c333(i,j,k) - 0.5e0*gup32(i,j,k)*g332(i,j,k)

            end do
         end do
      end do

      return
      end
c-----------------------------------------------------------------------
      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, nx , ny , nz)

      implicit none

      integer  nx, ny, nz
      real*8   g11(nx,ny,nz), g12(nx,ny,nz), g13(nx,ny,nz),
     *           g22(nx,ny,nz), g23(nx,ny,nz), g33(nx,ny,nz),
     *           gup11(nx,ny,nz), gup12(nx,ny,nz), gup13(nx,ny,nz),
     *           gup22(nx,ny,nz), gup23(nx,ny,nz), gup33(nx,ny,nz)

      real*8   c111(nx,ny,nz), c112(nx,ny,nz), c113(nx,ny,nz),
     *           c122(nx,ny,nz), c123(nx,ny,nz), c133(nx,ny,nz),
     *           c211(nx,ny,nz), c212(nx,ny,nz), c213(nx,ny,nz),
     *           c222(nx,ny,nz), c223(nx,ny,nz), c233(nx,ny,nz),
     *           c311(nx,ny,nz), c312(nx,ny,nz), c313(nx,ny,nz),
     *           c322(nx,ny,nz), c323(nx,ny,nz), c333(nx,ny,nz)

      real*8 g113(nx,ny,nz),g123(nx,ny,nz),g133(nx,ny,nz)
      real*8 g223(nx,ny,nz),g233(nx,ny,nz),g333(nx,ny,nz)
      real*8 g131(nx,ny,nz),g132(nx,ny,nz),g222(nx,ny,nz)

      real*8 h

      integer i,j,k

      do k = 1 , nz
         do j = 1 ,ny
            do i = 1 , nx

               c111(i,j,k)= c111(i,j,k) - 0.5e0*gup13(i,j,k)*g113(i,j,k)
               c112(i,j,k)= c112(i,j,k) - 0.5e0*gup13(i,j,k)*g123(i,j,k)
               c113(i,j,k)= c113(i,j,k) - 0.5e0*gup13(i,j,k)*g133(i,j,k)
               c122(i,j,k)= c122(i,j,k) - 0.5e0*gup13(i,j,k)*g223(i,j,k)
               c123(i,j,k)= c123(i,j,k) - 0.5e0*gup13(i,j,k)*g233(i,j,k)
               c133(i,j,k)= c133(i,j,k) - 0.5e0*gup13(i,j,k)*g333(i,j,k)

               c211(i,j,k)= c211(i,j,k) - 0.5e0*gup23(i,j,k)*g113(i,j,k)
               c212(i,j,k)= c212(i,j,k) - 0.5e0*gup23(i,j,k)*g123(i,j,k)
               c213(i,j,k)= c213(i,j,k) - 0.5e0*gup23(i,j,k)*g133(i,j,k)
               c222(i,j,k)= c222(i,j,k) - 0.5e0*gup23(i,j,k)*g223(i,j,k)
               c223(i,j,k)= c223(i,j,k) - 0.5e0*gup23(i,j,k)*g233(i,j,k)
               c233(i,j,k)= c233(i,j,k) - 0.5e0*gup23(i,j,k)*g333(i,j,k)

               c311(i,j,k)= c311(i,j,k) - 0.5e0*gup33(i,j,k)*g113(i,j,k)
               c312(i,j,k)= c312(i,j,k) - 0.5e0*gup33(i,j,k)*g123(i,j,k)
               c313(i,j,k)= c313(i,j,k) - 0.5e0*gup33(i,j,k)*g133(i,j,k)
               c322(i,j,k)= c322(i,j,k) - 0.5e0*gup33(i,j,k)*g223(i,j,k)
               c323(i,j,k)= c323(i,j,k) - 0.5e0*gup33(i,j,k)*g233(i,j,k)
               c333(i,j,k)= c333(i,j,k) - 0.5e0*gup33(i,j,k)*g333(i,j,k)

            end do
         end do
      end do

      return
      end
