      subroutine t2_initdata(lb,ub,shape,g11,g12,
     &                       g13,g22,g23,g33,
     &                       k11,k12,k13,k21,k22,k23,k31,k32,k33,
     &                       bha,lbr,ubr,shaper,bhr,temp1,
     &                       bb,dx,dy,dz)
      implicit none
      integer lb(3),ub(3),shape(3)
      real*8   g11(shape(1),shape(2),shape(3)),
     &         g12(shape(1),shape(2),shape(3)),
     &         g13(shape(1),shape(2),shape(3)),
     &         g22(shape(1),shape(2),shape(3)),
     &         g23(shape(1),shape(2),shape(3)),
     &         g33(shape(1),shape(2),shape(3))
c     extrinsic curvature:  K^i_j
      real*8   k11(shape(1),shape(2),shape(3)),
     &         k12(shape(1),shape(2),shape(3)),
     &         k13(shape(1),shape(2),shape(3)),
     &         k21(shape(1),shape(2),shape(3)),
     &         k22(shape(1),shape(2),shape(3)),
     &         k23(shape(1),shape(2),shape(3)),
     &         k31(shape(1),shape(2),shape(3)),
     &         k32(shape(1),shape(2),shape(3)),
     &         k33(shape(1),shape(2),shape(3))
      real*8   temp1(shape(1),shape(2),shape(3))
      real*8   dx,dy,dz
      real*8   bb(6)
      integer lbr(2),ubr(2),shaper(2)

      integer maxn
      parameter ( maxn = 2049 )
      real*8    x(maxn), y(maxn),z(maxn),bha(2),bhr(shaper(1),shaper(2))
      integer nx,ny,nz

      nx = shape(1)
      ny = shape(2)
      nz = shape(3)

      if (nx.gt.maxn.or.ny.gt.maxn.or.nz.gt.maxn) then
         write(0,*) '******** increase maxn in all routines ***********'
         stop
      end if

      call makegrid(nx,ny,nz,bb(1),bb(2),bb(3),bb(4),bb(5),bb(6),
     &              dx,dy,dz,x,y,z,lb)

      call initdata(nx,ny,nz,g11,g12,g13,g22,g23,g33,
     $     k11,k12,k13,k21,k22,k23,k31,k32,k33,
     $     x,y,z,bha,bhr,temp1)

      return
      end
C-----------------------------------------------------------------------------
      subroutine t2_horizon(lb,ub,shape,g11,g12,g13,g22,g23,g33,
     &           k11,k12,k13,k21,k22,k23,k31,k32,k33,
     &           bha,lbr,ubr,shaper,bhr,mask,bb,dx,dy,dz)

      implicit none
      integer lb(3),ub(3),shape(3)
      integer lbr(2),ubr(2),shaper(2)

      integer  nx,ny,nz
c     metric:  g_ij
      real*8   g11(shape(1),shape(2),shape(3)),
     &         g12(shape(1),shape(2),shape(3)),
     &         g13(shape(1),shape(2),shape(3)),
     &         g22(shape(1),shape(2),shape(3)),
     &         g23(shape(1),shape(2),shape(3)),
     &         g33(shape(1),shape(2),shape(3))
      real*8   k11(shape(1),shape(2),shape(3)),
     &         k12(shape(1),shape(2),shape(3)),
     &         k13(shape(1),shape(2),shape(3)),
     &         k21(shape(1),shape(2),shape(3)),
     &         k22(shape(1),shape(2),shape(3)),
     &         k23(shape(1),shape(2),shape(3)),
     &         k31(shape(1),shape(2),shape(3)),
     &         k32(shape(1),shape(2),shape(3)),
     &         k33(shape(1),shape(2),shape(3))

      integer maxn
      parameter ( maxn = 2049 )

      real*8   x(maxn), y(maxn), z(maxn), bha(2),
     .         bhr(shaper(1),shaper(2))
      integer  mask(shape(1),shape(2),shape(3))
      real*8   dx,dy,dz
      real*8   bb(6)

      nx = shape(1)
      ny = shape(2)
      nz = shape(3)

      call makegrid(nx,ny,nz,bb(1),bb(2),bb(3),bb(4),bb(5),bb(6),
     &              dx,dy,dz,x,y,z,lb)

      call horizon(nx,ny,nz,g11,g12,g13,g22,g23,g33,
     $     k11,k12,k13,k21,k22,k23,k31,k32,k33,
     $     x,y,z,bha,bhr,mask)

      return
      end
