!/*@@
!
! $Id: evolve_k.f90,v 1.2 1997/02/10 16:38:09 scott Exp $
!
!@@*/

!/*@@
!   @routine k_dot
!   @author R.Correl M.Huq S.Klasky R.Matzner
!   @desc
!   driver for computing k_dot
!   @enddesc
!   @calledby tstep
!   @calls diff
!   @calls k_dot0
!   @calls k_dot1
!   @calls k_dot2
!   @calls k_dot3
!   @comment
!   uses modules tensors and diffinterface
!   @endcomment
!@@*/
!-----------------------------------------------------------------------
! evolve_k.f90
!-----------------------------------------------------------------------
!  g_     := lower metric 
!  gup_   := upper metric
!  k_     := extrinsic curvature
!  alpha  := lapse function
!  beta_  := shift vector
!  c_     := connections
!  r_     := ricci
!  t_     := temporaries
!

subroutine k_dot(g_,params,k_,alpha,beta_,c_,gup_,r_,t_,h,m,d,nl,kdot_)
use input_params
use tensors
real*8,intent(in)                                  :: h
type(GridFunction),dimension(3,3),intent(in)       :: g_,gup_,k_, r_
type(GridFunction),dimension(3,3),intent(inout)    :: kdot_
type(GridFunction),dimension(3,3,3),intent(in)     :: c_
type(GridFunction),dimension(3),intent(in)         :: beta_
type(GridFunction),dimension(3),intent(inout)      :: t_
real*8, dimension(:,:,:), intent(in)               :: alpha
integer                              :: nl
real*8,dimension(:)                  :: d,params
integer,dimension(:)                 :: m


integer :: i,j

INTERFACE
subroutine k_dot0(params,g_,k_,alpha,beta_,c_,gup_,r_,t_,h,m,d,nl,kdot_)
use input_params
  use tensors
  use diffinterface
  real*8,intent(in)                                  :: h
  type(GridFunction),dimension(3,3),intent(in)       :: g_,gup_,k_, r_
  type(GridFunction),dimension(3,3),intent(inout)    :: kdot_
  type(GridFunction),dimension(3,3,3),intent(in)     :: c_
  type(GridFunction),dimension(3),intent(in)         :: beta_
  type(GridFunction),dimension(3),intent(inout)      :: t_
  real*8,dimension(:,:,:),intent(in)                 :: alpha
  integer                              :: nl
  real*8,dimension(:)                  :: d,params
  integer,dimension(:)                 :: m
  end subroutine k_dot0

subroutine k_dot1(params,g_,k_,alpha,beta_,c_,gup_,r_,t_,h,m,d,nl,kdot_)
use input_params
  use tensors
  use diffinterface
  real*8,intent(in)                                  :: h
  type(GridFunction),dimension(3,3),intent(in)       :: g_,gup_,k_, r_
  type(GridFunction),dimension(3,3),intent(inout)    :: kdot_
  type(GridFunction),dimension(3,3,3),intent(in)     :: c_
  type(GridFunction),dimension(3),intent(in)         :: beta_
  type(GridFunction),dimension(3),intent(inout)      :: t_
  real*8,dimension(:,:,:),intent(in)                 :: alpha
  integer                              :: nl
  real*8,dimension(:)                  :: d,params
  integer,dimension(:)                 :: m
  end subroutine k_dot1

subroutine k_dot2(params,g_,k_,alpha,beta_,c_,gup_,r_,t_,h,m,d,nl,kdot_)
use input_params
  use tensors
  use diffinterface
  real*8,intent(in)                                  :: h
  type(GridFunction),dimension(3,3),intent(in)       :: g_,gup_,k_, r_
  type(GridFunction),dimension(3,3),intent(inout)    :: kdot_
  type(GridFunction),dimension(3,3,3),intent(in)     :: c_
  type(GridFunction),dimension(3),intent(in)         :: beta_
  type(GridFunction),dimension(3),intent(inout)      :: t_
  real*8,dimension(:,:,:),intent(in)                 :: alpha
  integer                              :: nl
  real*8,dimension(:)                  :: d,params
  integer,dimension(:)                 :: m
  end subroutine k_dot2

