! $Header$

#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Functions.h"
#include "cctk_Parameters.h"

subroutine CalcK (CCTK_ARGUMENTS)
  implicit none
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_FUNCTIONS
  DECLARE_CCTK_PARAMETERS
  
  CCTK_INT, parameter :: idummy=0
  integer,  parameter :: ik=kind(idummy)
  
  CCTK_INT, parameter :: bndwidth = 1
  
  CCTK_INT        :: len_boundary, len_options
  character(1000) :: str_boundary, str_options
  integer         :: options
  
  CCTK_REAL :: dt, dx(3)
  
  CCTK_REAL :: gama(3,3), gama_dot(3,3), dgama(3,3,3)
  CCTK_REAL :: alfa
  CCTK_REAL :: beta(3), dbeta(3,3)
  CCTK_REAL :: kk(3,3)
  
  integer :: imin(3), imax(3)
  integer :: i, j, k
  integer :: a, b, c
  integer :: ierr
  
  dt = CCTK_DELTA_TIME
  dx(:) = CCTK_DELTA_SPACE(:)
  
  call CCTK_FortranString &
       (len_boundary, int(extcurv_boundary,ik), str_boundary)
  call CCTK_FortranString &
       (len_options, int(extcurv_boundary_options,ik), str_options)
  call Util_TableCreateFromString (options, str_options)
  if (options<0) then
     call CCTK_WARN (0, "Parameter ""extcurv_boundary_options"" has an illegal syntax")
  end if
  
  imin(:) = 1+cctk_nghostzones(:)
  imax(:) = cctk_lsh(:)-cctk_nghostzones(:)
  where (cctk_bbox(1::2)/=0) imin(:) = 1+bndwidth
  where (cctk_bbox(2::2)/=0) imax(:) = cctk_lsh(:)-bndwidth
  
  do k = imin(3), imax(3)
     do j = imin(2), imax(2)
        do i = imin(1), imax(1)
           
           gama(1,1) = gxx(i,j,k)
           gama(1,2) = gxy(i,j,k)
           gama(1,3) = gxz(i,j,k)
           gama(2,2) = gyy(i,j,k)
           gama(2,3) = gyz(i,j,k)
           gama(3,3) = gzz(i,j,k)
           gama(2,1) = gama(1,2)
           gama(3,1) = gama(1,3)
           gama(3,2) = gama(2,3)
           
#if 0
           gama_dot(1,1) = (-3*gxx(i,j,k) + 4*gxx_prev(i,j,k) - gxx_prev2(i,j,k)) / (2*dt)
           gama_dot(1,2) = (-3*gxy(i,j,k) + 4*gxy_prev(i,j,k) - gxy_prev2(i,j,k)) / (2*dt)
           gama_dot(1,3) = (-3*gxz(i,j,k) + 4*gxz_prev(i,j,k) - gxz_prev2(i,j,k)) / (2*dt)
           gama_dot(2,2) = (-3*gyy(i,j,k) + 4*gyy_prev(i,j,k) - gyy_prev2(i,j,k)) / (2*dt)
           gama_dot(2,3) = (-3*gyz(i,j,k) + 4*gyz_prev(i,j,k) - gyz_prev2(i,j,k)) / (2*dt)
           gama_dot(3,3) = (-3*gzz(i,j,k) + 4*gzz_prev(i,j,k) - gzz_prev2(i,j,k)) / (2*dt)
#else
           gama_dot(1,1) = (gxx_next(i,j,k) - gxx_prev(i,j,k)) / (2*dt)
           gama_dot(1,2) = (gxy_next(i,j,k) - gxy_prev(i,j,k)) / (2*dt)
           gama_dot(1,3) = (gxz_next(i,j,k) - gxz_prev(i,j,k)) / (2*dt)
           gama_dot(2,2) = (gyy_next(i,j,k) - gyy_prev(i,j,k)) / (2*dt)
           gama_dot(2,3) = (gyz_next(i,j,k) - gyz_prev(i,j,k)) / (2*dt)
           gama_dot(3,3) = (gzz_next(i,j,k) - gzz_prev(i,j,k)) / (2*dt)