C-----------------------------------------------------------------------------
      subroutine t2_remask(lb,ub,shape,g11,g12,g13,g22,g23,g33,
     &           k11,k12,k13,k21,k22,k23,k31,k32,k33,
     &           mask,maskold)
      implicit none
      integer lb(3),ub(3),shape(3)
      integer nx,ny,nz

      real*8   g11(shape(1),shape(2),shape(3)),
     &         g12(shape(1),shape(2),shape(3)),
     &         g13(shape(1),shape(2),shape(3)),
     &         g22(shape(1),shape(2),shape(3)),
     &         g23(shape(1),shape(2),shape(3)),
     &         g33(shape(1),shape(2),shape(3))
      real*8   k11(shape(1),shape(2),shape(3)),
     &         k12(shape(1),shape(2),shape(3)),
     &         k13(shape(1),shape(2),shape(3)),
     &         k21(shape(1),shape(2),shape(3)),
     &         k22(shape(1),shape(2),shape(3)),
     &         k23(shape(1),shape(2),shape(3)),
     &         k31(shape(1),shape(2),shape(3)),
     &         k32(shape(1),shape(2),shape(3)),
     &         k33(shape(1),shape(2),shape(3))
      real*8   mask(shape(1),shape(2),shape(3))
      real*8   maskold(shape(1),shape(2),shape(3))

      nx = shape(1)
      ny = shape(2)
      nz = shape(3)

      call remask(nx,ny,nz,g11,g12,g13,g22,g23,g33,
     $     k11,k12,k13,k21,k22,k23,k31,k32,k33,
     $     mask,maskold)

      return
      end
C-----------------------------------------------------------------------------
      subroutine t2_flags(lb,ub,shape,nlmax,xi,yj,zk,
     .                    fx,fy,fz,mask,lbm,ubm,shapem,m,d,nl)
      implicit none
      integer lb(3),ub(3),shape(3)
		integer lbm(3),ubm(3),shapem(3)
      real*8 mask(shape(1),shape(2),shape(3))
      integer nlmax,nl
      real*8 d(*)
      integer m(*)
      integer  xi(nlmax),yj(nlmax),zk(nlmax)
      real*8   fx(nlmax),fy(nlmax),fz(nlmax)


      integer nx,ny,nz
		integer l

      nx = shape(1)
      ny = shape(2)
      nz = shape(3)

      call flags(nx,ny,nz,nlmax,xi,yj,zk,fx,fy,fz,mask,m,d,nl)

      return
      end
C-----------------------------------------------------------------------------
      subroutine t2_setgauge(lb,ub,shape,mask,bha,
     .                          lbr,ubr,shaper,bhr,
     .                          alpha,beta1,beta2,beta3,
     .                          bb,dx,dy,dz)
      implicit none
      integer lb(3),ub(3),shape(3)

      real*8 mask(shape(1),shape(2),shape(3))
      real*8 alpha(shape(1),shape(2),shape(3))
      real*8 beta1(shape(1),shape(2),shape(3))
      real*8 beta2(shape(1),shape(2),shape(3))
      real*8 beta3(shape(1),shape(2),shape(3))
      integer maxn
      parameter ( maxn = 2049 )
      real*8  x(maxn), y(maxn),z(maxn),bha(2)
      integer nx,ny,nz
      real*8 dx,dy,dz
      real*8 bb(6)

      integer lbr(2),ubr(2),shaper(2)
      real*8 bhr(shaper(1),shaper(2))

      nx = shape(1)
      ny = shape(2)
      nz = shape(3)

      call makegrid(nx,ny,nz,bb(1),bb(2),bb(3),bb(4),bb(5),bb(6),
     &              dx,dy,dz,x,y,z,lb)

      call setgauge(nx,ny,nz,mask,
     $     x,y,z,bha,bhr,
     $     alpha,beta1,beta2,beta3)

      return
      end
