/*@@
  @file      DoubleLeap.F
  @desc
     Double 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_DoubleLeap
  @desc
     Double leapfrog stepper for the ADM Evolution.
  @enddesc
  @calls
  @calledby
  @history
  @endhistory
@@*/

      subroutine ADM_DoubleLeap(CCTK_ARGUMENTS)

      USE ADM_Scalars

      implicit none

      DECLARE_CCTK_ARGUMENTS
      DECLARE_CCTK_PARAMETERS
      DECLARE_CCTK_FUNCTIONS

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

      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 "CactusEinstein/ADM/src/macro/ALPSOURCES_declare.h"
#include "CactusEinstein/ADM/src/macro/KSOURCES_declare.h"
#include "CactusEinstein/ADM/src/macro/GSOURCES_declare.h"
#include "CactusEinstein/ADM/src/macro/DOUBLELEAP_ALP_declare.h"
#include "CactusEinstein/ADM/src/macro/DOUBLELEAP_G_declare.h"
#include "CactusEinstein/ADM/src/macro/DOUBLELEAP_K_declare.h"
#include "CactusEinstein/ADM/src/macro/FTCSSTEP_ALP_declare.h"
#include "CactusEinstein/ADM/src/macro/FTCSSTEP_G_declare.h"
#include "CactusEinstein/ADM/src/macro/FTCSSTEP_K_declare.h"

      CCTK_REAL :: pi,tem
      CCTK_REAL :: dx,dy,dz,dt


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

      pi = ACOS(-1.0D0)

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 (cartoon_active) then
         stencil_size_y = (cctk_lsh(2)-1)/2
      endif
#endif

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

      if (cctk_iteration == 1) then

         adms_alp = 0D0

         adms_gxx = 0D0; adms_gxy = 0D0; adms_gxz = 0D0
         adms_gyy = 0D0; adms_gyz = 0D0; adms_gzz = 0D0

         adms_kxx = 0D0; adms_kxy = 0D0; adms_kxz = 0D0
         adms_kyy = 0D0; adms_kyz = 0D0; adms_kzz = 0D0

      end if


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

c     [3-metric, extrinsic curvature, lapse]
c     Store them in the arrays: adms_kzz, adms_gxx, adms_alp

      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 "CactusEinstein/ADM/src/macro/ALPSOURCES_guts.h"
               end do
            end do
         end do
      end if

#include "CactusEinstein/ADM/src/macro/ALPSOURCES_undefine.h"

      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 "CactusEinstein/ADM/src/macro/GSOURCES_guts.h"
               end do
            end do
         end do

#include "CactusEinstein/ADM/src/macro/GSOURCES_undefine.h"

      else

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

#include "CactusEinstein/ADM/src/macro/GSOURCES_undefine.h"

      end if

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

#include "CactusEinstein/ADM/src/macro/KSOURCES_undefine.h"


c     ========================================================================
c
c     FIRST EVOLUTION STEP
c
c     ========================================================================

      if (cctk_iteration == 1) then

c        Special treatment for the first timestep

         if((CCTK_EQUALS(leapfrog_firststep,"ftcs")) .or. 
     &       CCTK_EQUALS(leapfrog_firststep,"predcorr")) then

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
#include "CactusEinstein/ADM/src/macro/FTCSSTEP_ALP_guts.h"
#include "CactusEinstein/ADM/src/macro/FTCSSTEP_G_guts.h"
#include "CactusEinstein/ADM/src/macro/FTCSSTEP_K_guts.h"
                  end do
               end do
            end do

#include "CactusEinstein/ADM/src/macro/FTCSSTEP_ALP_undefine.h"
#include "CactusEinstein/ADM/src/macro/FTCSSTEP_G_undefine.h"
#include "CactusEinstein/ADM/src/macro/FTCSSTEP_K_undefine.h"

         end if

         if(CCTK_EQUALS(leapfrog_firststep,"predcorr")) then

c           Do a corrector step

            write (*,*) "Using a corrector step"

c           Sync evolved grid functions and apply boundary conditions
c           Note that we do not need to sync the extrinsic curvature

            if (ADM_evolveLapse) then
               call CCTK_SyncGroup(ierror,cctkGH,"einstein::lapse")
               call ADM_Boundary(cctkGH,"einstein::lapse")   
            end if

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

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

c           Recalculate the sources
c           -----------------------

c           Get coefficients for 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 "CactusEinstein/ADM/src/macro/ALPSOURCES_guts.h"
                     end do
                  end do
               end do
            end if

#include "CactusEinstein/ADM/src/macro/ALPSOURCES_undefine.h"

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

#include "CactusEinstein/ADM/src/macro/GSOURCES_undefine.h"

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

#include "CactusEinstein/ADM/src/macro/KSOURCES_undefine.h"

c           Now do the corrector step.

            if (ADM_evolveLapse) then
               alp = 0.5D0*(alp + ADM_alp_p + dt*adms_alp)
            end if

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

            kxx = 0.5D0*(kxx + ADM_kxx_p + dt*adms_kxx)
            kyy = 0.5D0*(kyy + ADM_kyy_p + dt*adms_kyy)
            kzz = 0.5D0*(kzz + ADM_kzz_p + dt*adms_kzz)
            kxy = 0.5D0*(kxy + ADM_kxy_p + dt*adms_kxy)
            kxz = 0.5D0*(kxz + ADM_kxz_p + dt*adms_kxz)
            kyz = 0.5D0*(kyz + ADM_kyz_p + dt*adms_kyz)

	end if

      end if


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

      if (cctk_iteration > 1) then

         if (ADM_evolveLapse) then
            do k=1,nz
               do j=1,ny
                  do i=1,nx
                     if ((ADM_flatNablaAlpCoeff.eq.0).and.
     .                   (ADM_nablaAlpCoeff.eq.0)) then
#include "CactusEinstein/ADM/src/macro/DOUBLELEAP_ALP_guts.h"
                     else
#include "CactusEinstein/ADM/src/macro/FTCSSTEP_ALP_guts.h"
                     endif
                  end do
               end do
            end do
         end if

#include "CactusEinstein/ADM/src/macro/DOUBLELEAP_ALP_undefine.h"
#include "CactusEinstein/ADM/src/macro/FTCSSTEP_ALP_undefine.h"

         do k=1,nz
            do j=1,ny
               do i=1,nx
#include "CactusEinstein/ADM/src/macro/DOUBLELEAP_G_guts.h"
               end do
            end do
         end do

#include "CactusEinstein/ADM/src/macro/DOUBLELEAP_G_undefine.h"

         do k=2,nz-1
            do j=2,ny-1
               do i=2,nx-1
#include "CactusEinstein/ADM/src/macro/DOUBLELEAP_K_guts.h"
               end do
            end do
         end do

#include "CactusEinstein/ADM/src/macro/DOUBLELEAP_K_undefine.h"

      end if


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

c     Synchronize the evolved grid functions.

      if (ADM_evolveLapse) then
         call CCTK_SyncGroup(ierror,cctkGH,"einstein::lapse")
      end if

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

c     Boundaries.

      if (ADM_evolveLapse) then
         call ADM_Boundary(cctkGH,"einstein::lapse")
      end if

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

      call ADM_Boundary(cctkGH,"einstein::curv")

      end subroutine ADM_DoubleLeap