c-----------------------------------------------------------------------
*     Randall R. Correll
*     Mijan F. Huq
*     Scott A. Klasky
*     Richard A. Matzner
*
*     COPYRIGHT 1995-1996

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

      subroutine riccidd(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, 
     *           m,d,nl, h, nx, ny, nz)
        implicit none

        integer  nx, ny, nz, nl

        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 m(3*nl)

        real*8 d(3*nl)
        real*8   h

c		call subsidiary routines to calculate ricci tensor piecewise.

        call ricci1(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, nx, ny, nz)
        call ricci2(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, nx, ny, nz)
        call ricci3(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, nx, ny, nz)

        call ricci4(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, nx, ny, nz)

        call ricci5(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, nx, ny, nz)

        call ricci6(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, nx, ny, nz)
        call 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, nx, ny, nz)

        return 
      end
c-----------------------------------------------------------------------
c  even though c{ijk}d{l} is symmetric in jk we keep other terms for
c  clean interfaces b/w cases where we have non-symmetric components later
c
c  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, nx, ny, nz)
        implicit none

        integer     nx, ny, nz, nl

        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   c111d1(nx,ny,nz), c112d1(nx,ny,nz), c113d1(nx,ny,nz), 
     *           c121d1(nx,ny,nz), c122d1(nx,ny,nz), c123d1(nx,ny,nz), 
     *           c131d1(nx,ny,nz), c132d1(nx,ny,nz), c133d1(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)

        real*8   h
        integer m(3*nl)

        real*8 d(3*nl)

c    Local Variables

        integer  i, j, k

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

        call diff(nx, ny, nz, h, c111, c111d1, 'x',m,d,nl)

        call diff(nx, ny, nz, h, c112, c112d1, 'x',m,d,nl)

        call diff(nx, ny, nz, h, c113, c113d1, 'x',m,d,nl)

        call diff(nx, ny, nz, h, c122, c122d1, 'x',m,d,nl)

        call diff(nx, ny, nz, h, c123, c123d1, 'x',m,d,nl)

        call diff(nx, ny, nz, h, c133, c133d1, 'x',m,d,nl)

c		Evaluate first term 
        do k = 1, nz
          do j = 1, ny
            do i = 1, nx
              r11(i,j,k) = c111d1(i,j,k)
              r12(i,j,k) = c112d1(i,j,k)
              r13(i,j,k) = c113d1(i,j,k)
              r22(i,j,k) = c122d1(i,j,k)
              r23(i,j,k) = c123d1(i,j,k)
              r33(i,j,k) = c133d1(i,j,k)
            end do
          end do
         end do

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

        integer     nx, ny, nz, nl

        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   c211d2(nx,ny,nz), c212d2(nx,ny,nz), c213d2(nx,ny,nz), 
     *           c221d2(nx,ny,nz), c222d2(nx,ny,nz), c223d2(nx,ny,nz), 
     *           c231d2(nx,ny,nz), c232d2(nx,ny,nz), c233d2(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)

        real*8   h
        integer m(3*nl)

        real*8 d(3*nl)

c    Local Variables

        integer  i, j, k

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

        call diff(nx, ny, nz, h, c211, c211d2, 'y',m,d,nl)

        call diff(nx, ny, nz, h, c212, c212d2, 'y',m,d,nl)

        call diff(nx, ny, nz, h, c213, c213d2, 'y',m,d,nl)

        call diff(nx, ny, nz, h, c222, c222d2, 'y',m,d,nl)

        call diff(nx, ny, nz, h, c223, c223d2, 'y',m,d,nl)

        call diff(nx, ny, nz, h, c233, c233d2, 'y',m,d,nl)

c		Evaluate second term 
        do k = 1, nz
          do j = 1, ny
            do i = 1, nx
              r11(i,j,k) = r11(i,j,k) + c211d2(i,j,k)
              r12(i,j,k) = r12(i,j,k) + c212d2(i,j,k)
              r13(i,j,k) = r13(i,j,k) + c213d2(i,j,k)
              r22(i,j,k) = r22(i,j,k) + c222d2(i,j,k)
              r23(i,j,k) = r23(i,j,k) + c223d2(i,j,k)
              r33(i,j,k) = r33(i,j,k) + c233d2(i,j,k)
            end do
          end do
         end do

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

        integer     nx, ny, nz, nl

        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   c311d3(nx,ny,nz), c312d3(nx,ny,nz), c313d3(nx,ny,nz), 
     *           c321d3(nx,ny,nz), c322d3(nx,ny,nz), c323d3(nx,ny,nz), 
     *           c331d3(nx,ny,nz), c332d3(nx,ny,nz), c333d3(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)

        real*8   h
        integer m(3*nl)

        real*8 d(3*nl)