C-----------------------------------------------------------------------------
      subroutine t2_constraint(lb,ub,shape,mask,
     .                         k11,k12,k13,k21,k22,k23,
     .                         k31,k32,k33,
     .                         hamnorm,mom1norm,
     .                         mom2norm,mom3norm,
     .                         c111, c112, c113,
     .                         c122, c123, c133,
     .                         c211, c212, c213,
     .                         c222, c223, c233,
     .                         c311, c312, c313,
     .                         c322, c323, c333,
     .                         r11, r12, r13,
     .                         r22, r23, r33,
     .                         gup11,gup12,gup13,
     .                         gup22, gup23, gup33,
     .                         t11, t12, t13,
     .                         t21, t22, t23,
     .                         t31, t32, t33,
     .                         temp1,dx,nl,lbm,ubm,shapem,m,d)


      implicit none
      integer lb(3),ub(3),shape(3)
		integer lbm(3),ubm(3),shapem(3)
      real*8 dx
      integer nx,ny,nz
      real*8   k11(shape(1),shape(2),shape(3)),
     &         k12(shape(1),shape(2),shape(3)),
     &         k13(shape(1),shape(2),shape(3)),
     &         k21(shape(1),shape(2),shape(3)),
     &         k22(shape(1),shape(2),shape(3)),
     &         k23(shape(1),shape(2),shape(3)),
     &         k31(shape(1),shape(2),shape(3)),
     &         k32(shape(1),shape(2),shape(3)),
     &         k33(shape(1),shape(2),shape(3))
		real*8   mask(shape(1),shape(2),shape(3))
		real*8 hamnorm,mom1norm,mom2norm,mom3norm
      real*8 c111(shape(1),shape(2),shape(3))
      real*8 c112(shape(1),shape(2),shape(3))
      real*8 c113(shape(1),shape(2),shape(3))
      real*8 c122(shape(1),shape(2),shape(3))
      real*8 c123(shape(1),shape(2),shape(3))
      real*8 c133(shape(1),shape(2),shape(3))
      real*8 c211(shape(1),shape(2),shape(3))
      real*8 c212(shape(1),shape(2),shape(3))
      real*8 c213(shape(1),shape(2),shape(3))
      real*8 c222(shape(1),shape(2),shape(3))
      real*8 c223(shape(1),shape(2),shape(3))
      real*8 c233(shape(1),shape(2),shape(3))
      real*8 c311(shape(1),shape(2),shape(3))
      real*8 c312(shape(1),shape(2),shape(3))
      real*8 c313(shape(1),shape(2),shape(3))
      real*8 c322(shape(1),shape(2),shape(3))
      real*8 c323(shape(1),shape(2),shape(3))
      real*8 c333(shape(1),shape(2),shape(3))
      real*8 r11(shape(1),shape(2),shape(3))
      real*8 r12(shape(1),shape(2),shape(3))
      real*8 r13(shape(1),shape(2),shape(3))
      real*8 r22(shape(1),shape(2),shape(3))
      real*8 r23(shape(1),shape(2),shape(3))
      real*8 r33(shape(1),shape(2),shape(3))
      real*8 gup11(shape(1),shape(2),shape(3))
      real*8 gup12(shape(1),shape(2),shape(3))
      real*8 gup13(shape(1),shape(2),shape(3))
      real*8 gup22(shape(1),shape(2),shape(3))
      real*8 gup23(shape(1),shape(2),shape(3))
      real*8 gup33(shape(1),shape(2),shape(3))
		real*8 t11(shape(1),shape(2),shape(3))
		real*8 t12(shape(1),shape(2),shape(3))
		real*8 t13(shape(1),shape(2),shape(3))
		real*8 t21(shape(1),shape(2),shape(3))
		real*8 t22(shape(1),shape(2),shape(3))
		real*8 t23(shape(1),shape(2),shape(3))
		real*8 t31(shape(1),shape(2),shape(3))
		real*8 t32(shape(1),shape(2),shape(3))
		real*8 t33(shape(1),shape(2),shape(3))
		real*8 temp1(shape(1),shape(2),shape(3))
		integer nl

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

      nx = shape(1)
      ny = shape(2)
      nz = shape(3)

      call constraint(nx,ny,nz,mask,
     $     k11,k12,k13,k21,k22,k23,k31,k32,k33,
     $     hamnorm,mom1norm,mom2norm,mom3norm,
     $     c111,c112,c113,c122,c123,c133,
     $     c211,c212,c213,c222,c223,c233,
     $     c311,c312,c313,c322,c323,c333,
     $     r11,r12,r13,r22,r23,r33,
     $     gup11,gup12,gup13,gup22,gup23,gup33,
     $     t11,t12,t13,t21,t22,t23,t31,t32,t33,temp1,
     $     dx,m,d,nl)

      return
      end
C-----------------------------------------------------------------------------
      subroutine t2_inverse(lb,ub,shape,g11,g12,
     .                   g13,g22,g23,
     .                   g33,gup11,gup12,
     .                   gup13,gup22,gup23,
     .                   gup33,temp1,mask)
      implicit none
      integer lb(3),ub(3),shape(3)
      real*8 g11(shape(1),shape(2),shape(3))
      real*8 g12(shape(1),shape(2),shape(3))
      real*8 g13(shape(1),shape(2),shape(3))
      real*8 g22(shape(1),shape(2),shape(3))
      real*8 g23(shape(1),shape(2),shape(3))
      real*8 g33(shape(1),shape(2),shape(3))
      real*8 gup11(shape(1),shape(2),shape(3))
      real*8 gup12(shape(1),shape(2),shape(3))
      real*8 gup13(shape(1),shape(2),shape(3))
      real*8 gup22(shape(1),shape(2),shape(3))
      real*8 gup23(shape(1),shape(2),shape(3))
      real*8 gup33(shape(1),shape(2),shape(3))
      real*8 temp1(shape(1),shape(2),shape(3))
      real*8 mask(shape(1),shape(2),shape(3))
		integer nx,ny,nz

      nx = shape(1)
      ny = shape(2)
      nz = shape(3)

      call inverse(g11,g12,g13,g22,g23,g33,
     *     gup11,gup12,gup13,gup22, gup23, gup33,temp1,
     *     mask,nx, ny, nz)

      return
      end
