c-----------------------------------------------------------------------
c evolve_k.f
*     Randall R. Correll
*     Mijan F. Huq
*     Scott A. Klasky
*     Richard A. Matzner
*
*     COPYRIGHT 1995-1996
c 
c ADDED temporary variables temp1, temp2, temp3
c-----------------------------------------------------------------------
      subroutine evolve_k(g11, g12, g13, g22, g23, g33,
     $           gup11, gup12, gup13, gup22, gup23, gup33,
     $           k11m1,k12m1,k13m1,
     $           k21m1,k22m1,k23m1,
     $           k31m1,k32m1,k33m1,
     $           k11p1,k12p1,k13p1,
     $           k21p1,k22p1,k23p1,
     $           k31p1,k32p1,k33p1,
     $           k11, k12, k13,
     $           k21, k22, k23,
     $           k31, k32, k33,
     $           c111, c112, c113, c122, c123, c133,
     $           c211, c212, c213, c222, c223, c233,
     $           c311, c312, c313, c322, c323, c333,
     $           beta1, beta2, beta3, alpha,
     $           t11, t12, t13, t21, t22, t23, t31, t32, t33,
     $           temp1, temp2, temp3,
     $           r11, r12, r13, r22, r23, r33,
     $           m,d,nl, h, dt,nx, ny, nz)
      implicit none

      integer  nx, ny, nz, nl
      integer  i, j, k

      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)

      real*8   k11m1(nx,ny,nz), k12m1(nx,ny,nz), k13m1(nx,ny,nz),
     $         k21m1(nx,ny,nz), k22m1(nx,ny,nz), k23m1(nx,ny,nz),
     $         k31m1(nx,ny,nz), k32m1(nx,ny,nz), k33m1(nx,ny,nz)

      real*8   k11p1(nx,ny,nz), k12p1(nx,ny,nz), k13p1(nx,ny,nz),
     $         k21p1(nx,ny,nz), k22p1(nx,ny,nz), k23p1(nx,ny,nz),
     $         k31p1(nx,ny,nz), k32p1(nx,ny,nz), k33p1(nx,ny,nz)

      real*8   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   k11(nx,ny,nz), k12(nx,ny,nz), k13(nx,ny,nz),
     $         k21(nx,ny,nz), k22(nx,ny,nz), k23(nx,ny,nz),
     $         k31(nx,ny,nz), k32(nx,ny,nz), k33(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   t11(nx,ny,nz),t12(nx,ny,nz),t13(nx,ny,nz),t21(nx,ny,nz),
     $         t22(nx,ny,nz),t23(nx,ny,nz),t31(nx,ny,nz),t32(nx,ny,nz),
     $         t33(nx,ny,nz),
     $         temp1(nx,ny,nz), temp2(nx,ny,nz), temp3(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 beta1(nx,ny,nz),beta2(nx,ny,nz),beta3(nx,ny,nz)
      real*8 alpha(nx,ny,nz)
      integer m(3*nl)

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

c     Compute D^a(D_b(alpha)) terms

		 call evolve_k0(gup11,gup12,gup13,gup22,gup23,gup33,
     $                     c111,c112,c113,c122,c123,c133,
     $                     c211,c212,c213,c222,c223,c233,
     $                     c311,c312,c313,c322,c323,c333,
     $                     alpha,
     $                     t11,t12,t13,
     $                     t21,t22,t23,
     $                     t31,t32,t33,temp1,temp2,temp3,
     $                     m,d,nl,h,nx,ny,nz)

		 call evolve_k1(gup11,gup12,gup13,gup22,gup23,gup33,
     $                     c111,c112,c113,c122,c123,c133,
     $                     c211,c212,c213,c222,c223,c233,
     $                     c311,c312,c313,c322,c323,c333,
     $                     alpha,
     $                     t11,t12,t13,
     $                     t21,t22,t23,
     $                     t31,t32,t33,temp1,temp2,temp3,
     $                     m,d,nl,h,nx,ny,nz)

c     Compute Lie_beta(K^a_b) terms
      call evolve_k2(gup11,gup12,gup13,gup22,gup23,gup33,
     $               k11,k12,k13,
     $               k21,k22,k23,
     $               k31,k32,k33,
     $               beta1,beta2,beta3,
     $               t11,t12,t13,
     $               t21,t22,t23,
     $               t31,t32,t33,temp1,
     $               m,d,nl,h,nz,ny,nz)

c     Compute alpha*(R^a_b + K*K^a_b) terms
      call evolve_k3(gup11,gup12,gup13,gup22,gup23,gup33,
     $               r11,r12,r13,r22,r23,r33,
     $               k11,k12,k13,
     $               k21,k22,k23,
     $               k31,k32,k33,
     $               alpha,
     $               t11,t12,t13,
     $               t21,t22,t23,
     $               t31,t32,t33,temp1,
     $               nx,ny,nz)


c     Evolve one timestep via leapfrog method ( i.e. time-centered)
      call leapfrog(nx,ny,nz,k11m1,k11p1,t11,dt)
      call leapfrog(nx,ny,nz,k12m1,k12p1,t12,dt)
      call leapfrog(nx,ny,nz,k13m1,k13p1,t13,dt)
      call leapfrog(nx,ny,nz,k21m1,k21p1,t21,dt)
      call leapfrog(nx,ny,nz,k22m1,k22p1,t22,dt)
      call leapfrog(nx,ny,nz,k23m1,k23p1,t23,dt)
      call leapfrog(nx,ny,nz,k31m1,k31p1,t31,dt)
      call leapfrog(nx,ny,nz,k32m1,k32p1,t32,dt)
      call leapfrog(nx,ny,nz,k33m1,k33p1,t33,dt)

      return
      end

************************************************************************
*     Computes D^a(D_b(alpha)) terms
*     INCLUDED tempd in calling sequence.
*     REMEMBER no variables to be initialized inside subroutines. DAGH
*              structure requires us to pass in arrays from outside...
************************************************************************
      subroutine evolve_k0(gup11,gup12,gup13,gup22,gup23,gup33,
     $                     c111,c112,c113,c122,c123,c133,
     $                     c211,c212,c213,c222,c223,c233,
     $                     c311,c312,c313,c322,c323,c333,
     $                     alpha,
     $                     t11,t12,t13,
     $                     t21,t22,t23,
     $                     t31,t32,t33,temp1,temp2,temp3,
     $                     m,d,nl,h,nx,ny,nz)

      implicit none

      integer  nx, ny, nz, nl
      real*8   h

      real*8   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   alpha(nx,ny,nz)

      real*8   t11(nx,ny,nz),t12(nx,ny,nz),t13(nx,ny,nz),
     $         t21(nx,ny,nz),t22(nx,ny,nz),t23(nx,ny,nz),
     $         t31(nx,ny,nz),t32(nx,ny,nz),t33(nx,ny,nz),
     $         temp1(nx,ny,nz),temp2(nx,ny,nz),temp3(nx,ny,nz)

      integer  m(3*nl)

      real*8   d(3*nl)

      integer	i, j, k
C
c	Calculate first derivatives of alpha
      call diff(nx,ny,nz,h,alpha,temp1,'x',m,d,nl)
      call diff(nx,ny,nz,h,alpha,temp2,'y',m,d,nl)
      call diff(nx,ny,nz,h,alpha,temp3,'z',m,d,nl)
		return
		end

************************************************************************
*     Computes D^a(D_b(alpha)) terms
*     INCLUDED tempd in calling sequence.
*     REMEMBER no variables to be initialized inside subroutines. DAGH
*              structure requires us to pass in arrays from outside...
************************************************************************
      subroutine evolve_k1(gup11,gup12,gup13,gup22,gup23,gup33,
     $                     c111,c112,c113,c122,c123,c133,
     $                     c211,c212,c213,c222,c223,c233,
     $                     c311,c312,c313,c322,c323,c333,
     $                     alpha,
     $                     t11,t12,t13,
     $                     t21,t22,t23,
     $                     t31,t32,t33,temp1,temp2,temp3,
     $                     m,d,nl,h,nx,ny,nz)

      implicit none

      integer  nx, ny, nz, nl
      real*8   h

      real*8   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   alpha(nx,ny,nz)

      real*8   t11(nx,ny,nz),t12(nx,ny,nz),t13(nx,ny,nz),
     $         t21(nx,ny,nz),t22(nx,ny,nz),t23(nx,ny,nz),
     $         t31(nx,ny,nz),t32(nx,ny,nz),t33(nx,ny,nz),
     $         temp1(nx,ny,nz),temp2(nx,ny,nz),temp3(nx,ny,nz)

      integer  m(3*nl)

      real*8   d(3*nl)

      integer	i, j, k

c	Calculate second derivatives of alpha

      call diff(nx,ny,nz,h,temp1,t11,'x',m,d,nl)
      call diff(nx,ny,nz,h,temp1,t12,'y',m,d,nl)
      call diff(nx,ny,nz,h,temp1,t13,'z',m,d,nl)

		call copy(nx, ny, nz, t12, t21)

      call diff(nx,ny,nz,h,temp2,t22,'y',m,d,nl)
      call diff(nx,ny,nz,h,temp2,t23,'z',m,d,nl)

		call copy(nx, ny, nz, t13, t31)

		call copy(nx, ny, nz, t23, t32)

      call diff(nx,ny,nz,h,temp3,t33,'z',m,d,nl)

c form covariant derivative D_a D_b alpha
c D_a D_b alpha = t{ab} - c1{ab}*temp1 - c2{ab}*temp2 - c3{ab}*temp3

c   11
		do k = 1, nz
		  do j = 1, ny
		    do i = 1, nx
		      t11(i,j,k) = t11(i,j,k) - c111(i,j,k) * temp1(i,j,k)
     *                   - c211(i,j,k) * temp2(i,j,k)
     *                   - c311(i,j,k) * temp3(i,j,k)
		    end do
		  end do
		end do

c   12
		do k = 1, nz
		  do j = 1, ny
		    do i = 1, nx
		      t12(i,j,k) = t12(i,j,k) - c112(i,j,k) * temp1(i,j,k)
     *                   - c212(i,j,k) * temp2(i,j,k)
     *                   - c312(i,j,k) * temp3(i,j,k)
		    end do
		  end do
		end do

c   13
		do k = 1, nz
		  do j = 1, ny
		    do i = 1, nx
		      t13(i,j,k) = t13(i,j,k) - c113(i,j,k) * temp1(i,j,k)
     *                   - c213(i,j,k) * temp2(i,j,k)
     *                   - c313(i,j,k) * temp3(i,j,k)
		    end do
		  end do
		end do

c	21
c			copy since symmetric in {ab}

		call copy(nx, ny, nz, t12, t21)

c   22
		do k = 1, nz
		  do j = 1, ny
		    do i = 1, nx
		      t22(i,j,k) = t22(i,j,k) - c122(i,j,k) * temp1(i,j,k)
     *                   - c222(i,j,k) * temp2(i,j,k)
     *                   - c322(i,j,k) * temp3(i,j,k)
		    end do
		  end do
		end do

c   23
		do k = 1, nz
		  do j = 1, ny
		    do i = 1, nx
		      t23(i,j,k) = t23(i,j,k) - c123(i,j,k) * temp1(i,j,k)
     *                   - c223(i,j,k) * temp2(i,j,k)
     *                   - c323(i,j,k) * temp3(i,j,k)
		    end do
		  end do
		end do

c	31
		call copy(nx, ny, nz, t13, t31)

c	32
		call copy(nx, ny, nz, t23, t32)

c   33
		do k = 1, nz
		  do j = 1, ny
		    do i = 1, nx
		      t33(i,j,k) = t33(i,j,k) - c133(i,j,k) * temp1(i,j,k)
     *                   - c233(i,j,k) * temp2(i,j,k)
     *                   - c333(i,j,k) * temp3(i,j,k)
		    end do
		  end do
		end do

c	raise the first index of t{ab}

		call rais02i1(nx, ny, nz,
     *           gup11, gup12, gup13, gup22, gup23, gup33,
     *           t11, t12, t13, t21, t22, t23, t31, t32, t33,
     *           temp1, temp2, temp3)


c mulitply by -1 since term is negative in equation
		call ms_mult(nx,ny,nz,-1.0e0,t11,t11)
		call ms_mult(nx,ny,nz,-1.0e0,t12,t12)
		call ms_mult(nx,ny,nz,-1.0e0,t13,t13)

		call ms_mult(nx,ny,nz,-1.0e0,t21,t21)
		call ms_mult(nx,ny,nz,-1.0e0,t22,t22)
		call ms_mult(nx,ny,nz,-1.0e0,t23,t23)

		call ms_mult(nx,ny,nz,-1.0e0,t31,t31)
		call ms_mult(nx,ny,nz,-1.0e0,t32,t32)
		call ms_mult(nx,ny,nz,-1.0e0,t33,t33)
		
      
      return
      end

************************************************************************
*     Computes Lie_beta(K^a_b) terms
************************************************************************
      subroutine evolve_k2(gup11,gup12,gup13,gup22,gup23,gup33,
     $                     k11,k12,k13,
     $                     k21,k22,k23,
     $                     k31,k32,k33,
     $                     beta1,beta2,beta3,
     $                     t11,t12,t13,
     $                     t21,t22,t23,
     $                     t31,t32,t33,temp,
     $                     m,d,nl,h,nx,ny,nz)

      implicit none

      integer  nx, ny, nz, nl
      real*8   h

      real*8   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   k11(nx,ny,nz), k12(nx,ny,nz), k13(nx,ny,nz),
     $         k21(nx,ny,nz), k22(nx,ny,nz), k23(nx,ny,nx),
     $         k31(nx,ny,nz), k32(nx,ny,nz), k33(nx,ny,nz)

      real*8   beta1(nx,ny,nz), beta2(nx,ny,nz), beta3(nx,ny,nz)

      real*8   t11(nx,ny,nz),t12(nx,ny,nz),t13(nx,ny,nz),
     $         t21(nx,ny,nz),t22(nx,ny,nz),t23(nx,ny,nz),
     $         t31(nx,ny,nz),t32(nx,ny,nz),t33(nx,ny,nz),
     $         temp(nx,ny,nz)

      integer  m(3*nl)

      real*8   d(3*nl)

*     Computes Lie_beta(K^a_b) terms

*     Lie_beta(K^1_1)
      call diff(nx,ny,nz,h,k11,temp,'x',m,d,nl)
      call mm_mult(nx,ny,nz,beta1,temp,temp)
      call mm_add(nx,ny,nz,t11,temp,t11)
      call diff(nx,ny,nz,h,k11,temp,'y',m,d,nl)
      call mm_mult(nx,ny,nz,beta2,temp,temp)
      call mm_add(nx,ny,nz,t11,temp,t11)
      call diff(nx,ny,nz,h,k11,temp,'z',m,d,nl)
      call mm_mult(nx,ny,nz,beta3,temp,temp)
      call mm_add(nx,ny,nz,t11,temp,t11)
      call diff(nx,ny,nz,h,beta1,temp,'x',m,d,nl)
      call mm_mult(nx,ny,nz,k11,temp,temp)
      call mm_sub(nx,ny,nz,t11,temp,t11)
      call diff(nx,ny,nz,h,beta1,temp,'y',m,d,nl)
      call mm_mult(nx,ny,nz,k21,temp,temp)
      call mm_sub(nx,ny,nz,t11,temp,t11)
      call diff(nx,ny,nz,h,beta1,temp,'z',m,d,nl)
      call mm_mult(nx,ny,nz,k31,temp,temp)
      call mm_sub(nx,ny,nz,t11,temp,t11)
      call diff(nx,ny,nz,h,beta1,temp,'x',m,d,nl)
      call mm_mult(nx,ny,nz,k11,temp,temp)
      call mm_add(nx,ny,nz,t11,temp,t11)
      call diff(nx,ny,nz,h,beta2,temp,'x',m,d,nl)
      call mm_mult(nx,ny,nz,k12,temp,temp)
      call mm_add(nx,ny,nz,t11,temp,t11)
      call diff(nx,ny,nz,h,beta3,temp,'x',m,d,nl)
      call mm_mult(nx,ny,nz,k13,temp,temp)
      call mm_add(nx,ny,nz,t11,temp,t11)

*     Lie_beta(K^1_2)
      call diff(nx,ny,nz,h,k12,temp,'x',m,d,nl)
      call mm_mult(nx,ny,nz,beta1,temp,temp)
      call mm_add(nx,ny,nz,t12,temp,t12)
      call diff(nx,ny,nz,h,k12,temp,'y',m,d,nl)
      call mm_mult(nx,ny,nz,beta2,temp,temp)
      call mm_add(nx,ny,nz,t12,temp,t12)
      call diff(nx,ny,nz,h,k12,temp,'z',m,d,nl)
      call mm_mult(nx,ny,nz,beta3,temp,temp)
      call mm_add(nx,ny,nz,t12,temp,t12)
      call diff(nx,ny,nz,h,beta1,temp,'x',m,d,nl)
      call mm_mult(nx,ny,nz,k12,temp,temp)
      call mm_sub(nx,ny,nz,t12,temp,t12)
      call diff(nx,ny,nz,h,beta1,temp,'y',m,d,nl)
      call mm_mult(nx,ny,nz,k22,temp,temp)
      call mm_sub(nx,ny,nz,t12,temp,t12)
      call diff(nx,ny,nz,h,beta1,temp,'z',m,d,nl)
      call mm_mult(nx,ny,nz,k32,temp,temp)
      call mm_sub(nx,ny,nz,t12,temp,t12)
      call diff(nx,ny,nz,h,beta1,temp,'y',m,d,nl)
      call mm_mult(nx,ny,nz,k11,temp,temp)
      call mm_add(nx,ny,nz,t12,temp,t12)
      call diff(nx,ny,nz,h,beta2,temp,'y',m,d,nl)
      call mm_mult(nx,ny,nz,k12,temp,temp)
      call mm_add(nx,ny,nz,t12,temp,t12)
      call diff(nx,ny,nz,h,beta3,temp,'y',m,d,nl)
      call mm_mult(nx,ny,nz,k13,temp,temp)
      call mm_add(nx,ny,nz,t12,temp,t12)

*     Lie_beta(K^1_3)
      call diff(nx,ny,nz,h,k13,temp,'x',m,d,nl)
      call mm_mult(nx,ny,nz,beta1,temp,temp)
      call mm_add(nx,ny,nz,t13,temp,t13)
      call diff(nx,ny,nz,h,k13,temp,'y',m,d,nl)
      call mm_mult(nx,ny,nz,beta2,temp,temp)
      call mm_add(nx,ny,nz,t13,temp,t13)
      call diff(nx,ny,nz,h,k13,temp,'z',m,d,nl)
      call mm_mult(nx,ny,nz,beta3,temp,temp)
      call mm_add(nx,ny,nz,t13,temp,t13)
      call diff(nx,ny,nz,h,beta1,temp,'x',m,d,nl)
      call mm_mult(nx,ny,nz,k13,temp,temp)
      call mm_sub(nx,ny,nz,t13,temp,t13)
      call diff(nx,ny,nz,h,beta1,temp,'y',m,d,nl)
      call mm_mult(nx,ny,nz,k23,temp,temp)
      call mm_sub(nx,ny,nz,t13,temp,t13)
      call diff(nx,ny,nz,h,beta1,temp,'z',m,d,nl)
      call mm_mult(nx,ny,nz,k33,temp,temp)
      call mm_sub(nx,ny,nz,t13,temp,t13)
      call diff(nx,ny,nz,h,beta1,temp,'z',m,d,nl)
      call mm_mult(nx,ny,nz,k11,temp,temp)
      call mm_add(nx,ny,nz,t13,temp,t13)
      call diff(nx,ny,nz,h,beta2,temp,'z',m,d,nl)
      call mm_mult(nx,ny,nz,k12,temp,temp)
      call mm_add(nx,ny,nz,t13,temp,t13)
      call diff(nx,ny,nz,h,beta3,temp,'z',m,d,nl)
      call mm_mult(nx,ny,nz,k13,temp,temp)
      call mm_add(nx,ny,nz,t13,temp,t13)

*     Lie_beta(K^2_1)
      call diff(nx,ny,nz,h,k21,temp,'x',m,d,nl)
      call mm_mult(nx,ny,nz,beta1,temp,temp)
      call mm_add(nx,ny,nz,t21,temp,t21)
      call diff(nx,ny,nz,h,k21,temp,'y',m,d,nl)
      call mm_mult(nx,ny,nz,beta2,temp,temp)
      call mm_add(nx,ny,nz,t21,temp,t21)
      call diff(nx,ny,nz,h,k21,temp,'z',m,d,nl)
      call mm_mult(nx,ny,nz,beta3,temp,temp)
      call mm_add(nx,ny,nz,t21,temp,t21)
      call diff(nx,ny,nz,h,beta2,temp,'x',m,d,nl)
      call mm_mult(nx,ny,nz,k11,temp,temp)
      call mm_sub(nx,ny,nz,t21,temp,t21)
      call diff(nx,ny,nz,h,beta2,temp,'y',m,d,nl)
      call mm_mult(nx,ny,nz,k21,temp,temp)
      call mm_sub(nx,ny,nz,t21,temp,t21)
      call diff(nx,ny,nz,h,beta2,temp,'z',m,d,nl)
      call mm_mult(nx,ny,nz,k31,temp,temp)
      call mm_sub(nx,ny,nz,t21,temp,t21)
      call diff(nx,ny,nz,h,beta1,temp,'x',m,d,nl)
      call mm_mult(nx,ny,nz,k21,temp,temp)
      call mm_add(nx,ny,nz,t21,temp,t21)
      call diff(nx,ny,nz,h,beta2,temp,'x',m,d,nl)
      call mm_mult(nx,ny,nz,k22,temp,temp)
      call mm_add(nx,ny,nz,t21,temp,t21)
      call diff(nx,ny,nz,h,beta3,temp,'x',m,d,nl)
      call mm_mult(nx,ny,nz,k23,temp,temp)
      call mm_add(nx,ny,nz,t21,temp,t21)

*     Lie_beta(K^2_2)
      call diff(nx,ny,nz,h,k22,temp,'x',m,d,nl)
      call mm_mult(nx,ny,nz,beta1,temp,temp)
      call mm_add(nx,ny,nz,t22,temp,t22)
      call diff(nx,ny,nz,h,k22,temp,'y',m,d,nl)
      call mm_mult(nx,ny,nz,beta2,temp,temp)
      call mm_add(nx,ny,nz,t22,temp,t22)
      call diff(nx,ny,nz,h,k22,temp,'z',m,d,nl)
      call mm_mult(nx,ny,nz,beta3,temp,temp)
      call mm_add(nx,ny,nz,t22,temp,t22)
      call diff(nx,ny,nz,h,beta2,temp,'x',m,d,nl)
      call mm_mult(nx,ny,nz,k12,temp,temp)
      call mm_sub(nx,ny,nz,t22,temp,t22)
      call diff(nx,ny,nz,h,beta2,temp,'y',m,d,nl)
      call mm_mult(nx,ny,nz,k22,temp,temp)
      call mm_sub(nx,ny,nz,t22,temp,t22)
      call diff(nx,ny,nz,h,beta2,temp,'z',m,d,nl)
      call mm_mult(nx,ny,nz,k32,temp,temp)
      call mm_sub(nx,ny,nz,t22,temp,t22)
      call diff(nx,ny,nz,h,beta1,temp,'y',m,d,nl)
      call mm_mult(nx,ny,nz,k21,temp,temp)
      call mm_add(nx,ny,nz,t22,temp,t22)
      call diff(nx,ny,nz,h,beta2,temp,'y',m,d,nl)
      call mm_mult(nx,ny,nz,k22,temp,temp)
      call mm_add(nx,ny,nz,t22,temp,t22)
      call diff(nx,ny,nz,h,beta3,temp,'y',m,d,nl)
      call mm_mult(nx,ny,nz,k23,temp,temp)
      call mm_add(nx,ny,nz,t22,temp,t22)

*     Lie_beta(K^2_3)
      call diff(nx,ny,nz,h,k23,temp,'x',m,d,nl)
      call mm_mult(nx,ny,nz,beta1,temp,temp)
      call mm_add(nx,ny,nz,t23,temp,t23)
      call diff(nx,ny,nz,h,k23,temp,'y',m,d,nl)
      call mm_mult(nx,ny,nz,beta2,temp,temp)
      call mm_add(nx,ny,nz,t23,temp,t23)
      call diff(nx,ny,nz,h,k23,temp,'z',m,d,nl)
      call mm_mult(nx,ny,nz,beta3,temp,temp)
      call mm_add(nx,ny,nz,t23,temp,t23)
      call diff(nx,ny,nz,h,beta2,temp,'x',m,d,nl)
      call mm_mult(nx,ny,nz,k13,temp,temp)
      call mm_sub(nx,ny,nz,t23,temp,t23)
      call diff(nx,ny,nz,h,beta2,temp,'y',m,d,nl)
      call mm_mult(nx,ny,nz,k23,temp,temp)
      call mm_sub(nx,ny,nz,t23,temp,t23)
      call diff(nx,ny,nz,h,beta2,temp,'z',m,d,nl)
      call mm_mult(nx,ny,nz,k33,temp,temp)
      call mm_sub(nx,ny,nz,t23,temp,t23)
      call diff(nx,ny,nz,h,beta1,temp,'z',m,d,nl)
      call mm_mult(nx,ny,nz,k21,temp,temp)
      call mm_add(nx,ny,nz,t23,temp,t23)
      call diff(nx,ny,nz,h,beta2,temp,'z',m,d,nl)
      call mm_mult(nx,ny,nz,k22,temp,temp)
      call mm_add(nx,ny,nz,t23,temp,t23)
      call diff(nx,ny,nz,h,beta3,temp,'z',m,d,nl)
      call mm_mult(nx,ny,nz,k23,temp,temp)
      call mm_add(nx,ny,nz,t23,temp,t23)

*     Lie_beta(K^3_1)
      call diff(nx,ny,nz,h,k31,temp,'x',m,d,nl)
      call mm_mult(nx,ny,nz,beta1,temp,temp)
      call mm_add(nx,ny,nz,t31,temp,t31)
      call diff(nx,ny,nz,h,k31,temp,'y',m,d,nl)
      call mm_mult(nx,ny,nz,beta2,temp,temp)
      call mm_add(nx,ny,nz,t31,temp,t31)
      call diff(nx,ny,nz,h,k31,temp,'z',m,d,nl)
      call mm_mult(nx,ny,nz,beta3,temp,temp)
      call mm_add(nx,ny,nz,t31,temp,t31)
      call diff(nx,ny,nz,h,beta3,temp,'x',m,d,nl)
      call mm_mult(nx,ny,nz,k11,temp,temp)
      call mm_sub(nx,ny,nz,t31,temp,t31)
      call diff(nx,ny,nz,h,beta3,temp,'y',m,d,nl)
      call mm_mult(nx,ny,nz,k21,temp,temp)
      call mm_sub(nx,ny,nz,t31,temp,t31)
      call diff(nx,ny,nz,h,beta3,temp,'z',m,d,nl)
      call mm_mult(nx,ny,nz,k31,temp,temp)
      call mm_sub(nx,ny,nz,t31,temp,t31)
      call diff(nx,ny,nz,h,beta1,temp,'x',m,d,nl)
      call mm_mult(nx,ny,nz,k31,temp,temp)
      call mm_add(nx,ny,nz,t31,temp,t31)
      call diff(nx,ny,nz,h,beta2,temp,'x',m,d,nl)
      call mm_mult(nx,ny,nz,k32,temp,temp)
      call mm_add(nx,ny,nz,t31,temp,t31)
      call diff(nx,ny,nz,h,beta3,temp,'x',m,d,nl)
      call mm_mult(nx,ny,nz,k33,temp,temp)
      call mm_add(nx,ny,nz,t31,temp,t31)

*     Lie_beta(K^3_2)
      call diff(nx,ny,nz,h,k32,temp,'x',m,d,nl)
      call mm_mult(nx,ny,nz,beta1,temp,temp)
      call mm_add(nx,ny,nz,t32,temp,t32)
      call diff(nx,ny,nz,h,k32,temp,'y',m,d,nl)
      call mm_mult(nx,ny,nz,beta2,temp,temp)
      call mm_add(nx,ny,nz,t32,temp,t32)
      call diff(nx,ny,nz,h,k32,temp,'z',m,d,nl)
      call mm_mult(nx,ny,nz,beta3,temp,temp)
      call mm_add(nx,ny,nz,t32,temp,t32)
      call diff(nx,ny,nz,h,beta3,temp,'x',m,d,nl)
      call mm_mult(nx,ny,nz,k12,temp,temp)
      call mm_sub(nx,ny,nz,t32,temp,t32)
      call diff(nx,ny,nz,h,beta3,temp,'y',m,d,nl)
      call mm_mult(nx,ny,nz,k22,temp,temp)
      call mm_sub(nx,ny,nz,t32,temp,t32)
      call diff(nx,ny,nz,h,beta3,temp,'z',m,d,nl)
      call mm_mult(nx,ny,nz,k32,temp,temp)
      call mm_sub(nx,ny,nz,t32,temp,t32)
      call diff(nx,ny,nz,h,beta1,temp,'y',m,d,nl)
      call mm_mult(nx,ny,nz,k31,temp,temp)
      call mm_add(nx,ny,nz,t32,temp,t32)
      call diff(nx,ny,nz,h,beta2,temp,'y',m,d,nl)
      call mm_mult(nx,ny,nz,k32,temp,temp)
      call mm_add(nx,ny,nz,t32,temp,t32)
      call diff(nx,ny,nz,h,beta3,temp,'y',m,d,nl)
      call mm_mult(nx,ny,nz,k33,temp,temp)
      call mm_add(nx,ny,nz,t32,temp,t32)

*     Lie_beta(K^3_3)
      call diff(nx,ny,nz,h,k33,temp,'x',m,d,nl)
      call mm_mult(nx,ny,nz,beta1,temp,temp)
      call mm_add(nx,ny,nz,t33,temp,t33)
      call diff(nx,ny,nz,h,k33,temp,'y',m,d,nl)
      call mm_mult(nx,ny,nz,beta2,temp,temp)
      call mm_add(nx,ny,nz,t33,temp,t33)
      call diff(nx,ny,nz,h,k33,temp,'z',m,d,nl)
      call mm_mult(nx,ny,nz,beta3,temp,temp)
      call mm_add(nx,ny,nz,t33,temp,t33)
      call diff(nx,ny,nz,h,beta3,temp,'x',m,d,nl)
      call mm_mult(nx,ny,nz,k13,temp,temp)
      call mm_sub(nx,ny,nz,t33,temp,t33)
      call diff(nx,ny,nz,h,beta3,temp,'y',m,d,nl)
      call mm_mult(nx,ny,nz,k23,temp,temp)
      call mm_sub(nx,ny,nz,t33,temp,t33)
      call diff(nx,ny,nz,h,beta3,temp,'z',m,d,nl)
      call mm_mult(nx,ny,nz,k33,temp,temp)
      call mm_sub(nx,ny,nz,t33,temp,t33)
      call diff(nx,ny,nz,h,beta1,temp,'z',m,d,nl)
      call mm_mult(nx,ny,nz,k31,temp,temp)
      call mm_add(nx,ny,nz,t33,temp,t33)
      call diff(nx,ny,nz,h,beta2,temp,'z',m,d,nl)
      call mm_mult(nx,ny,nz,k32,temp,temp)
      call mm_add(nx,ny,nz,t33,temp,t33)
      call diff(nx,ny,nz,h,beta3,temp,'z',m,d,nl)
      call mm_mult(nx,ny,nz,k33,temp,temp)
      call mm_add(nx,ny,nz,t33,temp,t33)

      return
      end
      
************************************************************************
*     Computes alpha*(R^a_b + K*K^a_b) terms
************************************************************************
      subroutine evolve_k3(gup11,gup12,gup13,gup22,gup23,gup33,
     $                     r11,r12,r13,r22,r23,r33,
     $                     k11,k12,k13,
     $                     k21,k22,k23,
     $                     k31,k32,k33,
     $                     alpha,
     $                     t11,t12,t13,
     $                     t21,t22,t23,
     $                     t31,t32,t33,temp,
     $                     nx,ny,nz)

      implicit none

      integer  nx, ny, nz
      integer  i, j, k

      real*8   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   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   k11(nx,ny,nz), k12(nx,ny,nz), k13(nx,ny,nz),
     $         k21(nx,ny,nz), k22(nx,ny,nz), k23(nx,ny,nz),
     $         k31(nx,ny,nz), k32(nx,ny,nz), k33(nx,ny,nz)

      real*8   temp(nx,ny,nz), alpha(nx,ny,nz)

      real*8   t11(nx,ny,nz), t12(nx,ny,nz), t13(nx,ny,nz),
     $         t21(nx,ny,nz), t22(nx,ny,nz), t23(nx,ny,nz),
     $         t31(nx,ny,nz), t32(nx,ny,nz), t33(nx,ny,nz)

      call copy(nx,ny,nz,k11,temp)
      call mm_add(nx,ny,nz,k22,temp,temp)
      call mm_add(nx,ny,nz,k33,temp,temp)

c     alpha*(R^a_b + K*K^a_b)
c     do i = 1, nx
c     do j = 1, ny
c     do k = 1, nz

c        tab(i,j,k) = tab(i,j,k) + alpha(i,j,k)*(
c    $                (gupa1(i,j,k)*r1b(i,j,k)
c    $                +gupa2(i,j,k)*r2b(i,j,k)
c    $                +gupa3(i,j,k)*r3b(i,j,k))
c    $                +temp(i,j,k)*kab(i,j,k)   )

c     end do
c     end do
c     end do

c     alpha*(R^1_1 + K*K^1_1)
      do i = 1, nx
      do j = 1, ny
      do k = 1, nz

         t11(i,j,k) = t11(i,j,k) + alpha(i,j,k)*(
     $                (gup11(i,j,k)*r11(i,j,k)
     $                +gup12(i,j,k)*r12(i,j,k)
     $                +gup13(i,j,k)*r13(i,j,k))
     $                +temp(i,j,k)*k11(i,j,k)   )

      end do
      end do
      end do

c     alpha*(R^1_2 + K*K^1_2)
      do i = 1, nx
      do j = 1, ny
      do k = 1, nz

         t12(i,j,k) = t12(i,j,k) + alpha(i,j,k)*(
     $                (gup11(i,j,k)*r12(i,j,k)
     $                +gup12(i,j,k)*r22(i,j,k)
     $                +gup13(i,j,k)*r23(i,j,k))
     $                +temp(i,j,k)*k12(i,j,k)   )

      end do
      end do
      end do

c     alpha*(R^1_3 + K*K^1_3)
      do i = 1, nx
      do j = 1, ny
      do k = 1, nz

         t13(i,j,k) = t13(i,j,k) + alpha(i,j,k)*(
     $                (gup11(i,j,k)*r13(i,j,k)
     $                +gup12(i,j,k)*r23(i,j,k)
     $                +gup13(i,j,k)*r33(i,j,k))
     $                +temp(i,j,k)*k13(i,j,k)   )

      end do
      end do
      end do

c     alpha*(R^2_1 + K*K^2_1)
      do i = 1, nx
      do j = 1, ny
      do k = 1, nz

         t21(i,j,k) = t21(i,j,k) + alpha(i,j,k)*(
     $                (gup12(i,j,k)*r11(i,j,k)
     $                +gup22(i,j,k)*r12(i,j,k)
     $                +gup23(i,j,k)*r13(i,j,k))
     $                +temp(i,j,k)*k21(i,j,k)   )

      end do
      end do
      end do

c     alpha*(R^2_2 + K*K^2_2)
      do i = 1, nx
      do j = 1, ny
      do k = 1, nz

         t22(i,j,k) = t22(i,j,k) + alpha(i,j,k)*(
     $                (gup12(i,j,k)*r12(i,j,k)
     $                +gup22(i,j,k)*r22(i,j,k)
     $                +gup23(i,j,k)*r23(i,j,k))
     $                +temp(i,j,k)*k22(i,j,k)   )

      end do
      end do
      end do

c     alpha*(R^2_3 + K*K^2_3)
      do i = 1, nx
      do j = 1, ny
      do k = 1, nz

         t23(i,j,k) = t23(i,j,k) + alpha(i,j,k)*(
     $                (gup12(i,j,k)*r13(i,j,k)
     $                +gup22(i,j,k)*r23(i,j,k)
     $                +gup23(i,j,k)*r33(i,j,k))
     $                +temp(i,j,k)*k23(i,j,k)   )

      end do
      end do
      end do

c     alpha*(R^3_1 + K*K^3_1)
      do i = 1, nx
      do j = 1, ny
      do k = 1, nz

         t31(i,j,k) = t31(i,j,k) + alpha(i,j,k)*(
     $                (gup13(i,j,k)*r11(i,j,k)
     $                +gup23(i,j,k)*r12(i,j,k)
     $                +gup33(i,j,k)*r13(i,j,k))
     $                +temp(i,j,k)*k31(i,j,k)   )

      end do
      end do
      end do

c     alpha*(R^3_2 + K*K^3_2)
      do i = 1, nx
     		do j = 1, ny
      		do k = 1, nz

         		t32(i,j,k) = t32(i,j,k) + alpha(i,j,k)*(
     $                      (gup13(i,j,k)*r12(i,j,k)
     $                      +gup23(i,j,k)*r22(i,j,k)
     $                      +gup33(i,j,k)*r23(i,j,k))
     $                      +temp(i,j,k)*k32(i,j,k)   )

      		end do
      	end do
      end do

c     alpha*(R^3_3 + K*K^3_3)
      do i = 1, nx
      	do j = 1, ny
      		do k = 1, nz

         		t33(i,j,k) = t33(i,j,k) + alpha(i,j,k)*(
     $                   (gup13(i,j,k)*r13(i,j,k)
     $                   +gup23(i,j,k)*r23(i,j,k)
     $                   +gup33(i,j,k)*r33(i,j,k))
     $                   +temp(i,j,k)*k33(i,j,k)   )

      		end do
      	end do
      end do

      return
      end
