/*@@
  @file      StaggeredLeapfrog1.F
  @desc
     First step of staggered leapfrog stepper for the ADM Evolution.
  @enddesc
@@*/

c --------------------------------------------------------------

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

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

c --------------------------------------------------------------


/*@@
  @routine    ADM_StaggeredLeapfrog1
  @desc
     First step of 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     ==================================================================
c     ==================================================================
c
c     STAGGERED LEAPFROG initialization using predictor corrector
c     method :
c
c     GHiteration = 1: The SL routine calls ADM_predcorr to
c     do a simple predictor corrector 2nd order accurate evolution step
c
c     GHiteration = 2: the SL routine calls ADM_predcorr to
c     do a second predictor corrector step.  The resulting 3-tiered
c     extrinsic curvature is used to set up staggered data
c     for the GHiteration = 2 SL step.  Whew!
c     ==================================================================

      subroutine ADM_StaggeredLeapfrog1(CCTK_ARGUMENTS)

      USE ADM_Scalars

      implicit none

      DECLARE_CCTK_ARGUMENTS
      DECLARE_CCTK_PARAMETERS

      integer CCTK_Equals
      integer :: i,j,k,stencil_size_y
      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    Temporaries for 3-tiered extrinsic curvature juggling
      CCTK_REAL ::
     &     kxx_t1,kxy_t1,kxz_t1,kyy_t1,kyz_t1,kzz_t1,
     &     kxx_t2,kxy_t2,kxz_t2,kyy_t2,kyz_t2,kzz_t2,
     &     kxx_t3,kxy_t3,kxz_t3,kyy_t3,kyz_t3,kzz_t3

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

      CCTK_REAL zero,one,pi
      CCTK_REAL dth


c     ==================================================================
c
c     INITIALISATION
c
c     ==================================================================

c     Numbers.

      zero = 0.0D0
      one  = 1.0D0

      pi = acos(-one)

c     Grid.

      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     cartoon has special stencil size in y-direction (BB 2/99).

      stencil_size_y = cctk_nghostzones(2)
#ifdef THORN_CARTOON_2D
      if (CONTAINS('cartoon_active','yes')) then
         stencil_size_y = (cctk_lsh(2)-1)/2
      endif
#endif

c     On first iteration blank out the sources so we do not get FPEs
c     when we use (but clobber) them at boundaries later.

      if (cctk_iteration == 1) then

         adms_alp = zero

         adms_gxx = zero; adms_gxy = zero; adms_gxz = zero
         adms_gyy = zero; adms_gyz = zero; adms_gzz = zero

         adms_kxx = zero; adms_kxy = zero; adms_kxz = zero
         adms_kyy = zero; adms_kyz = zero; adms_kzz = zero

      end if


c     ==================================================================
c
c     FIRST TIME STEP
c
c     ==================================================================

      if (cctk_iteration == 1) then  

c        write(*,*) "INITIAL DATA for Staggered Leapfrog"

c        Set up sources of evolution eqs using g,K and matter at n=0
c        -----------------------------------------------------------

         CALL ADM_Lapse(CCTK_PASS_FTOF)

         if (ADM_evolveLapse) then
            do k=2,nz-1
               do j=2,ny-1
                  do i=2,nx-1
#include "macro/ALPSOURCES_guts.h"
#include "macro/ALPSOURCES_undefine.h"
                  end do
               end do
            end do
         end if

         if (evolveg_boundary.eq.1) then

            if ((shift_state.ne.SHIFT_INACTIVE)) then
               call CCTK_WARN(0,"evolveg_boundary = yes  is incompatible with a non-zero shift")
            end if

            do k=1,nz
               do j=1,ny
                  do i=1,nx
#include "macro/GSOURCES_guts.h"
#include "macro/GSOURCES_undefine.h"
                  end do
               end do
            end do

         else

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

         end if

         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        Evolve whole step variables from n=0 to n=-1
c        --------------------------------------------
c        Here we are differencing using FTCS. That is, it should
c        be only first order accurate in time.

         if (ADM_evolveLapse) then
            ADM_alp_p = alp - dt*adms_alp
         else
            ADM_alp_p = alp
         end if

         ADM_gxx_p = gxx - dt*adms_gxx
         ADM_gxy_p = gxy - dt*adms_gxy
         ADM_gxz_p = gxz - dt*adms_gxz
         ADM_gyy_p = gyy - dt*adms_gyy
         ADM_gyz_p = gyz - dt*adms_gyz
         ADM_gzz_p = gzz - dt*adms_gzz

c        Evolve extrinsic curvature from n=0 to n=1/2 and n=-1/2
c        -------------------------------------------------------
c        Here we are differencing using FTCS. That is, it should
c        be only first order accurate in time.

         dth = 0.5D0*dt

         ADM_kxx_p = kxx - dth*adms_kxx
         ADM_kxy_p = kxy - dth*adms_kxy
         ADM_kxz_p = kxz - dth*adms_kxz
         ADM_kyy_p = kyy - dth*adms_kyy
         ADM_kyz_p = kyz - dth*adms_kyz
         ADM_kzz_p = kzz - dth*adms_kzz

         ADM_kxx_stag = kxx + dth*adms_kxx
         ADM_kxy_stag = kxy + dth*adms_kxy
         ADM_kxz_stag = kxz + dth*adms_kxz
         ADM_kyy_stag = kyy + dth*adms_kyy
         ADM_kyz_stag = kyz + dth*adms_kyz
         ADM_kzz_stag = kzz + dth*adms_kzz

c        Synchronize evolved functions
c        -----------------------------

c        Note that the extrinsic curvature does not need
c        to be sychronized.

         call CCTK_SyncGroup(ierror,cctkGH,"adm::ADM_metric_prev")

c        Boundaries for evolved functions
c        --------------------------------

         if (ADM_evolveLapse) then
            call CCTK_SyncGroup(ierror,cctkGH,"adm::ADM_lapse_prev")
            call ADM_Boundary(cctkGH,"adm::ADM_lapse_prev")
         end if

         if (evolveg_boundary.eq.0) then
            call ADM_Boundary(cctkGH,"adm::ADM_metric_prev")
         end if

         call ADM_Boundary(cctkGH,"adm::ADM_curv_prev")
         call ADM_Boundary(cctkGH,"adm::ADM_curv_stag")

      end if

      if (cctk_iteration == 2) then
         if(CCTK_EQUALS(stagleap_firststep,"predcorr")) then

            write(*,*) "INITIAL DATA for SL use predcorr for step 2 setup"

c           copy h** original time step to ADM_h**_stag

            ADM_kxx_stag = ADM_kxx_p
            ADM_kxy_stag = ADM_kxy_p
            ADM_kxz_stag = ADM_kxz_p
            ADM_kyy_stag = ADM_kyy_p
            ADM_kyz_stag = ADM_kyz_p
            ADM_kzz_stag = ADM_kzz_p

c           copy g** original time step to adm_g**_extrap
c           Boy, it is a good thing stagleap is such a memory hog! :-)

            adm_gxx_extrap = ADM_gxx_p
            adm_gxy_extrap = ADM_gxy_p
            adm_gxz_extrap = ADM_gxz_p
            adm_gyy_extrap = ADM_gyy_p
            adm_gyz_extrap = ADM_gyz_p
            adm_gzz_extrap = ADM_gzz_p

c           now do another predictor corrector step

            call ADM_predcorr(CCTK_PASS_FTOF)

c           now we have h** on 3-tiered structure: average to
c           get staggered state for SL GHIteration=2

            do k=1,nz
               do j=1,ny
                  do i=1,nx

                     kxx_t1 = ADM_kxx_stag(i,j,k)
                     kxy_t1 = ADM_kxy_stag(i,j,k)
                     kxz_t1 = ADM_kxz_stag(i,j,k)
                     kyy_t1 = ADM_kyy_stag(i,j,k)
                     kyz_t1 = ADM_kyz_stag(i,j,k)
                     kzz_t1 = ADM_kzz_stag(i,j,k)

                     kxx_t2 = ADM_kxx_p(i,j,k)
                     kxy_t2 = ADM_kxy_p(i,j,k)
                     kxz_t2 = ADM_kxz_p(i,j,k)
                     kyy_t2 = ADM_kyy_p(i,j,k)
                     kyz_t2 = ADM_kyz_p(i,j,k)
                     kzz_t2 = ADM_kzz_p(i,j,k)

                     kxx_t3 = kxx(i,j,k)
                     kxy_t3 = kxy(i,j,k)
                     kxz_t3 = kxz(i,j,k)
                     kyy_t3 = kyy(i,j,k)
                     kyz_t3 = kyz(i,j,k)
                     kzz_t3 = kzz(i,j,k)

                     ADM_kxx_stag(i,j,k) = 0.5d0*(kxx_t3 + kxx_t2)
                     ADM_kxy_stag(i,j,k) = 0.5d0*(kxy_t3 + kxy_t2)
                     ADM_kxz_stag(i,j,k) = 0.5d0*(kxz_t3 + kxz_t2)
                     ADM_kyy_stag(i,j,k) = 0.5d0*(kyy_t3 + kyy_t2)
                     ADM_kyz_stag(i,j,k) = 0.5d0*(kyz_t3 + kyz_t2)
                     ADM_kzz_stag(i,j,k) = 0.5d0*(kzz_t3 + kzz_t2)

                     ADM_kxx_p(i,j,k) = 0.5d0*(kxx_t2 + kxx_t1)
                     ADM_kxy_p(i,j,k) = 0.5d0*(kxy_t2 + kxy_t1)
                     ADM_kxz_p(i,j,k) = 0.5d0*(kxz_t2 + kxz_t1)
                     ADM_kyy_p(i,j,k) = 0.5d0*(kyy_t2 + kyy_t1)
                     ADM_kyz_p(i,j,k) = 0.5d0*(kyz_t2 + kyz_t1)
                     ADM_kzz_p(i,j,k) = 0.5d0*(kzz_t2 + kzz_t1)

                     kxx(i,j,k) = kxx_t2
                     kxy(i,j,k) = kxy_t2
                     kxz(i,j,k) = kxz_t2
                     kyy(i,j,k) = kyy_t2
                     kyz(i,j,k) = kyz_t2
                     kzz(i,j,k) = kzz_t2

                  end do
               end do
            end do

c           now, put metric back to n=0 and n=1 timestep

            gxx = ADM_gxx_p
            gxy = ADM_gxy_p
            gxz = ADM_gxz_p
            gyy = ADM_gyy_p
            gyz = ADM_gyz_p
            gzz = ADM_gzz_p

            ADM_gxx_p = adm_gxx_extrap
            ADM_gxy_p = adm_gxy_extrap
            ADM_gxz_p = adm_gxz_extrap
            ADM_gyy_p = adm_gyy_extrap
            ADM_gyz_p = adm_gyz_extrap
            ADM_gzz_p = adm_gzz_extrap

         end if
      end if


c     ==================================================================
c
c     MAIN ITERATION LOOP 
c
c     ==================================================================
c
c     1. Evolve g and alpha to n+1
c
c        a) Extrapolate g and alpha to n+1/2 (using n, n-1)
c           (where we already have K = ADM_h_stag)
c        b) Relabel the new old timeslice
c        c) Rename metric and curvature for macros to work
c        d) Calculate sources for g and lapse evolution equations using 
c           K at n+1/2
c        e) Leapfrog to evolve g and lapse (boundaries and sychronise)
c
c     2. Evolve K to n+3/2
c
c        a) Extrapolate extrinsic curvature to n+1 (using n-1/2 and n+1/2)
c        b) Solve maximal slicing equation if required
c        c) Calculate sources for K evolution equations
c        d) Relabel the new old timeslice
c        e) Leapfrog to evolve K
c        f) Extrapolate K to n+1 (for analysis and output)
c
c     ==================================================================


c     EVOLVE 3-METRIC AND LAPSE
c     -------------------------

c     1a: Extrapolate g and alpha to n+1/2

      do k=1,nz
         do j=1,ny
            do i=1,nx
               adm_gxx_extrap(i,j,k) = 1.5D0*gxx(i,j,k) - 
     &                                 0.5D0*ADM_gxx_p(i,j,k)
               ADM_gxx_p(i,j,k)      = gxx(i,j,k)
               gxx(i,j,k)            = adm_gxx_extrap(i,j,k)
            end do
         end do
      end do

      do k=1,nz
         do j=1,ny
            do i=1,nx
               adm_gxy_extrap(i,j,k) = 1.5D0*gxy(i,j,k) 
     &                               - 0.5D0*ADM_gxy_p(i,j,k)
               ADM_gxy_p(i,j,k)      = gxy(i,j,k)
               gxy(i,j,k)            = adm_gxy_extrap(i,j,k)
            end do
         end do
      end do

      do k=1,nz
         do j=1,ny
            do i=1,nx
               adm_gxz_extrap(i,j,k) = 1.5D0*gxz(i,j,k) 
     &                               - 0.5D0*ADM_gxz_p(i,j,k)
               ADM_gxz_p(i,j,k)      = gxz(i,j,k)
               gxz(i,j,k)            = adm_gxz_extrap(i,j,k)
            end do
         end do
      end do

      do k=1,nz
         do j=1,ny
            do i=1,nx
               adm_gyy_extrap(i,j,k) = 1.5D0*gyy(i,j,k) 
     &                               - 0.5D0*ADM_gyy_p(i,j,k)
               ADM_gyy_p(i,j,k)      = gyy(i,j,k)
               gyy(i,j,k)            = adm_gyy_extrap(i,j,k)
            end do
         end do
      end do

      do k=1,nz
         do j=1,ny
            do i=1,nx
               adm_gyz_extrap(i,j,k) = 1.5D0*gyz(i,j,k) 
     &                               - 0.5D0*ADM_gyz_p(i,j,k)
               ADM_gyz_p(i,j,k)      = gyz(i,j,k)
               gyz(i,j,k)             = adm_gyz_extrap(i,j,k)
            end do
         end do
      end do

      do k=1,nz
         do j=1,ny
            do i=1,nx
               adm_gzz_extrap(i,j,k) = 1.5D0*gzz(i,j,k) 
     &                               - 0.5D0*ADM_gzz_p(i,j,k)
               ADM_gzz_p(i,j,k)      = gzz(i,j,k)
               gzz(i,j,k)            = adm_gzz_extrap(i,j,k)
            end do
         end do
      end do

      do k=1,nz
         do j=1,ny
            do i=1,nx
               adm_alp_extrap(i,j,k) = 1.5D0*alp(i,j,k) 
     &                               - 0.5D0*ADM_alp_p(i,j,k)
               ADM_alp_p(i,j,k)      = alp(i,j,k)
               alp(i,j,k)            = adm_alp_extrap(i,j,k)
            end do
         end do
      end do

      kxx = ADM_kxx_stag
      kxy = ADM_kxy_stag
      kxz = ADM_kxz_stag
      kyy = ADM_kyy_stag
      kyz = ADM_kyz_stag
      kzz = ADM_kzz_stag

c     1d: Calculate sources for metric and lapse evolution

      CALL ADM_Lapse(CCTK_PASS_FTOF)

      if (ADM_evolveLapse) then
         do k=2,nz-1
            do j=2,ny-1
               do i=2,nx-1
#include "macro/ALPSOURCES_guts.h"
#include "macro/ALPSOURCES_undefine.h"
               end do
            end do
         end do
      end if

      if (evolveg_boundary.eq.1) then

         if ((shift_state.ne.SHIFT_INACTIVE)) then
            call CCTK_WARN(0,"evolveg_boundary = yes  is incompatible with a non-zero shift")
         end if

         do k=1,nz
            do j=1,ny
               do i=1,nx
#include "macro/GSOURCES_guts.h"
#include "macro/GSOURCES_undefine.h"
               end do
            end do
         end do

      else

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

      end if

c     1e: Evolve metric and lapse.

      gxx = adms_gxx*dt + ADM_gxx_p 
      gxy = adms_gxy*dt + ADM_gxy_p
      gxz = adms_gxz*dt + ADM_gxz_p
      gyy = adms_gyy*dt + ADM_gyy_p
      gyz = adms_gyz*dt + ADM_gyz_p
      gzz = adms_gzz*dt + ADM_gzz_p

      if (ADM_evolveLapse) then
        alp = adms_alp*dt + ADM_alp_p
        call CCTK_SyncGroup(ierror,cctkGH,"einstein::lapse")
        call ADM_Boundary(cctkGH,"einstein::lapse")
      else
        alp = ADM_alp_p
      end if

      call CCTK_SyncGroup(ierror,cctkGH,"einstein::metric")

      if (evolveg_boundary.eq.0) then
         call ADM_Boundary(cctkGH,"einstein::metric")
      end if


c     EVOLVE EXTRINSIC CURVATURE
c     --------------------------

c     2a: Need the extrinsic curvature at n+1, so extrapolate from n-1/2 
c         and n+1/2

      kxx = 1.5D0*ADM_kxx_stag - 0.5D0*ADM_kxx_p
      kxy = 1.5D0*ADM_kxy_stag - 0.5D0*ADM_kxy_p
      kxz = 1.5D0*ADM_kxz_stag - 0.5D0*ADM_kxz_p
      kyy = 1.5D0*ADM_kyy_stag - 0.5D0*ADM_kyy_p
      kyz = 1.5D0*ADM_kyz_stag - 0.5D0*ADM_kyz_p
      kzz = 1.5D0*ADM_kzz_stag - 0.5D0*ADM_kzz_p

      end subroutine ADM_StaggeredLeapfrog1