C-----------------------------------------------------------------------------
		subroutine t2_connect(lb,ub,shape,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,
     .                      t11,t12,t13,
     .                      t21,t22,t23,
     .                      t31,t32,t33,
     .                      r11,r12,r13,
     .                      r22,r23,r33,
     .                      dx,nl,lbm,ubm,shapem,m,d)

		implicit none
		integer lb(3),ub(3),shape(3)
		integer lbm(3),ubm(3),shapem(3)

      real*8   g11(shape(1),shape(2),shape(3)),
     &         g12(shape(1),shape(2),shape(3)),
     &         g13(shape(1),shape(2),shape(3)),
     &         g22(shape(1),shape(2),shape(3)),
     &         g23(shape(1),shape(2),shape(3)),
     &         g33(shape(1),shape(2),shape(3))
      real*8   gup11(shape(1),shape(2),shape(3)),
     &         gup12(shape(1),shape(2),shape(3)),
     &         gup13(shape(1),shape(2),shape(3)),
     &         gup22(shape(1),shape(2),shape(3)),
     &         gup23(shape(1),shape(2),shape(3)),
     &         gup33(shape(1),shape(2),shape(3))

      real*8 c111(shape(1),shape(2),shape(3))
      real*8 c112(shape(1),shape(2),shape(3))
      real*8 c113(shape(1),shape(2),shape(3))
      real*8 c122(shape(1),shape(2),shape(3))
      real*8 c123(shape(1),shape(2),shape(3))
      real*8 c133(shape(1),shape(2),shape(3))
      real*8 c211(shape(1),shape(2),shape(3))
      real*8 c212(shape(1),shape(2),shape(3))
      real*8 c213(shape(1),shape(2),shape(3))
      real*8 c222(shape(1),shape(2),shape(3))
      real*8 c223(shape(1),shape(2),shape(3))
      real*8 c233(shape(1),shape(2),shape(3))
      real*8 c311(shape(1),shape(2),shape(3))
      real*8 c312(shape(1),shape(2),shape(3))
      real*8 c313(shape(1),shape(2),shape(3))
      real*8 c322(shape(1),shape(2),shape(3))
      real*8 c323(shape(1),shape(2),shape(3))
      real*8 c333(shape(1),shape(2),shape(3))

      real*8 r11(shape(1),shape(2),shape(3))
      real*8 r12(shape(1),shape(2),shape(3))
      real*8 r13(shape(1),shape(2),shape(3))
      real*8 r22(shape(1),shape(2),shape(3))
      real*8 r23(shape(1),shape(2),shape(3))
      real*8 r33(shape(1),shape(2),shape(3))

      real*8 t11(shape(1),shape(2),shape(3))
      real*8 t12(shape(1),shape(2),shape(3))
      real*8 t13(shape(1),shape(2),shape(3))
      real*8 t21(shape(1),shape(2),shape(3))
      real*8 t22(shape(1),shape(2),shape(3))
      real*8 t23(shape(1),shape(2),shape(3))
      real*8 t31(shape(1),shape(2),shape(3))
      real*8 t32(shape(1),shape(2),shape(3))
      real*8 t33(shape(1),shape(2),shape(3))

		integer nx,ny,nz
		real*8 dx
		integer nl
		integer m(*)
		real*8 d(*)
		integer l,i,j,k

      nx = shape(1)
      ny = shape(2)
      nz = shape(3)

      call connectv2(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,
     $     t11, t12, t13, t21, t22, t23, t31, t32, t33,
     $     r11, r12, r13, r22, r23, r33, dx, nx, ny, nz,nl,m,d)

		return
		end
C-----------------------------------------------------------------------------
		subroutine t2_riccidd(lb,ub,shape,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,
     .                      t11,t12,t13,
     .                      t21,t22,t23,
     .                      t31,t32,t33,
     .                      r11,r12,r13,
     .                      r22,r23,r33,
     .                      dx,nl,lbm,ubm,shapem,m,d)

		implicit none
		integer lb(3),ub(3),shape(3)
		integer lbm(3),ubm(3),shapem(3)

      real*8   g11(shape(1),shape(2),shape(3)),
     &         g12(shape(1),shape(2),shape(3)),
     &         g13(shape(1),shape(2),shape(3)),
     &         g22(shape(1),shape(2),shape(3)),
     &         g23(shape(1),shape(2),shape(3)),
     &         g33(shape(1),shape(2),shape(3))
      real*8   gup11(shape(1),shape(2),shape(3)),
     &         gup12(shape(1),shape(2),shape(3)),
     &         gup13(shape(1),shape(2),shape(3)),
     &         gup22(shape(1),shape(2),shape(3)),
     &         gup23(shape(1),shape(2),shape(3)),
     &         gup33(shape(1),shape(2),shape(3))

      real*8 c111(shape(1),shape(2),shape(3))
      real*8 c112(shape(1),shape(2),shape(3))
      real*8 c113(shape(1),shape(2),shape(3))
      real*8 c122(shape(1),shape(2),shape(3))
      real*8 c123(shape(1),shape(2),shape(3))
      real*8 c133(shape(1),shape(2),shape(3))
      real*8 c211(shape(1),shape(2),shape(3))
      real*8 c212(shape(1),shape(2),shape(3))
      real*8 c213(shape(1),shape(2),shape(3))
      real*8 c222(shape(1),shape(2),shape(3))
      real*8 c223(shape(1),shape(2),shape(3))
      real*8 c233(shape(1),shape(2),shape(3))
      real*8 c311(shape(1),shape(2),shape(3))
      real*8 c312(shape(1),shape(2),shape(3))
      real*8 c313(shape(1),shape(2),shape(3))
      real*8 c322(shape(1),shape(2),shape(3))
      real*8 c323(shape(1),shape(2),shape(3))
      real*8 c333(shape(1),shape(2),shape(3))

      real*8 r11(shape(1),shape(2),shape(3))
      real*8 r12(shape(1),shape(2),shape(3))
      real*8 r13(shape(1),shape(2),shape(3))
      real*8 r22(shape(1),shape(2),shape(3))
      real*8 r23(shape(1),shape(2),shape(3))
      real*8 r33(shape(1),shape(2),shape(3))

      real*8 t11(shape(1),shape(2),shape(3))
      real*8 t12(shape(1),shape(2),shape(3))
      real*8 t13(shape(1),shape(2),shape(3))
      real*8 t21(shape(1),shape(2),shape(3))
      real*8 t22(shape(1),shape(2),shape(3))
      real*8 t23(shape(1),shape(2),shape(3))
      real*8 t31(shape(1),shape(2),shape(3))
      real*8 t32(shape(1),shape(2),shape(3))
      real*8 t33(shape(1),shape(2),shape(3))

		integer nx,ny,nz
		real*8 dx
		integer nl
		integer m(*)
		real*8 d(*)

		nx = shape(1)
		ny = shape(2)
		nz = shape(3)

      call 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,
     $     t11, t12, t13, t21, t22, t23, t31, t32, t33,
     $     r11, r12, r13, r22, r23, r33,
     $     m,d,nl, dx, nx, ny, nz)

		return
		end