#endif
           gama_dot(2,1) = gama_dot(1,2)
           gama_dot(3,1) = gama_dot(1,3)
           gama_dot(3,2) = gama_dot(2,3)
           
           dgama(1,1,1) = (gxx(i+1,j,k) - gxx(i-1,j,k)) / (2*dx(1))
           dgama(1,2,1) = (gxy(i+1,j,k) - gxy(i-1,j,k)) / (2*dx(1))
           dgama(1,3,1) = (gxz(i+1,j,k) - gxz(i-1,j,k)) / (2*dx(1))
           dgama(2,2,1) = (gyy(i+1,j,k) - gyy(i-1,j,k)) / (2*dx(1))
           dgama(2,3,1) = (gyz(i+1,j,k) - gyz(i-1,j,k)) / (2*dx(1))
           dgama(3,3,1) = (gzz(i+1,j,k) - gzz(i-1,j,k)) / (2*dx(1))
           dgama(1,1,2) = (gxx(i,j+1,k) - gxx(i,j-1,k)) / (2*dx(2))
           dgama(1,2,2) = (gxy(i,j+1,k) - gxy(i,j-1,k)) / (2*dx(2))
           dgama(1,3,2) = (gxz(i,j+1,k) - gxz(i,j-1,k)) / (2*dx(2))
           dgama(2,2,2) = (gyy(i,j+1,k) - gyy(i,j-1,k)) / (2*dx(2))
           dgama(2,3,2) = (gyz(i,j+1,k) - gyz(i,j-1,k)) / (2*dx(2))
           dgama(3,3,2) = (gzz(i,j+1,k) - gzz(i,j-1,k)) / (2*dx(2))
           dgama(1,1,3) = (gxx(i,j,k+1) - gxx(i,j,k-1)) / (2*dx(3))
           dgama(1,2,3) = (gxy(i,j,k+1) - gxy(i,j,k-1)) / (2*dx(3))
           dgama(1,3,3) = (gxz(i,j,k+1) - gxz(i,j,k-1)) / (2*dx(3))
           dgama(2,2,3) = (gyy(i,j,k+1) - gyy(i,j,k-1)) / (2*dx(3))
           dgama(2,3,3) = (gyz(i,j,k+1) - gyz(i,j,k-1)) / (2*dx(3))
           dgama(3,3,3) = (gzz(i,j,k+1) - gzz(i,j,k-1)) / (2*dx(3))
           dgama(2,1,:) = dgama(1,2,:)
           dgama(3,1,:) = dgama(1,3,:)
           dgama(3,2,:) = dgama(2,3,:)
           
           alfa = alp(i,j,k)
           
           beta(1) = betax(i,j,k)
           beta(2) = betay(i,j,k)
           beta(3) = betaz(i,j,k)
           
           dbeta(1,1) = (betax(i+1,j,k) - betax(i-1,j,k)) / (2*dx(1))
           dbeta(2,1) = (betay(i+1,j,k) - betay(i-1,j,k)) / (2*dx(1))
           dbeta(3,1) = (betaz(i+1,j,k) - betaz(i-1,j,k)) / (2*dx(1))
           dbeta(1,2) = (betax(i,j+1,k) - betax(i,j-1,k)) / (2*dx(2))
           dbeta(2,2) = (betay(i,j+1,k) - betay(i,j-1,k)) / (2*dx(2))
           dbeta(3,2) = (betaz(i,j+1,k) - betaz(i,j-1,k)) / (2*dx(2))
           dbeta(1,3) = (betax(i,j,k+1) - betax(i,j,k-1)) / (2*dx(3))
           dbeta(2,3) = (betay(i,j,k+1) - betay(i,j,k-1)) / (2*dx(3))
           dbeta(3,3) = (betaz(i,j,k+1) - betaz(i,j,k-1)) / (2*dx(3))
           
           ! d/dt gamma_ij = - 2 alpha K_ij
           !                 + gamma_kj d_i beta^k
           !                 + gamma_kj d_i beta^k
           !                 + beta^k d_k gamma_ij
           
           do a=1,3
              do b=1,3
                 kk(a,b) = - gama_dot(a,b)
                 do c=1,3
                    kk(a,b) = kk(a,b) + gama(k,j) * dbeta(k,i) &
                         &            + gama(i,k) * dbeta(k,j) &
                         &            + beta(k) * dgama(i,j,k)
                 end do
                 kk(a,b) = kk(a,b) / (2*alfa)
              end do
           end do
           
           kxx(i,j,k) = kk(1,1)
           kxy(i,j,k) = kk(1,2)
           kxz(i,j,k) = kk(1,3)
           kyy(i,j,k) = kk(2,2)
           kyz(i,j,k) = kk(2,3)
           kzz(i,j,k) = kk(3,3)
           
        end do
     end do
  end do
  
  call CartSymGN (ierr, cctkGH, "ADMBase::curv")
  if (ierr /= 0) call CCTK_WARN (0, "internal error")
  ierr = Boundary_SelectGroupForBC (cctkGH, CCTK_ALL_FACES, bndwidth, &
       int(options,ik), "ADMBase::curv", str_boundary)
  if (ierr /= 0) call CCTK_WARN (0, "internal error")
  
end subroutine CalcK