c    Local Variables

        integer  i, j, k

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

        call diff(nx, ny, nz, h, c311, c311d3, 'z',m,d,nl)

        call diff(nx, ny, nz, h, c312, c312d3, 'z',m,d,nl)

        call diff(nx, ny, nz, h, c313, c313d3, 'z',m,d,nl)

        call diff(nx, ny, nz, h, c322, c322d3, 'z',m,d,nl)

        call diff(nx, ny, nz, h, c323, c323d3, 'z',m,d,nl)

        call diff(nx, ny, nz, h, c333, c333d3, 'z',m,d,nl)

c		Evaluate third term 
        do k = 1, nz
          do j = 1, ny
            do i = 1, nx
              r11(i,j,k) = r11(i,j,k) + c311d3(i,j,k)
              r12(i,j,k) = r12(i,j,k) + c312d3(i,j,k)
              r13(i,j,k) = r13(i,j,k) + c313d3(i,j,k)
              r22(i,j,k) = r22(i,j,k) + c322d3(i,j,k)
              r23(i,j,k) = r23(i,j,k) + c323d3(i,j,k)
              r33(i,j,k) = r33(i,j,k) + c333d3(i,j,k)
            end do
          end do
         end do

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

        integer     nx, ny, nz, nl

        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   c111d1(nx,ny,nz), c111d2(nx,ny,nz), c111d3(nx,ny,nz), 
     *           c121d1(nx,ny,nz), c121d2(nx,ny,nz), c121d3(nx,ny,nz), 
     *           c131d1(nx,ny,nz), c131d2(nx,ny,nz), c131d3(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)

        real*8   h
        integer m(3*nl)

        real*8 d(3*nl)

c    Local Variables

        integer  i, j, k

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

        call diff(nx, ny, nz, h, c111, c111d1, 'x',m,d,nl)

        call diff(nx, ny, nz, h, c111, c111d2, 'y',m,d,nl)

        call diff(nx, ny, nz, h, c111, c111d3, 'z',m,d,nl)

c        call diff(nx, ny, nz, h, c121, c121d1, 'x',m,d,nl)

        call diff(nx, ny, nz, h, c112, c121d2, 'y',m,d,nl)

        call diff(nx, ny, nz, h, c112, c121d3, 'z',m,d,nl)

c        call diff(nx, ny, nz, h, c131, c131d1, 'x',m,d,nl)

c        call diff(nx, ny, nz, h, c131, c131d2, 'y',m,d,nl)

        call diff(nx, ny, nz, h, c113, c131d3, 'z',m,d,nl)

c		Evaluate first term 
        do k = 1, nz
          do j = 1, ny
            do i = 1, nx
              r11(i,j,k) = r11(i,j,k) - c111d1(i,j,k)
              r12(i,j,k) = r12(i,j,k) - c111d2(i,j,k)
              r13(i,j,k) = r13(i,j,k) - c111d3(i,j,k)
              r22(i,j,k) = r22(i,j,k) - c121d2(i,j,k)
              r23(i,j,k) = r23(i,j,k) - c121d3(i,j,k)
              r33(i,j,k) = r33(i,j,k) - c131d3(i,j,k)
            end do
          end do
         end do

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

        integer     nx, ny, nz, nl

        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   c212d1(nx,ny,nz), c212d2(nx,ny,nz), c212d3(nx,ny,nz), 
     *           c222d1(nx,ny,nz), c222d2(nx,ny,nz), c222d3(nx,ny,nz), 
     *           c232d1(nx,ny,nz), c232d2(nx,ny,nz), c232d3(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)

        real*8   h

        integer m(3*nl)

        real*8 d(3*nl)
c    Local Variables

        integer  i, j, k

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

        call diff(nx, ny, nz, h, c212, c212d1, 'x',m,d,nl)

        call diff(nx, ny, nz, h, c212, c212d2, 'y',m,d,nl)

        call diff(nx, ny, nz, h, c212, c212d3, 'z',m,d,nl)

