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

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

      subroutine setupmaximal(CCTK_ARGUMENTS)

      implicit none

      DECLARE_CCTK_ARGUMENTS
      DECLARE_CCTK_PARAMETERS

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

      integer, dimension(3) :: sym

      integer ierr
      integer ione
      integer, dimension(3) :: stencil

      CCTK_REAL :: dx,dy,dz,dt
      CCTK_REAL :: idx,idy,idz
      CCTK_REAL :: zero,half,one,two,four,pi,fourpi

      CCTK_REAL :: Ttt,Ttx,Tty,Ttz,Txx,Txy,Txz,Tyy,Tyz,Tzz,trace_T
      CCTK_REAL :: det,uxx,uxy,uxz,uyy,uyz,uzz

      CCTK_REAL :: finf

      CCTK_REAL, DIMENSION(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) :: K_temp

      logical maxshift 

c     Declarations for macros.

#include "CactusEinstein/Einstein/src/macro/TRK_declare.h"
#include "CactusEinstein/Einstein/src/macro/TRKK_declare.h"
#include "CactusEinstein/Einstein/src/macro/DETG_declare.h"
#include "CactusEinstein/Einstein/src/macro/UPPERMET_declare.h"
#include "CactusEinstein/Einstein/src/macro/TRT_declare.h"

c     Numbers.

      ione = 1

      zero = 0.0D0
      half = 0.5D0
      one  = 1.0D0
      two  = 2.0D0
      four = 4.0D0

      pi = 3.141592653589793D0

      fourpi = four*pi

c     Grid parameters.

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

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

      idx = one/dx
      idy = one/dy
      idz = one/dz

c     Define the symmetries for the maximal GFs.

      sym(1) = 1
      sym(2) = 1
      sym(3) = 1

      call SetCartSymVN(ierr,cctkGH,sym,'maximal::Mlinear')
      call SetCartSymVN(ierr,cctkGH,sym,'maximal::Nsource')

c     Define the stencil for the boundary condition.

      stencil(1) = 1
      stencil(2) = 1
      stencil(3) = 1

c     Calculate the coefficient matrices M and N.

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

c              Compute source term: Mlinear = -K_{ij} K^{ij}.

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

               Mlinear(i,j,k) = - TRKK_TRKK

c              This may be needed for CalcTmunu.

#include "CactusEinstein/Einstein/src/macro/DETG_guts.h"
#include "CactusEinstein/Einstein/src/macro/UPPERMET_guts.h"

               det = DETG_DETG

               uxx = UPPERMET_UXX; uxy = UPPERMET_UXY; uxz = UPPERMET_UXZ
               uyy = UPPERMET_UYY; uyz = UPPERMET_UYZ; uzz = UPPERMET_UZZ

c              Matter terms (still missing).

               Ttt = 0.0D0; Ttx = 0.0D0; Tty = 0.0D0; Ttz = 0.0D0

               Txx = 0.0D0; Txy = 0.0D0; Txz = 0.0D0
               Tyy = 0.0D0; Tyz = 0.0D0; Tzz = 0.0D0

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

c              Add matter terms.

               trace_T = TRT_TRT

               Mlinear(i,j,k) = Mlinear(i,j,k) - fourpi
     .            *(trace_T + two*Ttt/alp(i,j,k)**2)

c              Shift terms for matter.

               if (shift_state.ne.SHIFT_INACTIVE) then

                  Mlinear(i,j,k) = Mlinear(i,j,k) + fourpi*
     .               ((betax(i,j,k)*Ttx
     .               + betay(i,j,k)*Tty
     .               + betaz(i,j,k)*Ttz)*two
     .               - betax(i,j,k)*betax(i,j,k)*Txx
     .               - betay(i,j,k)*betay(i,j,k)*Tyy
     .               - betaz(i,j,k)*betaz(i,j,k)*Tzz
     .               -(betax(i,j,k)*betay(i,j,k)*Txy
     .               - betax(i,j,k)*betaz(i,j,k)*Txz
     .               - betay(i,j,k)*betaz(i,j,k)*Tyz)*two)
     .               /alp(i,j,k)**2

               end if

c              Zero the source matrix.

               Nsource(i,j,k) = zero

            end do
         end do
      end do
      
c     Undefine the Einstein macros.

#include "CactusEinstein/Einstein/src/macro/TRKK_undefine.h"
#include "CactusEinstein/Einstein/src/macro/DETG_guts.h"
#include "CactusEinstein/Einstein/src/macro/UPPERMET_undefine.h"
#include "CactusEinstein/Einstein/src/macro/TRT_undefine.h"

c     Do we want shift terms?

      maxshift = (max_shift.eq.1)

c     Add the shift term: N = B^i D_i(trK).

      if ((maxshift).and.(shift_state.eq.SHIFT_ACTIVE)) then

         do k=1,nz
            do j=1,ny
               do i=1,nx
#include "CactusEinstein/Einstein/src/macro/TRK_guts.h"
                  K_temp(i,j,k) = TRK_TRK
	       end do
             end do
         end do

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

c     Note the loop boundaries, as we take derivs.

         do k=2,nz-1
            do j=2,ny-1
               do i=2,nx-1
                  Nsource(i,j,k) = Nsource(i,j,k) - half
     .              *(betax(i,j,k)*(K_temp(i+1,j,k)-K_temp(i-1,j,k))*idx
     .              - betay(i,j,k)*(K_temp(i,j+1,k)-K_temp(i,j-1,k))*idy
     .              - betaz(i,j,k)*(K_temp(i,j,k+1)-K_temp(i,j,k-1))*idz)
               end do
            end do
         end do

      end if

c     Synchronization and boundaries.

      call CCTK_SyncGroup(ierr,cctkGH,'maximal::maximal_coeff')
c     call ADM_Boundary(cctkGH,'maximal::maximal_coeff')

      call CartSymGN(ierr, cctkGH, 'maximal::maximal_coeff')
      finf = 0
      call BndRobinGN(ierr, cctkGH, stencil, finf, ione, 
     .     'maximal::maximal_coeff')
      call CCTK_SyncGroup(ierr,cctkGH,'maximal::maximal_coeff')

      return
      end