subroutine k_dot3(params,g_,k_,alpha,beta_,c_,gup_,r_,t_,h,m,d,nl,kdot_)
use input_params
  use tensors
  use diffinterface
  implicit none
  real*8,intent(in)                                  :: h
  type(GridFunction),dimension(3,3),intent(in)       :: g_,gup_,k_, r_
  type(GridFunction),dimension(3,3),intent(inout)    :: kdot_
  type(GridFunction),dimension(3,3,3),intent(in)     :: c_
  type(GridFunction),dimension(3),intent(in)         :: beta_
  type(GridFunction),dimension(3),intent(inout)      :: t_
  real*8,dimension(:,:,:),intent(in)                 :: alpha
  integer                              :: nl
  real*8,dimension(:)                  :: d,params
  integer,dimension(:)                 :: m

  end subroutine k_dot3

END INTERFACE


   do i = 1, 3
     do j = i, 3
       kdot_(i,j)%d = 0.0d0
     end do
   end do

!   Break down evolution equation for K_{ij} term by term 

!   Partial derivatives of the Lapse
   call k_dot0(params,g_,k_,alpha,beta_,c_,gup_,r_,t_,h,m,d,nl,kdot_)

!      DAGH: Sync 
   call k_dot1(params,g_,k_,alpha,beta_,c_,gup_,r_,t_,h,m,d,nl,kdot_)

!   Second partials and addition of connections => -D_i D_j \alpha
  call k_dot2(params,g_,k_,alpha,beta_,c_,gup_,r_,t_,h,m,d,nl,kdot_)

!      DAGH: Sync

!  alpha*(R_ab + K*K_ab -2 K_acK^c_b) terms
   call k_dot3(params,g_,k_,alpha,beta_,c_,gup_,r_,t_,h,m,d,nl,kdot_)

return
end

!/*@@
!   @routine k_dot0
!   @author R.Correl M.Huq S.Klasky R.Matzner
!   @desc
!   ***********************************************************************
!    k_dot0: \partial_j \alpha
!   Broken up this way for use with DAGH.
!   ***********************************************************************
!   @enddesc
!   @calledby k_dot
!   @calls diff
!   @comment
!   uses modules tensors and diffinterface
!   @endcomment
!@@*/

subroutine k_dot0(params,g_,k_,alpha,beta_,c_,gup_,r_,t_,h,m,d,nl,kdot_)
use input_params
  use tensors
  use diffinterface
  implicit none
  real*8,intent(in)                                  :: h

  type(GridFunction),dimension(3,3),intent(in)       :: g_,gup_,k_, r_
  type(GridFunction),dimension(3,3),intent(inout)    :: kdot_
  type(GridFunction),dimension(3,3,3),intent(in)     :: c_
  type(GridFunction),dimension(3),intent(in)         :: beta_
  type(GridFunction),dimension(3),intent(inout)      :: t_
  real*8,dimension(:,:,:),intent(in)                 :: alpha
  real*8                                             :: dummy
  integer                              :: nl
  real*8,dimension(:)                  :: d,params
  integer,dimension(:)                 :: m


! Local variables

  integer   :: i, j, x, y, z

  x = 1
  y = 2
  z = 3

! Derivative of lapse

  call diff(params,alpha,alpha,dummy,x,h,int(params(P_BNDRY_DERIV_ORDER)),3,m,d,nl,t_(x)%d)
  call diff(params,alpha,alpha,dummy,y,h,int(params(P_BNDRY_DERIV_ORDER)),3,m,d,nl,t_(y)%d)
  call diff(params,alpha,alpha,dummy,z,h,int(params(P_BNDRY_DERIV_ORDER)),3,m,d,nl,t_(z)%d)

  return