C--------------------------------------------------------------------------
		subroutine t2_evolve_g(lb,ub,shape,g11, g12, g13, g22, g23, g33,
     $                    gup11, gup12, gup13, gup22, gup23, gup33,
     $                    g11m1,g12m1,g13m1,g22m1,g23m1,g33m1,
     $                    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, dx,dt,nl,lbm,ubm,shapem,m,d) 

		implicit none
      integer lb(3),ub(3),shape(3)
		integer lbm(3),ubm(3),shapem(3)
      real*8   g11(shape(1),shape(2),shape(3)),
     &         g12(shape(1),shape(2),shape(3)),
     &         g13(shape(1),shape(2),shape(3)),
     &         g22(shape(1),shape(2),shape(3)),
     &         g23(shape(1),shape(2),shape(3)),
     &         g33(shape(1),shape(2),shape(3))
      real*8   g11m1(shape(1),shape(2),shape(3)),
     &         g12m1(shape(1),shape(2),shape(3)),
     &         g13m1(shape(1),shape(2),shape(3)),
     &         g22m1(shape(1),shape(2),shape(3)),
     &         g23m1(shape(1),shape(2),shape(3)),
     &         g33m1(shape(1),shape(2),shape(3))
      real*8   gup11(shape(1),shape(2),shape(3)),
     &         gup12(shape(1),shape(2),shape(3)),
     &         gup13(shape(1),shape(2),shape(3)),
     &         gup22(shape(1),shape(2),shape(3)),
     &         gup23(shape(1),shape(2),shape(3)),
     &         gup33(shape(1),shape(2),shape(3))
c     extrinsic curvature:  K^i_j
      real*8   k11(shape(1),shape(2),shape(3)),
     &         k12(shape(1),shape(2),shape(3)),
     &         k13(shape(1),shape(2),shape(3)),
     &         k21(shape(1),shape(2),shape(3)),
     &         k22(shape(1),shape(2),shape(3)),
     &         k23(shape(1),shape(2),shape(3)),
     &         k31(shape(1),shape(2),shape(3)),
     &         k32(shape(1),shape(2),shape(3)),
     &         k33(shape(1),shape(2),shape(3))


      real*8 c111(shape(1),shape(2),shape(3))
      real*8 c112(shape(1),shape(2),shape(3))
      real*8 c113(shape(1),shape(2),shape(3))
      real*8 c122(shape(1),shape(2),shape(3))
      real*8 c123(shape(1),shape(2),shape(3))
      real*8 c133(shape(1),shape(2),shape(3))
      real*8 c211(shape(1),shape(2),shape(3))
      real*8 c212(shape(1),shape(2),shape(3))
      real*8 c213(shape(1),shape(2),shape(3))
      real*8 c222(shape(1),shape(2),shape(3))
      real*8 c223(shape(1),shape(2),shape(3))
      real*8 c233(shape(1),shape(2),shape(3))
      real*8 c311(shape(1),shape(2),shape(3))
      real*8 c312(shape(1),shape(2),shape(3))
      real*8 c313(shape(1),shape(2),shape(3))
      real*8 c322(shape(1),shape(2),shape(3))
      real*8 c323(shape(1),shape(2),shape(3))
      real*8 c333(shape(1),shape(2),shape(3))
      real*8 r11(shape(1),shape(2),shape(3))
      real*8 r12(shape(1),shape(2),shape(3))
      real*8 r13(shape(1),shape(2),shape(3))
      real*8 r22(shape(1),shape(2),shape(3))
      real*8 r23(shape(1),shape(2),shape(3))
      real*8 r33(shape(1),shape(2),shape(3))
      real*8 alpha(shape(1),shape(2),shape(3))
      real*8 beta1(shape(1),shape(2),shape(3))
      real*8 beta2(shape(1),shape(2),shape(3))
      real*8 beta3(shape(1),shape(2),shape(3))
      real*8 temp1(shape(1),shape(2),shape(3))
      real*8 temp2(shape(1),shape(2),shape(3))
      real*8 temp3(shape(1),shape(2),shape(3))
      real*8 t11(shape(1),shape(2),shape(3))
      real*8 t12(shape(1),shape(2),shape(3))
      real*8 t13(shape(1),shape(2),shape(3))
      real*8 t21(shape(1),shape(2),shape(3))
      real*8 t22(shape(1),shape(2),shape(3))
      real*8 t23(shape(1),shape(2),shape(3))
      real*8 t31(shape(1),shape(2),shape(3))
      real*8 t32(shape(1),shape(2),shape(3))
      real*8 t33(shape(1),shape(2),shape(3))
		real*8 dx,dt

      integer nl
      integer m(*)
      real*8 d(*)
		integer nx,ny,nz

		nx = shape(1)
		ny = shape(2)
		nz = shape(3)

      call evolve_g(g11, g12, g13, g22, g23, g33,
     $              gup11, gup12, gup13, gup22, gup23, gup33,
     $              g11m1,g12m1,g13m1,g22m1,g23m1,g33m1,
     $              g11m1,g12m1,g13m1,g22m1,g23m1,g33m1,
     $              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, dx, dt,nx, ny, nz)

		return
		end
