!*********************************************************************** 
!     This code evolves binary black hole spacetimes using the ADM 3+1 
!  formalism.  A initial spacelike hypersurface is formed by using 
!  an exact form for a known spacetime, or by solving the time-independent 
!  Einstein equations, also known as the constraint equations.  This 
!  initial hypersurface is evolved using the time-dependent Einstein 
!  equations.  In this manner, a foliation of the full spacetime can 
!  be computed. 
!     The hypersurfaces are freely evolved, i.e.--the constraint 
!  equations are not used in the solution of hypersurfaces subsequent 
!  to the initial one.  The evolution is a three-level update scheme, 
!  using time-centered and space-centered derivatives.  At the boundaries 
!  of the computational domain, the outer edges of the grid and the 
!  boundary of the excised singularity, we use one-sided differences 
!  to compute spatial derivatives. 
!     Two important subroutines used are the apparent 
!  horizon solver and a generic interpolator/extrapolator.  The program 
!  is designed in a modular fasion to allow a variety of finite differencing 
!  schemes and conforms to a DAGH compatible data structure. 
! 
!     Randall R. Correll 
!     Mijan F. Huq 
!     Scott A. Klasky 
!     Richard A. Matzner 
! 
!     COPYRIGHT 1995 - 1996 
! 
!     Rewritten to use symmetric K_{ij} 1996 Mijan Huq 
!*********************************************************************** 
 
program t2 
implicit  none 
 
!  Problem parameters 
 
integer nx, ny, nz, nl, nlmax, nplot 
parameter (nlmax = 1000)
 
real*8    tmin, tmax 
real*8,dimension(100)    :: params 
 
!  Black hole parameters 
      
 
!  Grid specifications 
real*8,dimension(:),allocatable    :: x,y,z
real*8    h
real*8    xns

 
!  Main loop parameters 
 
integer   :: it,nt,itsave 
real*8    :: dt, time 
real*8    :: dtsave, timesave 
real*8    :: tol 
 
!  Dynamic field variables 
real*8,dimension(:,:,:,:),allocatable :: g, gm1, gup, k, km1
 
!  Gauge variables 
real*8,dimension(:,:,:,:),allocatable :: beta
real*8,dimension(:,:,:),allocatable   :: alpha 
 
!  Spatial tensors 
 
real*8,dimension(:,:,:,:),allocatable :: r
real*8,dimension(:,:,:,:),allocatable :: c

real*8,dimension(:,:,:,:),allocatable :: ca, rieman
 
!  Temporary arrays for calculations within subroutines 
real*8,dimension(:,:,:,:),allocatable :: t
 
!  Variables to output constraint monitors 
real*8                              :: hamnorm, mom1norm, mom2norm, mom3norm 
 
!  Characteristic functions 
real*8,dimension(:,:,:),allocatable :: mask,maskold
real*8,dimension(:)    ,allocatable :: d ,fx,fy,fz
integer,dimension(:)   ,allocatable :: m,xi,yj,zk