end
!/*@@
!   @routine k_dot1
!   @author R.Correl M.Huq S.Klasky R.Matzner
!   @desc
!   ***********************************************************************
!   k_dot1: -D_i D_j \alpha 
!   ***********************************************************************
!   @enddesc
!   @calledby k_dot
!   @calls diff
!   @comment
!   uses modules tensors and diffinterface
!   @endcomment
!@@*/

subroutine k_dot1(params,g_,k_,alpha,beta_,c_,gup_,r_,t_,h,m,d,nl,kdot_)
use input_params
  
use tensors
  use diffinterface
  implicit none

  real*8,intent(in)                                  :: h
  type(GridFunction),dimension(3,3),intent(in)       :: g_,gup_,k_, r_
  type(GridFunction),dimension(3,3),intent(inout)    :: kdot_
  type(GridFunction),dimension(3,3,3),intent(in)     :: c_
  type(GridFunction),dimension(3),intent(in)         :: beta_
  type(GridFunction),dimension(3),intent(inout)      :: t_
  real*8,dimension(:,:,:),intent(in)                 :: alpha
  integer                              :: nl
  real*8,dimension(:)                  :: d,params
  integer,dimension(:)                 :: m

  real*8                                             :: dummy
  integer,dimension(3,3)                             :: wd

! Local variables

  integer   :: i, j, k, x, y, z

  x = 1
  y = 2
  z = 3

! Second partials of \partial_j \alpha: store in kdot_
! \partial_j \alpha stored in t_ to begin with.
  wd(1,1) = 4
  wd(1,2) = 7
  wd(1,3) = 8
  wd(2,1) = 7
  wd(2,2) = 5
  wd(2,3) = 9
  wd(3,1) = 8
  wd(3,2) = 9
  wd(3,3) = 6

! Second partials of \partial_j \alpha: store in kdot_
! \partial_j \alpha stored in t_ to begin with.
   do i = 1 ,3
      do j = 1 ,3
        if (i <= j) then   ! Take symmetry into account
         call diff(params,alpha,alpha,dummy,wd(i,j),h,int(params(P_BNDRY_DERIV_ORDER)),3, &
                   m,d,nl,kdot_(i,j)%d)
        end if
      end do
   end do


! Form the covariant derivative.
! Multiply by -1 to get -D_i D_j \alpha : Since this is the first term
! in kdot_ it is fine to do it this way. If this is moved to later then
! the -1 multiplication must be done earlier.
  
  do i = 1 , 3
    do j = 1 , 3
      if (i <= j) then   ! Take symmetry into account

        kdot_(i,j)%d = -kdot_(i,j)%d + t_(1)%d * c_(1,i,j)%d + &
                                     t_(2)%d * c_(2,i,j)%d + &
                                     t_(3)%d * c_(3,i,j)%d 
      end if
    end do
  end do


  return
end

!/*@@
!   @routine k_dot2
!   @author R.Correl M.Huq S.Klasky R.Matzner
!   @desc
!   ***********************************************************************
!   k_dot2: Lie derivative terms
!   K_kj \partial_i \beta^k + K_ik \partial_j \beta^k
!   Note that the convective part beta^k \partial_k K_ij is missing because
!   we use causal differencing.
!   ***********************************************************************
!   @enddesc
!   @calledby k_dot
!   @calls diff
!   @comment
!   uses modules tensors and diffinterface
!   @endcomment
!@@*/

subroutine k_dot2(params,g_,k_,alpha,beta_,c_,gup_,r_,t_,h,m,d,nl,kdot_)
use input_params
  use tensors
  use diffinterface
  implicit none

  real*8,intent(in)                    :: h
  type(GridFunction),dimension(3,3),intent(in)       :: g_,gup_,k_, r_
  type(GridFunction),dimension(3,3),intent(inout)    :: kdot_
  type(GridFunction),dimension(3,3,3),intent(in)     :: c_
  type(GridFunction),dimension(3),intent(in)         :: beta_
  type(GridFunction),dimension(3),intent(inout)      :: t_
  real*8,dimension(:,:,:),intent(in)                 :: alpha
  real*8                                             :: dummy
  integer                              :: nl
  real*8,dimension(:)                  :: d,params
  integer,dimension(:)                 :: m


