c Simplistic evolution using geodesic slicing and double leapfrog

/*@@
  @file      Evol.F
  @desc
     Double leapfrog stepper for the ADM Evolution.
  @enddesc
@@*/

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


/*@@
  @routine    ADM_DoubleLeap
  @desc
     Double leapfrog stepper for the ADM Evolution.
  @enddesc
  @calls
  @calledby
  @history
  @endhistory
@@*/

      subroutine EvolSimple_Evol(CCTK_ARGUMENTS)

      implicit none

      DECLARE_CCTK_ARGUMENTS
      DECLARE_CCTK_PARAMETERS
      DECLARE_CCTK_FUNCTIONS

      integer i,j,k
      integer nx,ny,nz
      integer ierr
      integer sw(3)
      
      CCTK_REAL dx,dy,dz,dt

c     ========================================================================
c
c     INITIAL STUFF
c
c     ========================================================================

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)

      sw(1) = 1
      sw(2) = 1
      sw(3) = 1

c     ========================================================================
c
c     CONSTRUCT SOURCE ARRAYS FOR ALL EVOLUTION EQUATIONS
c
c     ========================================================================

c     Geodesic slicing: Set lapse to one
      alp = 1D0

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

               adms_gxx(i,j,k) = -2.0*alp(i,j,k)*kxx_p(i,j,k)
               adms_gxy(i,j,k) = -2.0*alp(i,j,k)*kxy_p(i,j,k)
               adms_gxz(i,j,k) = -2.0*alp(i,j,k)*kxz_p(i,j,k)
               adms_gyy(i,j,k) = -2.0*alp(i,j,k)*kyy_p(i,j,k)
               adms_gyz(i,j,k) = -2.0*alp(i,j,k)*kyz_p(i,j,k)
               adms_gzz(i,j,k) = -2.0*alp(i,j,k)*kzz_p(i,j,k)

            end do
         end do
      end do

      call EvolSimple_KSources(CCTK_ARGUMENTS)
 
c     ========================================================================
c
c     FIRST EVOLUTION STEP
c
c     ========================================================================

      if (cctk_iteration == 1) then

c        Special treatment for the first timestep

c        Evolve first timestep using FTCS

            call CCTK_INFO("Initializing Leapfrog with FTCS Step")

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

                     kxx(i,j,k) = 0.0d0
                     kxy(i,j,k) = 0.0d0
                     kxz(i,j,k) = 0.0d0
                     kyy(i,j,k) = 0.0d0
                     kyz(i,j,k) = 0.0d0
                     kzz(i,j,k) = 0.0d0

                     gxx(i,j,k) = gxx_p(i,j,k)+dt*adms_gxx(i,j,k)
                     gxy(i,j,k) = gxy_p(i,j,k)+dt*adms_gxy(i,j,k)
                     gxz(i,j,k) = gxz_p(i,j,k)+dt*adms_gxz(i,j,k)
                     gyy(i,j,k) = gyy_p(i,j,k)+dt*adms_gyy(i,j,k)
                     gyz(i,j,k) = gyz_p(i,j,k)+dt*adms_gyz(i,j,k)
                     gzz(i,j,k) = gzz_p(i,j,k)+dt*adms_gzz(i,j,k)

                  end do
               end do
            end do

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

                     kxx(i,j,k) = kxx_p(i,j,k)+dt*adms_kxx(i,j,k)
                     kxy(i,j,k) = kxy_p(i,j,k)+dt*adms_kxy(i,j,k)
                     kxz(i,j,k) = kxz_p(i,j,k)+dt*adms_kxz(i,j,k)
                     kyy(i,j,k) = kyy_p(i,j,k)+dt*adms_kyy(i,j,k)
                     kyz(i,j,k) = kyz_p(i,j,k)+dt*adms_kyz(i,j,k)
                     kzz(i,j,k) = kzz_p(i,j,k)+dt*adms_kzz(i,j,k)

                  end do
               end do
            end do

      end if

c     ========================================================================
c
c     STANDARD EVOLUTION STEP
c
c     ========================================================================

      if (cctk_iteration .gt. 1) then

         call CCTK_INFO("Standard Evolution Step")

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

                  kxx(i,j,k) = 0.0d0
                  kxy(i,j,k) = 0.0d0
                  kxz(i,j,k) = 0.0d0
                  kyy(i,j,k) = 0.0d0
                  kyz(i,j,k) = 0.0d0
                  kzz(i,j,k) = 0.0d0
    
                  gxx(i,j,k) = gxx_p_p(i,j,k) + 2D0*dt*adms_gxx(i,j,k)
                  gxy(i,j,k) = gxy_p_p(i,j,k) + 2D0*dt*adms_gxy(i,j,k)
                  gxz(i,j,k) = gxz_p_p(i,j,k) + 2D0*dt*adms_gxz(i,j,k)
                  gyy(i,j,k) = gyy_p_p(i,j,k) + 2D0*dt*adms_gyy(i,j,k)
                  gyz(i,j,k) = gyz_p_p(i,j,k) + 2D0*dt*adms_gyz(i,j,k)
                  gzz(i,j,k) = gzz_p_p(i,j,k) + 2D0*dt*adms_gzz(i,j,k)

               end do
            end do
         end do

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

                  kxx(i,j,k) = kxx_p_p(i,j,k) + 2D0*dt*adms_kxx(i,j,k)
                  kxy(i,j,k) = kxy_p_p(i,j,k) + 2D0*dt*adms_kxy(i,j,k)
                  kxz(i,j,k) = kxz_p_p(i,j,k) + 2D0*dt*adms_kxz(i,j,k)
                  kyy(i,j,k) = kyy_p_p(i,j,k) + 2D0*dt*adms_kyy(i,j,k)
                  kyz(i,j,k) = kyz_p_p(i,j,k) + 2D0*dt*adms_kyz(i,j,k)
                  kzz(i,j,k) = kzz_p_p(i,j,k) + 2D0*dt*adms_kzz(i,j,k)

               end do
            end do
         end do

      end if


c     ========================================================================
c
c     END OF EVOLUTION STEP
c
c     ========================================================================

c     Apply flat boundaries
c      call BndFlatGN(ierr,cctkGH,sw,"admbase::curv")
      ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, -1,
     $     "admbase::curv", "Flat")
      if (ierr < 0) then
         call CCTK_WARN(1, "Error selecting admbase::curv for flat boundary condition")
      end if
     
      end