!real*8                              :: n2i
!--------------------------------------------------------------------------
interface
   subroutine input(params)
      implicit none
      real*8,dimension(:)           :: params      
      integer                       :: nx, ny, nz, nt, nplot
      real*8                        :: xmin, xmax, ymin, ymax,zmin, zmax, tol
      integer                       :: hdfflag, asciiflag
      real*8                        :: bha(2),bhr(2,3),bhp(2,3),bhs(2,3)
      real*8,parameter              :: UNDEF = -1.0d0
      integer                       :: nxi, nyi, nzi
   end subroutine
   subroutine initialize(g,gm1,gup,km1,k,c,r,ca,rieman,alpha,beta,t,x,y,z,xi,&
                         yj,zk,fx,fy,fz,mask,maskold,m,d,nl,params,h,dt)
      implicit none
      integer                       :: nl
      real*8                        :: h,dt
      real*8,dimension(:,:,:,:)     :: g,gup,gm1,k,km1,c,t,r,beta,ca,rieman
      real*8,dimension(:,:,:)       :: alpha,mask,maskold
      real*8,dimension(:)           :: params
      real*8,dimension(:)           :: d
      integer,dimension(:)          :: m
      real*8,dimension(:)          :: fx,fy,fz
      integer,dimension(:)         :: xi,yj,zk
      real*8,dimension(:)           :: x
      real*8,dimension(:)           :: y
      real*8,dimension(:)           :: z
      real*8                        :: dtsave,timesave
   end subroutine
   subroutine swap(a,b,temp)
      implicit none
      real*8,dimension(:,:,:,:)    :: a,b,temp
   end subroutine
   subroutine h_constraint(g,k,r,gup,t,mask,hc)
      implicit none
      real*8,dimension(:,:,:)             :: mask ,hc
      real*8,dimension(:,:,:,:),target    :: g,gup,k,r,t
   end subroutine
   subroutine k_constraint(g,k,gup,c,h,t,mask,m,d,nl,kc)
      implicit none
      integer                             :: nl
      real*8,dimension(:,:,:)             :: mask
      real*8,dimension(:,:,:,:),target    :: g,gup,k,t,c,kc
      real*8,dimension(:)              :: d
      integer,dimension(:)             :: m
      real*8                              :: h
   end subroutine
   subroutine takestep(g,gup,gm1,k,km1,h,t,r,c,ca,rieman,alpha,beta,      &
                       mask,m,d,nl,dt,gp1,kp1)
      implicit none
      integer                       :: nl
      real*8                        :: h,dt
      real*8,dimension(:,:,:,:)     :: g,gup,gm1,k,km1
      real*8,dimension(:,:,:,:)     :: t,r,c,ca,rieman,beta
      real*8,dimension(:,:,:,:)     :: gp1,kp1
      real*8,dimension(:,:,:)       :: alpha,mask
      real*8,dimension(:)        :: d
      integer,dimension(:)       :: m
   end subroutine
   subroutine horizon(g,k,xx,yy,zz,t,params,mask)
      implicit none
      real*8,dimension(:,:,:,:),target     :: g,k,t
      real*8,dimension(:,:,:)              :: mask
      real*8,dimension(:)                  :: params
      real*8,dimension(:)                  :: xx,yy,zz
      real*8                               :: test,h
      integer                              :: i
      real*8,dimension(:,:,:),pointer      :: x,y,z,r
   end subroutine
   subroutine lapse(x,y,z,g,k,gup,t,mask,params,alpha)
      implicit none
      integer                              :: nl
      real*8,dimension(:,:,:,:),target     :: g,k,gup,t
      real*8,dimension(:,:,:)              :: mask
      real*8,dimension(:,:,:)              :: alpha
      real*8,dimension(:)                  :: x,y,z
      real*8,dimension(:)                  :: params
      real*8,dimension(:,:,:),pointer      :: t1,t2,t3,t4
      real*8                               :: bha(2),bhr(3,2)
   end subroutine
   subroutine shift(x,y,z,g,k,gup,t,mask,params,beta)
      implicit none
      real*8,dimension(:,:,:,:),target     :: g,k,gup,t
      real*8,dimension(:,:,:)              :: mask
      real*8,dimension(:,:,:,:)            :: beta
      real*8,dimension(:)                  :: x,y,z
      real*8,dimension(:)                  :: params
      real*8,dimension(:,:,:),pointer      :: t1,t2,t3,t4
      real*8                               :: bha(2),bhr(3,2)
   end subroutine

   real*8 function n2i(n,a,mask,m,nl)
      implicit none
      integer                       :: n,nl
      real*8,dimension(:,:,:)       :: a,mask
      integer                       :: i,j,k,l
      integer,dimension(:)       :: m