C--------------------------------------------------------------------------
		subroutine t2_evolve_ka(lb,ub,shape,g11, g12, g13, g22, g23, g33,
     $                    gup11, gup12, gup13, gup22, gup23, gup33,
     $                    k11m1, k12m1, k13m1,
     $                    k21m1, k22m1, k23m1,
     $                    k31m1, k32m1, k33m1,
     $                    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, dx,dt,nl,lbm,ubm,shapem,m,d) 

		implicit none
		integer nx,ny,nz
      integer lb(3),ub(3),shape(3)
		integer lbm(3),ubm(3),shapem(3)
      real*8   g11(shape(1),shape(2),shape(3)),
     &         g12(shape(1),shape(2),shape(3)),
     &         g13(shape(1),shape(2),shape(3)),
     &         g22(shape(1),shape(2),shape(3)),
     &         g23(shape(1),shape(2),shape(3)),
     &         g33(shape(1),shape(2),shape(3))

      real*8   k11m1(shape(1),shape(2),shape(3)),
     &         k12m1(shape(1),shape(2),shape(3)),
     &         k13m1(shape(1),shape(2),shape(3)),
     &         k21m1(shape(1),shape(2),shape(3)),
     &         k22m1(shape(1),shape(2),shape(3)),
     &         k23m1(shape(1),shape(2),shape(3)),
     &         k31m1(shape(1),shape(2),shape(3)),
     &         k32m1(shape(1),shape(2),shape(3)),
     &         k33m1(shape(1),shape(2),shape(3))

      real*8   gup11(shape(1),shape(2),shape(3)),
     &         gup12(shape(1),shape(2),shape(3)),
     &         gup13(shape(1),shape(2),shape(3)),
     &         gup22(shape(1),shape(2),shape(3)),
     &         gup23(shape(1),shape(2),shape(3)),
     &         gup33(shape(1),shape(2),shape(3))