! Local variables

  integer   :: i, j, k, x, y, z

  x = 1
  y = 2
  z = 3

! Next K_kj * \partial_i \beta^k

   do i = 1 ,3
      do j = 1 ,3
        if (i <= j) then   ! Take symmetry into account

          do k = 1, 3

            ! Derivative in the i direction of each k and store in t_
            call diff(params,beta_(k)%d,alpha,dummy,i,h,int(params(P_BNDRY_DERIV_ORDER)),3, &
                      m,d,nl,t_(i)%d)  ! beta^k,i

            ! Add K_kj * \partial_i \beta^k combination: 
            !           t_(i) contains \partial_i \beta^k for current k
            kdot_(i,j)%d = kdot_(i,j)%d + k_(k,j)%d * t_(i)%d

          end do

        end if 
      end do
   end do

! Next K_ik * \partial_j \beta^k

   do i = 1 ,3
      do j = 1 ,3
        if (i <= j) then   ! Take symmetry into account

          do k = 1, 3

            ! Derivative in the j direction of each k and store in t_
            call diff(params,beta_(k)%d,alpha,dummy,j,h,int(params(P_BNDRY_DERIV_ORDER)),3, &
                      m,d,nl,t_(j)%d) ! beta^k,j

            ! Add K_ik * \partial_j \beta^k combination: 
            !           t_(j) contains \partial_i \beta^k for current k
            kdot_(i,j)%d = kdot_(i,j)%d + k_(i,k)%d * t_(j)%d

          end do

        end if 
      end do
   end do

   return
end

!/*@@
!   @routine k_dot3
!   @author R.Correl M.Huq S.Klasky R.Matzner
!   @desc
!   ***********************************************************************
!   k_dot3: alpha*(R_ab + K*K_ab -2 K_acK^c_b) terms
!   ***********************************************************************
!   @enddesc
!   @calledby k_dot
!   @comment
!   uses modules tensors and diffinterface
!   @endcomment
!@@*/

subroutine k_dot3(params,g_,k_,alpha,beta_,c_,gup_,r_,t_,h,m,d,nl,kdot_)
use input_params
  use tensors
  use diffinterface
  implicit none

  real*8,intent(in)                                  :: h
  type(GridFunction),dimension(3,3),intent(in)       :: g_,gup_,k_, r_
  type(GridFunction),dimension(3,3),intent(inout)    :: kdot_
  type(GridFunction),dimension(3,3,3),intent(in)     :: c_
  type(GridFunction),dimension(3),intent(in)         :: beta_
  type(GridFunction),dimension(3),intent(inout)      :: t_
  real*8,dimension(:,:,:),intent(in)                 :: alpha
  integer                              :: nl
  real*8,dimension(:)                  :: d,params
  integer,dimension(:)                 :: m



! Local variables

  integer   :: i, j, k, l


!  \alpha*R_ij
   do i = 1 ,3
      do j = i ,3
          kdot_(i,j)%d =  kdot_(i,j)%d + alpha * r_(i,j)%d
      end do
   end do


!  \alpha*K*K_ij

   ! Use t_(1) to store trace
 
   t_(1)%d = 0.0d0
   do i = 1, 3
     do j = 1 ,3
       t_(1)%d = t_(1)%d  + gup_(i,j)%d * k_(i,j)%d
     end do
   end do

   do i = 1 ,3
      do j = i ,3
          kdot_(i,j)%d =  kdot_(i,j)%d + alpha * t_(1)%d * k_(i,j)%d
      end do
   end do

! -2 k_il k_mj g^{lm}
   do i = 1 ,3
      do j = i ,3
	 ! For current j get k^l_j and store in t_
	 do k = 1, 3
	    do l = 1, 3
	       kdot_(i,j)%d = kdot_(i,j)%d - 2.0d0 * alpha * k_(i,k)%d * k_(j,l)%d *gup_(k,l)%d
	    end do
	 end do
      end do
   end do


  return
end
!***********************************************************************
