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

c-----------------------------------------------------------------------
C This routine will be in the C++ routine in DAGH
      subroutine evolve_g(g11, g12, g13, g22, g23, g33,
     *           gup11, gup12, gup13, gup22, gup23, gup33,
     *           g11m1,g12m1,g13m1,g22m1,g23m1,g33m1,
     *           g11p1,g12p1,g13p1,g22p1,g23p1,g33p1,
     *           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

      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   g11m1(nx,ny,nz), g12m1(nx,ny,nz), g13m1(nx,ny,nz), 
     *           g22m1(nx,ny,nz), g23m1(nx,ny,nz), g33m1(nx,ny,nz)

      real*8   g11p1(nx,ny,nz), g12p1(nx,ny,nz), g13p1(nx,ny,nz), 
     *           g22p1(nx,ny,nz), g23p1(nx,ny,nz), g33p1(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   Computes the derivative of beta's 
      call evolve_g1(nx,ny,nz,h,beta1,beta2,beta3,
     .               t11,t12,t13,t21,t22,t23,t31,t32,t33,temp1,m,d,nl)


C     Sync(t11,t,l)
C     Sync(t12,t,l)
C     Sync(t13,t,l)
C     Sync(t22,t,l)
C     Sync(t23,t,l)
C     Sync(t31,t,l)
C     Sync(t32,t,l)
C     Sync(t33,t,l)

      call evolve_g2(nx,ny,nz,beta1,beta2,beta3,
     *           g11, g12, g13, g22, g23, g33,
     *           c111, c112, c113, c122, c123, c133,
     *           c211, c212, c213, c222, c223, c233,
     *           c311, c312, c313, c322, c323, c333,
     *           t11,t12,t13,t21,t22,t23,t31,t32,t33,
     *           temp1, temp2, temp3 ,m,d,nl) 


      call evolve_g34(nx,ny,nz,beta1,beta2,beta3,
     *                t11,t12,t13,t21,t22,t23,t31,t32,t33,temp1,m,d,nl)


      call evolve_g5(nx,ny,nz,g11,g12,g13,g22,g23,g33,
     *              k11,k12,k13,k21,k22,k23,k31,k32,k33,alpha,
     *              t11,t12,t13,t21,t22,t23,t31,t32,t33,temp1,m,d,nl)


      call evolve_g6(nx,ny,nz,dt,
     *               g11,g12,g13,g22,g23,g33,
     *               g11m1,g12m1,g13m1,g22m1,g23m1,g33m1,
     *               g11p1,g12p1,g13p1,g22p1,g23p1,g33p1,
     *               t11,t12,t13,t21,t22,t23,t31,t32,t33,temp1,m,d,nl)


      return 
      end
c-----------------------------------------------------------------------
c	routine checked to be 2nd order in h
c
      subroutine evolve_g1(nx,ny,nz,h,beta1,beta2,beta3,
     .                     beta11,beta12,beta13,
     .                     beta21,beta22,beta23,
     .                     beta31,beta32,beta33,temp,m,d,nl)
      implicit none
      integer nx,ny,nz,nl
      real*8  h
      real*8 beta1(nx,ny,nz),beta2(nx,ny,nz),beta3(nx,ny,nz)
      real*8 beta11(nx,ny,nz),beta12(nx,ny,nz),beta13(nx,ny,nz)
      real*8 beta21(nx,ny,nz),beta22(nx,ny,nz),beta23(nx,ny,nz)
      real*8 beta31(nx,ny,nz),beta32(nx,ny,nz),beta33(nx,ny,nz)
      real*8 temp(nx,ny,nz)
      integer m(3*nl)
      real*8  d(3*nl)

      call diff(nx,ny,nz,h,beta1,beta11,'x',m,d,nl)
      call diff(nx,ny,nz,h,beta1,beta12,'y',m,d,nl)
      call diff(nx,ny,nz,h,beta1,beta13,'z',m,d,nl)

      call diff(nx,ny,nz,h,beta2,beta21,'x',m,d,nl)
      call diff(nx,ny,nz,h,beta2,beta22,'y',m,d,nl)
      call diff(nx,ny,nz,h,beta2,beta23,'z',m,d,nl)

      call diff(nx,ny,nz,h,beta3,beta31,'x',m,d,nl)
      call diff(nx,ny,nz,h,beta3,beta32,'y',m,d,nl)
      call diff(nx,ny,nz,h,beta3,beta33,'z',m,d,nl)

      return
      end

c-----------------------------------------------------------------------
c found to be O(h^2) with convergence factor ~3.9
c
      subroutine evolve_g2(nx,ny,nz,beta1,beta2,beta3,
     .                     g11, g12, g13, g22, g23, g33,
     *           c111, c112, c113, c122, c123, c133,
     *           c211, c212, c213, c222, c223, c233,
     *           c311, c312, c313, c322, c323, c333,
     *           t11,t12,t13,t21,t22,t23,t31,t32,t33,
     *           temp1, temp2, temp3 ,m,d,nl) 
      implicit none
      integer nx,ny,nz,nl
      real*8 beta1(nx,ny,nz),beta2(nx,ny,nz),beta3(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)
      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)
      real*8   t21(nx,ny,nz),t22(nx,ny,nz),t23(nx,ny,nz)
      real*8   t31(nx,ny,nz),t32(nx,ny,nz),t33(nx,ny,nz)
      real*8   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 mm_mult will be matrix matrix multiplication t1 = c11 * beta1 last arg is 
C the storage aguement
C mm_add will be to add two matrixs   t1 = a + b
C  Tki = tki + cki1*b1 + cki2*b2 + cki3*b3
C T11:
      call mm_mult(nx,ny,nz,c111,beta1,temp1)
      call mm_add(nx,ny,nz,t11,temp1,t11)
      call mm_mult(nx,ny,nz,c112,beta2,temp1)
      call mm_add(nx,ny,nz,t11,temp1,t11)
      call mm_mult(nx,ny,nz,c113,beta3,temp1)
      call mm_add(nx,ny,nz,t11,temp1,t11)
C T12:
      call mm_mult(nx,ny,nz,c112,beta1,temp1)
      call mm_add(nx,ny,nz,t12,temp1,t12)
      call mm_mult(nx,ny,nz,c122,beta2,temp1)
      call mm_add(nx,ny,nz,t12,temp1,t12)
      call mm_mult(nx,ny,nz,c123,beta3,temp1)
      call mm_add(nx,ny,nz,t12,temp1,t12)
C T13:
      call mm_mult(nx,ny,nz,c113,beta1,temp1)
      call mm_add(nx,ny,nz,t13,temp1,t13)
      call mm_mult(nx,ny,nz,c123,beta2,temp1)
      call mm_add(nx,ny,nz,t13,temp1,t13)
      call mm_mult(nx,ny,nz,c133,beta3,temp1)
      call mm_add(nx,ny,nz,t13,temp1,t13)
C T21:
      call mm_mult(nx,ny,nz,c211,beta1,temp1)
      call mm_add(nx,ny,nz,t21,temp1,t21)
      call mm_mult(nx,ny,nz,c212,beta2,temp1)
      call mm_add(nx,ny,nz,t21,temp1,t21)
      call mm_mult(nx,ny,nz,c213,beta3,temp1)
      call mm_add(nx,ny,nz,t21,temp1,t21)
C T22:
      call mm_mult(nx,ny,nz,c212,beta1,temp1)
      call mm_add(nx,ny,nz,t22,temp1,t22)
      call mm_mult(nx,ny,nz,c222,beta2,temp1)
      call mm_add(nx,ny,nz,t22,temp1,t22)
      call mm_mult(nx,ny,nz,c223,beta3,temp1)
      call mm_add(nx,ny,nz,t22,temp1,t22)
C T23:
      call mm_mult(nx,ny,nz,c213,beta1,temp1)
      call mm_add(nx,ny,nz,t23,temp1,t23)
      call mm_mult(nx,ny,nz,c223,beta2,temp1)
      call mm_add(nx,ny,nz,t23,temp1,t23)
      call mm_mult(nx,ny,nz,c233,beta3,temp1)
      call mm_add(nx,ny,nz,t23,temp1,t23)
C T31:
      call mm_mult(nx,ny,nz,c311,beta1,temp1)
      call mm_add(nx,ny,nz,t31,temp1,t31)
      call mm_mult(nx,ny,nz,c312,beta2,temp1)
      call mm_add(nx,ny,nz,t31,temp1,t31)
      call mm_mult(nx,ny,nz,c313,beta3,temp1)
      call mm_add(nx,ny,nz,t31,temp1,t31)
C T32:
      call mm_mult(nx,ny,nz,c312,beta1,temp1)
      call mm_add(nx,ny,nz,t32,temp1,t32)
      call mm_mult(nx,ny,nz,c322,beta2,temp1)
      call mm_add(nx,ny,nz,t32,temp1,t32)
      call mm_mult(nx,ny,nz,c323,beta3,temp1)
      call mm_add(nx,ny,nz,t32,temp1,t32)
C T33:                                                      
      call mm_mult(nx,ny,nz,c313,beta1,temp1)
      call mm_add(nx,ny,nz,t33,temp1,t33)
      call mm_mult(nx,ny,nz,c323,beta2,temp1)
      call mm_add(nx,ny,nz,t33,temp1,t33)
      call mm_mult(nx,ny,nz,c333,beta3,temp1)
      call mm_add(nx,ny,nz,t33,temp1,t33)

c Lower second index of Tij

c    tij = tik * gkj => store tik = temp1k for every i,j

c  1k
	   call copy(nx, ny, nz, t11 , temp1)
	   call copy(nx, ny, nz, t12 , temp2)
	   call copy(nx, ny, nz, t13 , temp3)

      do k = 1, nz
        do j = 1, ny
          do i = 1, nx
            t11(i,j,k) = temp1(i,j,k) * g11(i,j,k) +
     .                   temp2(i,j,k) * g12(i,j,k) + 
     .                   temp3(i,j,k) * g13(i,j,k)
          end do
        end do
      end do

      do k = 1, nz
        do j = 1, ny
          do i = 1, nx
            t12(i,j,k) = temp1(i,j,k) * g12(i,j,k) +
     .                   temp2(i,j,k) * g22(i,j,k) + 
     .                   temp3(i,j,k) * g23(i,j,k)
          end do
        end do
      end do

      do k = 1, nz
        do j = 1, ny
          do i = 1, nx
            t13(i,j,k) = temp1(i,j,k) * g13(i,j,k) +
     .                   temp2(i,j,k) * g23(i,j,k) + 
     .                   temp3(i,j,k) * g33(i,j,k)
          end do
        end do
      end do

c 2k
	   call copy(nx, ny, nz, t21 , temp1)
	   call copy(nx, ny, nz, t22 , temp2)
	   call copy(nx, ny, nz, t23 , temp3)

      do k = 1, nz
        do j = 1, ny
          do i = 1, nx
            t21(i,j,k) = temp1(i,j,k) * g11(i,j,k) +
     .                   temp2(i,j,k) * g12(i,j,k) + 
     .                   temp3(i,j,k) * g13(i,j,k)
          end do
        end do
      end do

      do k = 1, nz
        do j = 1, ny
          do i = 1, nx
            t22(i,j,k) = temp1(i,j,k) * g12(i,j,k) +
     .                   temp2(i,j,k) * g22(i,j,k) + 
     .                   temp3(i,j,k) * g23(i,j,k)
          end do
        end do
      end do

      do k = 1, nz
        do j = 1, ny
          do i = 1, nx
            t23(i,j,k) = temp1(i,j,k) * g13(i,j,k) +
     .                   temp2(i,j,k) * g23(i,j,k) + 
     .                   temp3(i,j,k) * g33(i,j,k)
          end do
        end do
      end do

c 3k
	   call copy(nx, ny, nz, t31 , temp1)
	   call copy(nx, ny, nz, t32 , temp2)
	   call copy(nx, ny, nz, t33 , temp3)

      do k = 1, nz
        do j = 1, ny
          do i = 1, nx
            t31(i,j,k) = temp1(i,j,k) * g11(i,j,k) +
     .                   temp2(i,j,k) * g12(i,j,k) + 
     .                   temp3(i,j,k) * g13(i,j,k)
          end do
        end do
      end do

      do k = 1, nz
        do j = 1, ny
          do i = 1, nx
            t32(i,j,k) = temp1(i,j,k) * g12(i,j,k) +
     .                   temp2(i,j,k) * g22(i,j,k) + 
     .                   temp3(i,j,k) * g23(i,j,k)
          end do
        end do
      end do

      do k = 1, nz
        do j = 1, ny
          do i = 1, nx
            t33(i,j,k) = temp1(i,j,k) * g13(i,j,k) +
     .                   temp2(i,j,k) * g23(i,j,k) + 
     .                   temp3(i,j,k) * g33(i,j,k)
          end do
        end do
      end do

      return
      end
C-----------------------------------------------------------------------------
      subroutine evolve_g34(nx,ny,nz,beta1,beta2,beta3,
     *                t11,t12,t13,t21,t22,t23,t31,t32,t33,temp,m,d,nl)
      implicit none
      integer nx,ny,nz,nl
      real*8 beta1(nx,ny,nz)
      real*8 beta2(nx,ny,nz)
      real*8 beta3(nx,ny,nz)


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

      integer m(3*nl)
      real*8  d(3*nl)
C T11:
      call mm_add(nx,ny,nz,t11,t11,t11)

C T12: T12 = t12 + t21
      call copy(nx,ny,nz,t21,temp)
      call mm_add(nx,ny,nz,t12,temp,t12)

C T13: T13 = t13 + t31
      call copy(nx,ny,nz,t31,temp)
      call mm_add(nx,ny,nz,t13,temp,t13)

C T21: T21 = T12 
      call copy(nx,ny,nz,t12,t21)

C T22:
      call mm_add(nx,ny,nz,t22,t22,t22)

C T23: T23 = t23 + t32
      call copy(nx,ny,nz,t32,temp)
      call mm_add(nx,ny,nz,t23,temp,t23)

C T31: T31 = T13
      call copy(nx,ny,nz,t13,t31)

C T32: T32 = T23
      call copy(nx,ny,nz,t23,t32)

C T33:
      call mm_add(nx,ny,nz,t33,t33,t33)

      return
      end
C---------------------------------------------------------------------------
      subroutine evolve_g5(nx,ny,nz,g11,g12,g13,g22,g23,g33,
     *              k11,k12,k13,k21,k22,k23,k31,k32,k33,alpha,
     *              t11,t12,t13,t21,t22,t23,t31,t32,t33,temp,m,d,nl)
		implicit none
      integer nx,ny,nz,nl

      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   g11(nx,ny,nz), g12(nx,ny,nz), g13(nx,ny,nz)
      real*8   g22(nx,ny,nz), g23(nx,ny,nz), g33(nx,ny,nz)

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

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

      integer i,j,k


C    tij = tij - 2*alpha*(g1i*K1j+g2i*K2j+g3i*K3j)
C but allow for g21 = g12 , g31 = g13, g32 = g23

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

               t11(i,j,k) = t11(i,j,k) - 2.0e0*alpha(i,j,k)* (
     .                      g11(i,j,k)*k11(i,j,k) +
     .                      g12(i,j,k)*k21(i,j,k) +
     .                      g13(i,j,k)*k31(i,j,k) ) 

               t12(i,j,k) = t12(i,j,k) - 2.0e0*alpha(i,j,k)* (
     .                      g11(i,j,k)*k12(i,j,k) +
     .                      g12(i,j,k)*k22(i,j,k) +
     .                      g13(i,j,k)*k32(i,j,k) )

               t13(i,j,k) = t13(i,j,k) - 2.0e0*alpha(i,j,k)* (
     .                      g11(i,j,k)*k13(i,j,k) +
     .                      g12(i,j,k)*k23(i,j,k) +
     .                      g13(i,j,k)*k33(i,j,k) )

               t21(i,j,k) = t21(i,j,k) - 2.0e0*alpha(i,j,k)* (
     .                      g12(i,j,k)*k11(i,j,k) +
     .                      g22(i,j,k)*k21(i,j,k) +
     .                      g23(i,j,k)*k31(i,j,k) )

               t22(i,j,k) = t22(i,j,k) - 2.0e0*alpha(i,j,k)* (
     .                      g12(i,j,k)*k12(i,j,k) +
     .                      g22(i,j,k)*k22(i,j,k) +
     .                      g23(i,j,k)*k32(i,j,k) )

               t23(i,j,k) = t23(i,j,k) - 2.0e0*alpha(i,j,k)* (
     .                      g12(i,j,k)*k13(i,j,k) +
     .                      g22(i,j,k)*k23(i,j,k) +
     .                      g23(i,j,k)*k33(i,j,k) )

               t31(i,j,k) = t31(i,j,k) - 2.0e0*alpha(i,j,k)* (
     .                      g13(i,j,k)*k11(i,j,k) +
     .                      g23(i,j,k)*k21(i,j,k) +
     .                      g33(i,j,k)*k31(i,j,k) )

               t32(i,j,k) = t32(i,j,k) - 2.0e0*alpha(i,j,k)* (
     .                      g13(i,j,k)*k12(i,j,k) +
     .                      g23(i,j,k)*k22(i,j,k) +
     .                      g33(i,j,k)*k32(i,j,k) )

               t33(i,j,k) = t33(i,j,k) - 2.0e0*alpha(i,j,k)* (
     .                      g13(i,j,k)*k13(i,j,k) +
     .                      g23(i,j,k)*k23(i,j,k) +
     .                      g33(i,j,k)*k33(i,j,k) )

            end do
         end do
      end do

      return
      end
C----------------------------------------------------------------------------
      subroutine evolve_g6(nx,ny,nz,dt,
     *               g11,g12,g13,g22,g23,g33,
     *               g11m1,g12m1,g13m1,g22m1,g23m1,g33m1,
     *               g11p1,g12p1,g13p1,g22p1,g23p1,g33p1,
     *               t11,t12,t13,t21,t22,t23,t31,t32,t33,temp,m,d,nl)
		implicit none
		integer nx,ny,nz,nl

		real*8 g11(nx,ny,nz),g12(nx,ny,nz),g13(nx,ny,nz)
		real*8 g22(nx,ny,nz),g23(nx,ny,nz)
		real*8 g33(nx,ny,nz)

		real*8 g11m1(nx,ny,nz),g12m1(nx,ny,nz),g13m1(nx,ny,nz)
		real*8 g22m1(nx,ny,nz),g23m1(nx,ny,nz)
		real*8 g33m1(nx,ny,nz)

		real*8 g11p1(nx,ny,nz),g12p1(nx,ny,nz),g13p1(nx,ny,nz)
		real*8 g22p1(nx,ny,nz),g23p1(nx,ny,nz)
		real*8 g33p1(nx,ny,nz)

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

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

		real*8 dt

		integer i,j,k

		call leapfrog(nx,ny,nz,g11m1,g11p1,t11,dt) 
		call leapfrog(nx,ny,nz,g12m1,g12p1,t12,dt) 
		call leapfrog(nx,ny,nz,g13m1,g13p1,t13,dt) 

      call leapfrog(nx,ny,nz,g22m1,g22p1,t22,dt)
      call leapfrog(nx,ny,nz,g23m1,g23p1,t23,dt)

      call leapfrog(nx,ny,nz,g33m1,g33p1,t33,dt)

		return
		end
C----------------------------------------------------------------------------
		subroutine leapfrog(nx,ny,nz,gm1,gp1,lg,dt)
		implicit none
		integer nx,ny,nz
		real*8 gm1(nx,ny,nz),gp1(nx,ny,nz),lg(nx,ny,nz)
		real*8 dt

		integer i,j,k

		do k = 1 , nz
			do j = 1 , ny
				do i = 1 , nx
					gp1(i,j,k) = 2.0e0*dt*lg(i,j,k) + gm1(i,j,k)
				end do
			end do
		end do

		return
		end