end function

   subroutine remask(g,k,mask,maskold)
      implicit none
      integer                       :: nx, ny, nz
      real*8,dimension(:,:,:,:)     :: g,k
      real*8,dimension(:,:,:)       :: mask,maskold
   end subroutine

   subroutine flags(params,xi,yj,zk,fx,fy,fz,mask,m,d,nl)
      implicit none
      integer                       :: nx,ny,nz,nl
      real*8,dimension(:)           :: fx,fy,fz
      real*8,dimension(:)           :: params
      integer,dimension(:)          :: xi,yj,zk
      real*8,dimension(:)           :: d
      integer,dimension(:)          :: m
      real*8,dimension(:,:,:)       :: mask
   end subroutine

end interface
!--------------------------------------------------------------------------
 
 
! * * * * * * * * * * * * * Read Input file * * * * * * * * * * * * *  
   call input(params) 

   tol = params(22)
   nt  = int(params(23))
   nx  = int(params(1))
   ny  = int(params(2))
   nz  = int(params(3))

   allocate(x(nx),y(ny),z(nz))
   allocate(g(6,nx,ny,nz),gup(6,nx,ny,nz),gm1(6,nx,ny,nz))
   allocate(k(6,nx,ny,nz),km1(6,nx,ny,nz))
   allocate(r(6,nx,ny,nz),beta(6,nx,ny,nz))
   allocate(alpha(nx,ny,nz))
   allocate(c(18,nx,ny,nz),t(10,nx,ny,nz))
   allocate(mask(nx,ny,nz),maskold(nx,ny,nz))
   allocate (d(3*nlmax),m(3*nlmax),fx(nlmax),fy(nlmax),fz(nlmax))
   allocate (xi(nlmax),yj(nlmax),zk(nlmax))

! dummy arrays in this version of the adm code
   allocate(ca(1,1,1,1),rieman(1,1,1,1))


! * * * * * * * * * Initialize n and n-1 levels * * * * * * * * * * *  
   call initialize(g,gm1,gup,km1,k,c,r,ca,rieman,alpha,beta,t,x,y,z,xi,   &
                   yj,zk,fx,fy,fz,mask,maskold,m,d,nl,params,h,dt)

 
! * * * * * * * * * *Evolution Loop * * * * * * * * * * * * * * * * *  
   write(6,*) 'Begin TCODE run.' 
 
 
   time  = 0.0d0 
   it = 0 
 
   write(6,*) '...Begin main loop, forward in time' 
 
 
   write(6,*) 'Iter  Time  HamErr  Mom1Err  Mom2Err  Mom3Err' 
 
   call h_constraint(g,k,r,gup,t,mask,t(7,:,:,:))
   call k_constraint(g,k,gup,c,h,t,mask,m,d,nl,t(:,:,:,8:10))

write(6,*) 'before call'
write(6,*) t(7,3,3,3:5)

   hamnorm  = n2i(nx,t(7,:,:,:),mask,m,nl)
stop
   mom1norm = n2i(nx,t(8,:,:,:),mask,m,nl)
   mom2norm = n2i(nx,t(9,:,:,:),mask,m,nl)
   mom3norm = n2i(nx,t(10,:,:,:),mask,m,nl)

   write(6,'(i5,f9.5,1P,4e12.3)')it,time,hamnorm,mom1norm,mom2norm,mom3norm 
   write(2,'(i5,f9.5,1P,4e12.3)')it,time,hamnorm,mom1norm,mom2norm,mom3norm 
 
   do it = 1, nt 
 
      call lapse(x,y,z,g,k,gup,t,mask,params,alpha)
      call shift(x,y,z,g,k,gup,t,mask,params,beta)
      dt = h / 4.0d0
      time = time + dt 
 
      call  takestep(g,gup,gm1,k,km1,h,t,r,c,ca,rieman,alpha,beta,       &
                     mask,m,d,nl,dt,gm1,km1)

      call swap(g,gm1,t)
      call swap(k,km1,t)

      maskold = mask
 
      call horizon(g,k,x,y,z,t,params,mask)
 
      call remask(g,k,mask,maskold) 
      call remask(gm1,km1,mask,maskold)
 
      call flags(params,xi,yj,zk,fx,fy,fz,mask,m,d,nl) 
 
      call h_constraint(g,k,r,gup,t,mask,t(7,:,:,:))
      call k_constraint(g,k,gup,c,h,t,mask,m,d,nl,t(:,:,:,8:10))

      hamnorm  = n2i(nx,t(7,:,:,:),mask,m,nl)
      mom1norm = n2i(nx,t(8,:,:,:),mask,m,nl)
      mom2norm = n2i(nx,t(9,:,:,:),mask,m,nl)
      mom3norm = n2i(nx,t(10,:,:,:),mask,m,nl)
 
      write(6,'(i5,f9.5,1P,4e12.3)') it,time,hamnorm,mom1norm,mom2norm,mom3norm 
 
      if(hamnorm .gt. tol)goto 1919 
 
   end do 
 