c     extrinsic curvature:  K^i_j
      real*8   k11(shape(1),shape(2),shape(3)),
     &         k12(shape(1),shape(2),shape(3)),
     &         k13(shape(1),shape(2),shape(3)),
     &         k21(shape(1),shape(2),shape(3)),
     &         k22(shape(1),shape(2),shape(3)),
     &         k23(shape(1),shape(2),shape(3)),
     &         k31(shape(1),shape(2),shape(3)),
     &         k32(shape(1),shape(2),shape(3)),
     &         k33(shape(1),shape(2),shape(3))

      real*8 c111(shape(1),shape(2),shape(3))
      real*8 c112(shape(1),shape(2),shape(3))
      real*8 c113(shape(1),shape(2),shape(3))
      real*8 c122(shape(1),shape(2),shape(3))
      real*8 c123(shape(1),shape(2),shape(3))
      real*8 c133(shape(1),shape(2),shape(3))
      real*8 c211(shape(1),shape(2),shape(3))
      real*8 c212(shape(1),shape(2),shape(3))
      real*8 c213(shape(1),shape(2),shape(3))
      real*8 c222(shape(1),shape(2),shape(3))
      real*8 c223(shape(1),shape(2),shape(3))
      real*8 c233(shape(1),shape(2),shape(3))
      real*8 c311(shape(1),shape(2),shape(3))
      real*8 c312(shape(1),shape(2),shape(3))
      real*8 c313(shape(1),shape(2),shape(3))
      real*8 c322(shape(1),shape(2),shape(3))
      real*8 c323(shape(1),shape(2),shape(3))
      real*8 c333(shape(1),shape(2),shape(3))
      real*8 r11(shape(1),shape(2),shape(3))
      real*8 r12(shape(1),shape(2),shape(3))
      real*8 r13(shape(1),shape(2),shape(3))
      real*8 r22(shape(1),shape(2),shape(3))
      real*8 r23(shape(1),shape(2),shape(3))
      real*8 r33(shape(1),shape(2),shape(3))
      real*8 alpha(shape(1),shape(2),shape(3))
      real*8 beta1(shape(1),shape(2),shape(3))
      real*8 beta2(shape(1),shape(2),shape(3))
      real*8 beta3(shape(1),shape(2),shape(3))
      real*8 temp1(shape(1),shape(2),shape(3))
      real*8 temp2(shape(1),shape(2),shape(3))
      real*8 temp3(shape(1),shape(2),shape(3))
      real*8 t11(shape(1),shape(2),shape(3))
      real*8 t12(shape(1),shape(2),shape(3))
      real*8 t13(shape(1),shape(2),shape(3))
      real*8 t21(shape(1),shape(2),shape(3))
      real*8 t22(shape(1),shape(2),shape(3))
      real*8 t23(shape(1),shape(2),shape(3))
      real*8 t31(shape(1),shape(2),shape(3))
      real*8 t32(shape(1),shape(2),shape(3))
      real*8 t33(shape(1),shape(2),shape(3))

      integer nl
      integer m(*)
      real*8 d(*)
		real*8 dx,dt

		nx = shape(1)
		ny = shape(2)
		nz = shape(3)

      call evolve_ka(g11, g12, g13, g22, g23, g33,
     $              gup11, gup12, gup13, gup22, gup23, gup33,
     $              k11m1,k12m1,k13m1,
     $              k21m1,k22m1,k23m1,
     $              k31m1,k32m1,k33m1,
     $              k11m1,k12m1,k13m1,
     $              k21m1,k22m1,k23m1,
     $              k31m1,k32m1,k33m1,
     $              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, dx, dt,nx, ny, nz)

		return
		end
C--------------------------------------------------------------------------
		subroutine t2_evolve_kb(lb,ub,shape,g11, g12, g13, g22, g23, g33,
     $                    gup11, gup12, gup13, gup22, gup23, gup33,
     $                    k11m1, k12m1, k13m1,
     $                    k21m1, k22m1, k23m1,
     $                    k31m1, k32m1, k33m1,
     $                    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, dx,dt,nl,lbm,ubm,shapem,m,d) 

		implicit none
		integer nx,ny,nz
      integer lb(3),ub(3),shape(3)
		integer lbm(3),ubm(3),shapem(3)
      real*8   g11(shape(1),shape(2),shape(3)),
     &         g12(shape(1),shape(2),shape(3)),
     &         g13(shape(1),shape(2),shape(3)),
     &         g22(shape(1),shape(2),shape(3)),
     &         g23(shape(1),shape(2),shape(3)),
     &         g33(shape(1),shape(2),shape(3))

      real*8   k11m1(shape(1),shape(2),shape(3)),
     &         k12m1(shape(1),shape(2),shape(3)),
     &         k13m1(shape(1),shape(2),shape(3)),
     &         k21m1(shape(1),shape(2),shape(3)),
     &         k22m1(shape(1),shape(2),shape(3)),
     &         k23m1(shape(1),shape(2),shape(3)),
     &         k31m1(shape(1),shape(2),shape(3)),
     &         k32m1(shape(1),shape(2),shape(3)),
     &         k33m1(shape(1),shape(2),shape(3))

      real*8   gup11(shape(1),shape(2),shape(3)),
     &         gup12(shape(1),shape(2),shape(3)),
     &         gup13(shape(1),shape(2),shape(3)),
     &         gup22(shape(1),shape(2),shape(3)),
     &         gup23(shape(1),shape(2),shape(3)),
     &         gup33(shape(1),shape(2),shape(3))
