/*@@
  @file      StaggeredLeapfrog2.F
  @desc
     Second step of Staggered leapfrog stepper for the ADM Evolution.
  @enddesc
@@*/

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

#include "CactusEinstein/Einstein/src/Einstein.h"

/*@@
  @routine    ADM_StaggeredLeapfrog2
  @desc
     Staggered leapfrog stepper for the ADM Evolution.
  @enddesc
  @calls
  @calledby
  @history

  @endhistory

@@*/

c     ==================================================================
c
c     STAGGERED LEAPFROG:
c
c     metric and lapse on the n (full) timesteps
c     extrinsic curvature on the n+1/2 (half) timesteps
c     the GF kxx,... are interpolated onto the FULL timesteps
c     but not used in the evolution scheme
c
c     all in GFList declared:
c     g_p: n-1   ADM_h_stag:n+1/2
c     g  : n     ADM_h_p :n-1/2
c                h   :n   <-- this is shown to the outside, has same
c                            time slice as g and alp
c
c     alp_p: n-1
c     alp  : n
c
c     ==================================================================


      subroutine ADM_StaggeredLeapfrog2(CCTK_ARGUMENTS)

      USE ADM_Scalars

      implicit none

      DECLARE_CCTK_ARGUMENTS
      DECLARE_CCTK_PARAMETERS
      DECLARE_CCTK_FUNCTIONS

      integer :: i,j,k
      integer :: nx,ny,nz
      integer :: ierror

      CCTK_REAL :: dx,dy,dz,dt

      CCTK_REAL :: 
     &     adm_rhsK_metric_xx,adm_rhsK_metric_xy,adm_rhsK_metric_xz,
     &     adm_rhsK_metric_yy,adm_rhsK_metric_yz,adm_rhsK_metric_zz,
     &     adm_rhsK_matter_xx,adm_rhsK_matter_xy,adm_rhsK_matter_xz,
     &     adm_rhsK_matter_yy,adm_rhsK_matter_yz,adm_rhsK_matter_zz

c     Declarations for macros from ADM Evolution
#include "macro/KSOURCES_declare.h"

      CCTK_REAL zero,one,pi

c     ==================================================================

      zero = 0.0D0
      one  = 1.0D0
      pi = acos(-one)

      dt = cctk_delta_time
      dx = cctk_delta_space(1)
      dy = cctk_delta_space(2)
      dz = cctk_delta_space(3)

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


c     ==================================================================
c
c     Evolve extrinsic curvature 
c
c     ==================================================================

c     2c: Everything has the right names to find the source of the K 
c         evolution equation using the standard macros

      do k=2,nz-1
         do j=2,ny-1
            do i=2,nx-1
#include "macro/KSOURCES_guts.h"
#include "macro/KSOURCES_undefine.h"
            end do
         end do
      end do

c     2d: Relabel the new old timeslice (n+1/2)
c     2e: Evolve the extrinsic curvature. The result is saved
c     in the variables with no suffix in order for the boundary
c     conditions to work.

      do k=1,nz
         do j=1,ny
            do i=1,nx
               ADM_kxx_p(i,j,k) = ADM_kxx_stag(i,j,k)
               kxx(i,j,k)       = ADM_kxx_p(i,j,k) + adms_kxx(i,j,k)*dt
            end do
         end do
      end do

      do k=1,nz
         do j=1,ny
            do i=1,nx
               ADM_kxy_p(i,j,k) = ADM_kxy_stag(i,j,k)
               kxy(i,j,k)       = ADM_kxy_p(i,j,k) + adms_kxy(i,j,k)*dt
            end do
         end do
      end do

      do k=1,nz
         do j=1,ny
            do i=1,nx
               ADM_kxz_p(i,j,k) = ADM_kxz_stag(i,j,k)
               kxz(i,j,k)       = ADM_kxz_p(i,j,k) + adms_kxz(i,j,k)*dt
            end do
         end do
      end do

      do k=1,nz
         do j=1,ny
            do i=1,nx
               ADM_kyy_p(i,j,k) = ADM_kyy_stag(i,j,k)
               kyy(i,j,k)       = ADM_kyy_p(i,j,k) + adms_kyy(i,j,k)*dt
            end do
         end do
      end do

      do k=1,nz
         do j=1,ny
            do i=1,nx
               ADM_kyz_p(i,j,k) = ADM_kyz_stag(i,j,k)
               kyz(i,j,k)       = ADM_kyz_p(i,j,k) + adms_kyz(i,j,k)*dt
            end do
         end do
      end do

      do k=1,nz
         do j=1,ny
            do i=1,nx
               ADM_kzz_p(i,j,k) = ADM_kzz_stag(i,j,k)
               kzz(i,j,k)       = ADM_kzz_p(i,j,k) + adms_kzz(i,j,k)*dt
            end do
         end do
      end do

      if (minimal_communications.ne.1) then
         call CCTK_SyncGroup(ierror,cctkGH,"einstein::curv")
      end if
      call ADM_Boundary(cctkGH,"einstein::curv")

c     Rename variables.
c     2f: Use the extrinsic curvature at n+3/2 and n+1/2 to 
c         interpolate to n+1

      do k=1,nz
         do j=1,ny
            do i=1,nx
                ADM_kxx_stag(i,j,k) = kxx(i,j,k)
                kxx(i,j,k) = 0.5D0*(ADM_kxx_p(i,j,k) + ADM_kxx_stag(i,j,k))
            end do
         end do
      end do

      do k=1,nz
         do j=1,ny
            do i=1,nx
                ADM_kxy_stag(i,j,k) = kxy(i,j,k)
                kxy(i,j,k) = 0.5D0*(ADM_kxy_p(i,j,k) + ADM_kxy_stag(i,j,k))
            end do
         end do
      end do

      do k=1,nz
         do j=1,ny
            do i=1,nx
                ADM_kxz_stag(i,j,k) = kxz(i,j,k)
                kxz(i,j,k) = 0.5D0*(ADM_kxz_p(i,j,k) + ADM_kxz_stag(i,j,k))
            end do
         end do
      end do

      do k=1,nz
         do j=1,ny
            do i=1,nx
                ADM_kyy_stag(i,j,k) = kyy(i,j,k)
                kyy(i,j,k) = 0.5D0*(ADM_kyy_p(i,j,k) + ADM_kyy_stag(i,j,k))
            end do
         end do
      end do

      do k=1,nz
         do j=1,ny
            do i=1,nx
                ADM_kyz_stag(i,j,k) = kyz(i,j,k)
                kyz(i,j,k) = 0.5D0*(ADM_kyz_p(i,j,k) + ADM_kyz_stag(i,j,k))
            end do
         end do
      end do

      do k=1,nz
         do j=1,ny
            do i=1,nx
                ADM_kzz_stag(i,j,k) = kzz(i,j,k)
                kzz(i,j,k) = 0.5D0*(ADM_kzz_p(i,j,k) + ADM_kzz_stag(i,j,k))
            end do
         end do
      end do

      end subroutine ADM_StaggeredLeapfrog2