1919   continue 
 
stop 
end 
!----------------------------------------------------------------------- 
subroutine initialize(g,gm1,gup,km1,k,c,r,ca,rieman,alpha,beta,t,x,y,z,xi,  &
                      yj,zk,fx,fy,fz,mask,maskold,m,d,nl,params,h,dt) 
implicit none
integer                       :: nl
real*8                        :: h,dt 
real*8,dimension(:,:,:,:)     :: g,gup,gm1,k,km1,c,t,r,beta,ca,rieman
real*8,dimension(:,:,:)       :: alpha,mask,maskold
real*8,dimension(:)           :: params

real*8,dimension(:)           :: d
integer,dimension(:)          :: m
real*8,dimension(:)          :: fx,fy,fz 
integer,dimension(:)         :: xi,yj,zk 
real*8,dimension(:)           :: x
real*8,dimension(:)           :: y
real*8,dimension(:)           :: z

real*8                        :: dtsave,timesave 
!--------------------------------------------------------------------------
interface
   subroutine makegrid(params,h,x,y,z)
      implicit none
      real*8                       :: h
      real*8,dimension(:)          :: params
      real*8,dimension(:)          :: x,y,z
      integer                       :: i
   end subroutine
   subroutine initial(x,y,z,params,t,g,k)
      implicit none
      real*8,dimension(:,:,:,:)    :: g,k,t
      real*8,dimension(:)          :: x,y,z
      real*8,dimension(:)          :: params
   end subroutine
   subroutine remask(g,k,mask,maskold)
      implicit none
      integer                       :: nx, ny, nz
      real*8,dimension(:,:,:,:)     :: g,k
      real*8,dimension(:,:,:)       :: mask,maskold
   end subroutine
   subroutine flags(params,xi,yj,zk,fx,fy,fz,mask,m,d,nl)
      implicit none
      integer                       :: nx,ny,nz,nl
      real*8,dimension(:)           :: fx,fy,fz
      real*8,dimension(:)           :: params
      integer,dimension(:)          :: xi,yj,zk
      real*8,dimension(:)           :: d
      integer,dimension(:)          :: m
      real*8,dimension(:,:,:)       :: mask
      integer                      :: i,j,k,l,nlxmax,nlymax,nlzmax
      real*8 ,parameter            :: zero = 0.0d0
      real*8 ,parameter            :: one  = 1.0d0
   end subroutine
   subroutine lapse(x,y,z,g,k,gup,t,mask,params,alpha)
      implicit none
      real*8,dimension(:,:,:,:),target     :: g,k,gup,t
      real*8,dimension(:,:,:)              :: mask
      real*8,dimension(:,:,:)              :: alpha
      real*8,dimension(:)                  :: x,y,z
      real*8,dimension(:)                  :: params
      real*8,dimension(:,:,:),pointer      :: t1,t2,t3,t4
      real*8                               :: bha(2),bhr(3,2)
   end subroutine
   subroutine shift(x,y,z,g,k,gup,t,mask,params,beta)
      implicit none
      real*8,dimension(:,:,:,:),target     :: g,k,gup,t
      real*8,dimension(:,:,:)              :: mask
      real*8,dimension(:,:,:,:)            :: beta
      real*8,dimension(:)                  :: x,y,z
      real*8,dimension(:)                  :: params
      real*8,dimension(:,:,:),pointer      :: t1,t2,t3,t4
      real*8                               :: bha(2),bhr(3,2)
   end subroutine
   subroutine takestep(g,gup,gm1,k,km1,h,t,r,c,ca,rieman,alpha,beta,      &
                       mask,m,d,nl,dt,gp1,kp1)
      implicit none
      integer                       :: nl
      real*8                        :: h,dt
      real*8,dimension(:,:,:,:)     :: g,gup,gm1,k,km1
      real*8,dimension(:,:,:,:)     :: t,r,c,ca,rieman,beta
      real*8,dimension(:,:,:,:)     :: gp1,kp1
      real*8,dimension(:,:,:)       :: alpha,mask
      real*8,dimension(:)        :: d
      integer,dimension(:)       :: m
   end subroutine
   subroutine h_constraint(g,k,r,gup,t,mask,hc)
      implicit none
      real*8,dimension(:,:,:)             :: mask ,hc
      real*8,dimension(:,:,:,:),target    :: g,gup,k,r,t
   end subroutine
   subroutine k_constraint(g,k,gup,c,h,t,mask,m,d,nl,kc)
      implicit none
      integer                             :: nl
      real*8,dimension(:,:,:)             :: mask
      real*8,dimension(:,:,:,:),target    :: g,gup,k,t,c,kc
      real*8,dimension(:)              :: d
      integer,dimension(:)             :: m
      real*8                              :: h
   end subroutine
   subroutine swap(a,b,temp)
      implicit none
      real*8,dimension(:,:,:,:)    :: a,b,temp
   end subroutine
   subroutine horizon(g,k,xx,yy,zz,t,params,mask)
      implicit none
      real*8,dimension(:,:,:,:),target     :: g,k,t
      real*8,dimension(:,:,:)              :: mask
      real*8,dimension(:)                  :: params
      real*8,dimension(:)                  :: xx,yy,zz
      real*8                               :: test,h
      integer                              :: i
      real*8,dimension(:,:,:),pointer      :: x,y,z,r
   end subroutine