c     extrinsic curvature:  K^i_j
      real*8   k11(shape(1),shape(2),shape(3)),
     &         k12(shape(1),shape(2),shape(3)),
     &         k13(shape(1),shape(2),shape(3)),
     &         k21(shape(1),shape(2),shape(3)),
     &         k22(shape(1),shape(2),shape(3)),
     &         k23(shape(1),shape(2),shape(3)),
     &         k31(shape(1),shape(2),shape(3)),
     &         k32(shape(1),shape(2),shape(3)),
     &         k33(shape(1),shape(2),shape(3))

      real*8 c111(shape(1),shape(2),shape(3))
      real*8 c112(shape(1),shape(2),shape(3))
      real*8 c113(shape(1),shape(2),shape(3))
      real*8 c122(shape(1),shape(2),shape(3))
      real*8 c123(shape(1),shape(2),shape(3))
      real*8 c133(shape(1),shape(2),shape(3))
      real*8 c211(shape(1),shape(2),shape(3))
      real*8 c212(shape(1),shape(2),shape(3))
      real*8 c213(shape(1),shape(2),shape(3))
      real*8 c222(shape(1),shape(2),shape(3))
      real*8 c223(shape(1),shape(2),shape(3))
      real*8 c233(shape(1),shape(2),shape(3))
      real*8 c311(shape(1),shape(2),shape(3))
      real*8 c312(shape(1),shape(2),shape(3))
      real*8 c313(shape(1),shape(2),shape(3))
      real*8 c322(shape(1),shape(2),shape(3))
      real*8 c323(shape(1),shape(2),shape(3))
      real*8 c333(shape(1),shape(2),shape(3))
      real*8 r11(shape(1),shape(2),shape(3))
      real*8 r12(shape(1),shape(2),shape(3))
      real*8 r13(shape(1),shape(2),shape(3))
      real*8 r22(shape(1),shape(2),shape(3))
      real*8 r23(shape(1),shape(2),shape(3))
      real*8 r33(shape(1),shape(2),shape(3))
      real*8 alpha(shape(1),shape(2),shape(3))
      real*8 beta1(shape(1),shape(2),shape(3))
      real*8 beta2(shape(1),shape(2),shape(3))
      real*8 beta3(shape(1),shape(2),shape(3))
      real*8 temp1(shape(1),shape(2),shape(3))
      real*8 temp2(shape(1),shape(2),shape(3))
      real*8 temp3(shape(1),shape(2),shape(3))
      real*8 t11(shape(1),shape(2),shape(3))
      real*8 t12(shape(1),shape(2),shape(3))
      real*8 t13(shape(1),shape(2),shape(3))
      real*8 t21(shape(1),shape(2),shape(3))
      real*8 t22(shape(1),shape(2),shape(3))
      real*8 t23(shape(1),shape(2),shape(3))
      real*8 t31(shape(1),shape(2),shape(3))
      real*8 t32(shape(1),shape(2),shape(3))
      real*8 t33(shape(1),shape(2),shape(3))

      integer nl
      integer m(*)
      real*8 d(*)
		real*8 dx,dt

		nx = shape(1)
		ny = shape(2)
		nz = shape(3)

      call evolve_kb(g11, g12, g13, g22, g23, g33,
     $              gup11, gup12, gup13, gup22, gup23, gup33,
     $              k11m1,k12m1,k13m1,
     $              k21m1,k22m1,k23m1,
     $              k31m1,k32m1,k33m1,
     $              k11m1,k12m1,k13m1,
     $              k21m1,k22m1,k23m1,
     $              k31m1,k32m1,k33m1,
     $              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, dx, dt,nx, ny, nz)

		return
		end
c-----------------------------------------------------------------------
c  Frontend to Teukolsky wave initial data
c
		subroutine t2_linearwave(lb,ub,shape,
     *      g11, g12, g13, g22, g23, g33,
     *      k11, k12, k13, k21, k22, k23, k31, k32, k33,
     *      xx, yy, zz, rr, A, B, C, 
     *      bb, dx, dy, dz, amp, width)
		implicit none

      integer lb(3),ub(3),shape(3)

      real*8   g11(shape(1),shape(2),shape(3)),
     &         g12(shape(1),shape(2),shape(3)),
     &         g13(shape(1),shape(2),shape(3)),
     &         g22(shape(1),shape(2),shape(3)),
     &         g23(shape(1),shape(2),shape(3)),
     &         g33(shape(1),shape(2),shape(3))

c     extrinsic curvature:  K^i_j
      real*8   k11(shape(1),shape(2),shape(3)),
     &         k12(shape(1),shape(2),shape(3)),
     &         k13(shape(1),shape(2),shape(3)),
     &         k21(shape(1),shape(2),shape(3)),
     &         k22(shape(1),shape(2),shape(3)),
     &         k23(shape(1),shape(2),shape(3)),
     &         k31(shape(1),shape(2),shape(3)),
     &         k32(shape(1),shape(2),shape(3)),
     &         k33(shape(1),shape(2),shape(3))

      real*8   xx(shape(1),shape(2),shape(3)),
     &         yy(shape(1),shape(2),shape(3)),
     &         zz(shape(1),shape(2),shape(3)),
     &         rr(shape(1),shape(2),shape(3))

      real*8   A(shape(1),shape(2),shape(3)),
     &         B(shape(1),shape(2),shape(3)),
     &         C(shape(1),shape(2),shape(3))

		real*8	bb(6), dx, dy, dz, amp, width
		  call makegrid3d(shape(1),shape(2),shape(3),
     &           bb(1),bb(2),bb(3),bb(4),bb(5),bb(6),
     &           dx,dy,dz,xx,yy,zz,rr,lb)

		  call linearwave(shape(1), shape(2), shape(3), amp, width,
     *           xx, yy, zz, rr, g11, g12, g13, g22, g23, g33,
     *           k11, k12, k13, k21, k22, k23, k31, k32, k33,
     *           A, B, C)

		  return
		end