c        call diff(nx, ny, nz, h, c222, c222d1, 'x',m,d,nl)

        call diff(nx, ny, nz, h, c222, c222d2, 'y',m,d,nl)

        call diff(nx, ny, nz, h, c222, c222d3, 'z',m,d,nl)

c        call diff(nx, ny, nz, h, c232, c232d1, 'x',m,d,nl)

c        call diff(nx, ny, nz, h, c232, c232d2, 'y',m,d,nl)

        call diff(nx, ny, nz, h, c223, c232d3, 'z',m,d,nl)

c		Evaluate first term 
        do k = 1, nz
          do j = 1, ny
            do i = 1, nx
              r11(i,j,k) = r11(i,j,k) - c212d1(i,j,k)
              r12(i,j,k) = r12(i,j,k) - c212d2(i,j,k)
              r13(i,j,k) = r13(i,j,k) - c212d3(i,j,k)
              r22(i,j,k) = r22(i,j,k) - c222d2(i,j,k)
              r23(i,j,k) = r23(i,j,k) - c222d3(i,j,k)
              r33(i,j,k) = r33(i,j,k) - c232d3(i,j,k)
            end do
          end do
         end do

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

        integer     nx, ny, nz, nl

        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   c313d1(nx,ny,nz), c313d2(nx,ny,nz), c313d3(nx,ny,nz), 
     *           c323d1(nx,ny,nz), c323d2(nx,ny,nz), c323d3(nx,ny,nz), 
     *           c333d1(nx,ny,nz), c333d2(nx,ny,nz), c333d3(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)

        real*8   h
        integer m(3*nl)

        real*8 d(3*nl)

c    Local Variables

        integer  i, j, k

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

        call diff(nx, ny, nz, h, c313, c313d1, 'x',m,d,nl)

        call diff(nx, ny, nz, h, c313, c313d2, 'y',m,d,nl)

        call diff(nx, ny, nz, h, c313, c313d3, 'z',m,d,nl)

        call diff(nx, ny, nz, h, c323, c323d2, 'y',m,d,nl)

        call diff(nx, ny, nz, h, c323, c323d3, 'z',m,d,nl)

        call diff(nx, ny, nz, h, c333, c333d3, 'z',m,d,nl)

c		Evaluate first term 
        do k = 1, nz
          do j = 1, ny
            do i = 1, nx
              r11(i,j,k) = r11(i,j,k) - c313d1(i,j,k)
              r12(i,j,k) = r12(i,j,k) - c313d2(i,j,k)
              r13(i,j,k) = r13(i,j,k) - c313d3(i,j,k)
              r22(i,j,k) = r22(i,j,k) - c323d2(i,j,k)
              r23(i,j,k) = r23(i,j,k) - c323d3(i,j,k)
              r33(i,j,k) = r33(i,j,k) - c333d3(i,j,k)
            end do
          end do
         end do

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

        integer     nx, ny, nz, nl

        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)

        real*8   h

        integer m(3*nl)

        real*8 d(3*nl)
c    Local Variables

        integer  i, j, k

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

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