end interface
!--------------------------------------------------------------------------
 
 
   h = (params(16)-params(15))/(params(1)-1) 
 
   write(6,*) 'Generating coordinates' 
   call makegrid(params,h,x,y,z) 
 
   write(6,*)'Set up initial data' 
   call initial(x,y,z,params,t,g,k)
 
   write(6,*) 'Generate initial horizon mask' 

   call horizon(g,k,x,y,z,t,params,mask)

   maskold = mask
 
   call remask(g,k,mask,maskold) 
 
   write(6,*) '...Compute excise mask and boundary points' 
   
   call flags(params,xi,yj,zk,fx,fy,fz,mask,m,d,nl)

   write(6,*) 'nl = ',nl
 
   call lapse(x,y,z,g,k,gup,t,mask,params,alpha)
   call shift(x,y,z,g,k,gup,t,mask,params,beta)
 
   dt = h / 4.0d0 
 
   dtsave = dt 
 
   dt = -0.5d0*dt 
 
!         Set up n+1/2 step 
   gm1 = g
   km1 = k
!         Evolve to n-1/2 step 
   call  takestep(g,gup,gm1,k,km1,h,t,r,c,ca,rieman,alpha,beta,       &
                  mask,m,d,nl,dt,gm1,km1)
 
   call swap(g,gm1,t)
   call swap(k,km1,t)
 
!         Evolve from n-1/2 to n-1 step 
   call  takestep(g,gup,gm1,k,km1,h,t,r,c,ca,rieman,alpha,beta,       &
                  mask,m,d,nl,dt,gm1,km1)
!            Leave m1 data in place. Do not swap 
 
!  Then continue with standard three-level scheme. 
 
   write(6,*) 'Completed initial backward timestep' 
 
!  Restore times and time levels 
 
   dt = dtsave 
 
   write(6,*) '****Re-initialize t=0 data****' 
   call initial(x,y,z,params,t,g,k)

   write(6,*)'****Done setting up n,n-1 data****' 
 
 
return 
end    
!----------------------------------------------------------------------- 
! 
! Routine to take one evolution time step 
! 
subroutine takestep(g,gup,gm1,k,km1,h,t,r,c,ca,rieman,alpha,beta,      &
                    mask,m,d,nl,dt,gp1,kp1)
implicit none 
integer                       :: nl
real*8                        :: h,dt 
real*8,dimension(:,:,:,:)     :: g,gup,gm1,k,km1
real*8,dimension(:,:,:,:)     :: t,r,c,ca,rieman,beta 
real*8,dimension(:,:,:,:)     :: gp1,kp1 
real*8,dimension(:,:,:)       :: alpha,mask 
real*8                        :: xns
real*8,dimension(:)        :: d
integer,dimension(:)       :: m

!--------------------------------------------------------------------------
interface

   subroutine ginv(g,detg,mask,gup)
      implicit none
      real*8,dimension(:,:,:,:),target :: g,gup
      real*8,dimension(:,:,:)          :: detg,mask
      real*8,parameter                 :: zero = 0.0d0
      real*8,parameter                 :: value  = 1.0d0
   end subroutine
   subroutine gamma(g,gup,h,t,nl,m,d,c,ca)
      implicit none
      real*8                                        :: h
      real*8,dimension(:,:,:,:),target,intent(in)   :: g,gup
      real*8,dimension(:,:,:,:),target,intent(out)  :: c,t,ca
      integer                                       :: nl
      real*8,dimension(:)                        :: d
      integer,dimension(:)                       :: m
   end subroutine
   subroutine ricci(g,gup,c,ca,h,t,rieman,m,d,nl,r)
      implicit none
      integer                              :: nl
      real*8                               :: h
      real*8,dimension(:,:,:,:),target     :: g,gup,c,t,r
      real*8,dimension(:,:,:,:),target     :: rieman,ca
      real*8,dimension(:  )                :: d
      integer,dimension(:)                 :: m
   end subroutine
   subroutine k_dot(g,gup,k,alpha,beta,r,c,h,t,m,d,nl,kdot)
      implicit none
      integer                              :: nl
      real*8                               :: h
      real*8,dimension(:,:,:,:),target     :: g,gup,k,beta,r,c,t,kdot
      real*8,dimension(:,:,:)              :: alpha
      real*8,dimension(  : )               :: d
      integer,dimension( :  )              :: m
   end subroutine
   subroutine g_dot(g,k,alpha,beta,c,gup,h,m,d,nl,gdot)
      implicit none
      integer                              :: nl
      real*8                               :: h
      real*8,dimension(:,:,:,:)            :: g,gup,k,c,beta,gdot
      real*8,dimension(:,:,:)              :: alpha
      real*8,dimension(:)                  :: d
      integer,dimension(:)                 :: m
   end subroutine
   subroutine leap(sn,snm1,dt,rhs,snp1)
      implicit none
      real*8,dimension(:,:,:,:),target     :: sn,snm1,snp1,rhs
      real*8                               :: dt
   end subroutine

   real*8 function n2i(n,a,mask,m,nl)
      implicit none
      integer                       :: n,nl
      real*8,dimension(:,:,:)       :: a,mask
      integer                       :: i,j,k,l
      integer,dimension(:)       :: m
	end function


end interface
!--------------------------------------------------------------------------


   call ginv(g,t(1,:,:,:),mask,gup)
   call gamma(g,gup,h,t,nl,m,d,c,ca) 
   call ricci(g,gup,c,ca,h,t,rieman,m,d,nl,r)

   write(0,*)'****gij norms****'
write(6,*) 'before call'
write(6,*) g(1,3,3,3:5)
   xns  = n2i(17,g(1,:,:,:),mask,m,nl)
   write(0,*) 'g11 ',xns
write(6,*) 'before call'
write(6,*) g(2,3,3,3:5)
   xns  = n2i(17,g(2,:,:,:),mask,m,nl)
   write(0,*)  'g12 ',xns
write(6,*) 'before call'
write(6,*) g(3,3,3,3:5)
   xns  = n2i(17,g(3,:,:,:),mask,m,nl)
   write(0,*)  'g13 ',xns