c              rij(i,j,k) = rij(i,j,k) + 
c     *          (c111(i,j,k) + c212(i,j,k) + c313(i,j,k))*c1ij(i,j,k) +
c     *          (c112(i,j,k) + c222(i,j,k) + c323(i,j,k))*c2ij(i,j,k) +
c     *          (c113(i,j,k) + c223(i,j,k) + c333(i,j,k))*c3ij(i,j,k)

              r11(i,j,k) = r11(i,j,k) + 
     *          (c111(i,j,k) + c212(i,j,k) + c313(i,j,k))*c111(i,j,k) +
     *          (c112(i,j,k) + c222(i,j,k) + c323(i,j,k))*c211(i,j,k) +
     *          (c113(i,j,k) + c223(i,j,k) + c333(i,j,k))*c311(i,j,k)

              r12(i,j,k) = r12(i,j,k) + 
     *          (c111(i,j,k) + c212(i,j,k) + c313(i,j,k))*c112(i,j,k) +
     *          (c112(i,j,k) + c222(i,j,k) + c323(i,j,k))*c212(i,j,k) +
     *          (c113(i,j,k) + c223(i,j,k) + c333(i,j,k))*c312(i,j,k)

              r13(i,j,k) = r13(i,j,k) + 
     *          (c111(i,j,k) + c212(i,j,k) + c313(i,j,k))*c113(i,j,k) +
     *          (c112(i,j,k) + c222(i,j,k) + c323(i,j,k))*c213(i,j,k) +
     *          (c113(i,j,k) + c223(i,j,k) + c333(i,j,k))*c313(i,j,k)

              r22(i,j,k) = r22(i,j,k) + 
     *          (c111(i,j,k) + c212(i,j,k) + c313(i,j,k))*c122(i,j,k) +
     *          (c112(i,j,k) + c222(i,j,k) + c323(i,j,k))*c222(i,j,k) +
     *          (c113(i,j,k) + c223(i,j,k) + c333(i,j,k))*c322(i,j,k)

              r23(i,j,k) = r23(i,j,k) + 
     *          (c111(i,j,k) + c212(i,j,k) + c313(i,j,k))*c123(i,j,k) +
     *          (c112(i,j,k) + c222(i,j,k) + c323(i,j,k))*c223(i,j,k) +
     *          (c113(i,j,k) + c223(i,j,k) + c333(i,j,k))*c323(i,j,k)

              r33(i,j,k) = r33(i,j,k) + 
     *          (c111(i,j,k) + c212(i,j,k) + c313(i,j,k))*c133(i,j,k) +
     *          (c112(i,j,k) + c222(i,j,k) + c323(i,j,k))*c233(i,j,k) +
     *          (c113(i,j,k) + c223(i,j,k) + c333(i,j,k))*c333(i,j,k)

            end do
          end do
        end do


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

        do k = 1, nz
          do j = 1, ny
            do i = 1, nx
c              rij(i,j,k)  = rij(i,j,k) - ( c1i1(i,j,k) * c11j(i,j,k) + 
c     *          c2i1(i,j,k) * c12j(i,j,k) + c3i1(i,j,k)*c13j(i,j,k)) -
c     *         ( c1i2(i,j,k) * c21j(i,j,k) + c2i2(i,j,k) * c22j(i,j,k) + 
c     *         c3i2(i,j,k)*c23j(i,j,k)) - ( c1i3(i,j,k) * c31j(i,j,k) + 
c     *         c2i3(i,j,k) * c32j(i,j,k) + c3i3(i,j,k)*c33j(i,j,k))

              r11(i,j,k)  = r11(i,j,k) - ( c111(i,j,k) * c111(i,j,k) + 
     *          c211(i,j,k) * c112(i,j,k) + c311(i,j,k)*c113(i,j,k)) -
     *         (c112(i,j,k) * c211(i,j,k) + c212(i,j,k) * c212(i,j,k) + 
     *         c312(i,j,k)*c213(i,j,k)) - ( c113(i,j,k) * c311(i,j,k) + 
     *         c213(i,j,k) * c312(i,j,k) + c313(i,j,k)*c313(i,j,k))

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

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

              r22(i,j,k)  = r22(i,j,k) - ( c112(i,j,k) * c112(i,j,k) + 
     *          c212(i,j,k) * c122(i,j,k) + c312(i,j,k)*c123(i,j,k)) -
     *         ( c122(i,j,k) * c212(i,j,k) + c222(i,j,k) * c222(i,j,k) + 
     *         c322(i,j,k)*c223(i,j,k)) - ( c123(i,j,k) * c312(i,j,k) + 
     *         c223(i,j,k) * c322(i,j,k) + c323(i,j,k)*c323(i,j,k))

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

              r33(i,j,k)  = r33(i,j,k) - ( c113(i,j,k) * c113(i,j,k) + 
     *          c213(i,j,k) * c123(i,j,k) + c313(i,j,k)*c133(i,j,k)) -
     *         ( c123(i,j,k) * c213(i,j,k) + c223(i,j,k) * c223(i,j,k) + 
     *         c323(i,j,k)*c233(i,j,k)) - ( c133(i,j,k) * c313(i,j,k) + 
     *         c233(i,j,k) * c323(i,j,k) + c333(i,j,k)*c333(i,j,k))

            end do
          end do
        end do

        return
      end
c-----------------------------------------------------------------------