write(6,*) 'before call'
write(6,*) g(4,3,3,3:5)
   xns  = n2i(17,g(4,:,:,:),mask,m,nl)
   write(0,*)  'g22 ',xns
write(6,*) 'before call'
write(6,*) g(5,3,3,3:5)
   xns  = n2i(17,g(5,:,:,:),mask,m,nl)
   write(0,*)  'g23 ',xns
write(6,*) 'before call'
write(6,*) g(6,3,3,3:5)
   xns  = n2i(17,g(6,:,:,:),mask,m,nl)
   write(0,*)  'g33 ',xns
stop
   write(0,*)'****c1ij norms****'
   xns  = n2i(17,c(1,:,:,:),mask,m,nl)
   write(0,*) 'c111 ',xns
   xns  = n2i(17,c(2,:,:,:),mask,m,nl)
   write(0,*)  'c112 ',xns
   xns  = n2i(17,c(3,:,:,:),mask,m,nl)
   write(0,*)  'c113 ',xns
   xns  = n2i(17,c(4,:,:,:),mask,m,nl)
   write(0,*)  'c122 ',xns
   xns  = n2i(17,c(5,:,:,:),mask,m,nl)
   write(0,*)  'c123 ',xns
   xns  = n2i(17,c(6,:,:,:),mask,m,nl)
   write(0,*)  'c133 ',xns

   write(0,*)'****rij norms****'
   xns  = n2i(17,r(1,:,:,:),mask,m,nl)
   write(0,*) 'r11 ',xns
   xns  = n2i(17,r(2,:,:,:),mask,m,nl)
   write(0,*)  'r12 ',xns
   xns  = n2i(17,r(3,:,:,:),mask,m,nl)
   write(0,*)  'r13 ',xns
   xns  = n2i(17,r(4,:,:,:),mask,m,nl)
   write(0,*)  'r22 ',xns
   xns  = n2i(17,r(5,:,:,:),mask,m,nl)
   write(0,*)  'r23 ',xns
   xns  = n2i(17,r(6,:,:,:),mask,m,nl)
   write(0,*)  'r33 ',xns

 
   call k_dot(g,gup,k,alpha,beta,r,c,h,t,m,d,nl,t)
   call leap(k,km1,dt,t,kp1)

   write(0,*)'****kij norms****'
   xns  = n2i(17,t(1,:,:,:),mask,m,nl)
   write(0,*) 't11 ',xns
   xns  = n2i(17,t(2,:,:,:),mask,m,nl)
   write(0,*) 't12 ',xns
   xns  = n2i(17,t(3,:,:,:),mask,m,nl)
   write(0,*) 't13 ',xns
   xns  = n2i(17,t(4,:,:,:),mask,m,nl)
   write(0,*) 't22 ',xns
   xns  = n2i(17,t(5,:,:,:),mask,m,nl)
   write(0,*) 't23 ',xns
   xns  = n2i(17,t(6,:,:,:),mask,m,nl)
   write(0,*) 't33 ',xns


   call g_dot(g,k,alpha,beta,c,gup,h,m,d,nl,t)
	write(6,*) 'calling leap'
   call leap(g,gm1,dt,t,gp1)

   write(0,*)'****gij norms****'
   xns  = n2i(17,t(1,:,:,:),mask,m,nl)
   write(0,*) 't11 ',xns
   xns  = n2i(17,t(2,:,:,:),mask,m,nl)
   write(0,*)  't12 ',xns
   xns  = n2i(17,t(3,:,:,:),mask,m,nl)
   write(0,*)  't13 ',xns
   xns  = n2i(17,t(4,:,:,:),mask,m,nl)
   write(0,*)  't22 ',xns
   xns  = n2i(17,t(5,:,:,:),mask,m,nl)
   write(0,*)  't23 ',xns
   xns  = n2i(17,t(6,:,:,:),mask,m,nl)
   write(0,*)  't33 ',xns


return 
end 
!----------------------------------------------------------------------- 
