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

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

      subroutine maximal(CCTK_ARGUMENTS)

      implicit none

      DECLARE_CCTK_ARGUMENTS
      DECLARE_CCTK_PARAMETERS

      integer i,j,k
      integer Mlin_index,Nsrc_index,field_index
      integer metpsi_index(7)
      integer CCTK_Equals
      integer Einstein_GetSlicingHandle
      integer ierr

      CCTK_REAL AbsTol(3),RelTol(3)

c     Do we really want to be here?

      if (mod(cctk_iteration,int(max_every)).ne.0) then
         return
      end if

c     Check if it is our turn to do maximal.

      if (active_slicing_handle.ne.Einstein_GetSlicingHandle("maximal")) then
         return
      end if

c     Perform some basic setup.

      call setupmaximal(CCTK_ARGUMENTS)

c     Prepare the call the linear elliptic, by deriving the
c     grid function indices.

      call CCTK_VarIndex (metpsi_index(1), "einstein::gxx")
      call CCTK_VarIndex (metpsi_index(2), "einstein::gxy")
      call CCTK_VarIndex (metpsi_index(3), "einstein::gxz")
      call CCTK_VarIndex (metpsi_index(4), "einstein::gyy")
      call CCTK_VarIndex (metpsi_index(5), "einstein::gyz")
      call CCTK_VarIndex (metpsi_index(6), "einstein::gzz")
      call CCTK_VarIndex (metpsi_index(7), "einstein::psi")
      call CCTK_VarIndex (field_index,     "einstein::alp")
      call CCTK_VarIndex (Mlin_index,      "maximal::Mlinear")
      call CCTK_VarIndex (Nsrc_index,      "maximal::Nsource")

      AbsTol(1)= maximal_thresh
      AbsTol(2)= -1
      AbsTol(3)= -1
      
      RelTol(1)= -1
      RelTol(2)= -1
      RelTol(3)= -1

c     Call the Linear Elliptic solver interface; depending on the 
c     chosen boudary for maximal, set the approriate values to the 
c     keys. We could do A LOT of error checking, which I do not do 
c     here: e.g. is a boundary provided?, is the solver provided?,
c     are the values supported?

c     Boundaries.

      if (CCTK_EQUALS(max_bound,"const")) then
         call Ell_SetRealKey(ierr,max_const_v0,
     $           "EllLinConfMetric::Bnd::Const::V0")
      end if

      if (CCTK_EQUALS(max_bound,"robin")) then
         call Ell_SetStrKey(ierr, "yes", "EllLinConfMetric::Bnd::Robin")
         if (ierr < 0) then
            call CCTK_WARN(0,"maximal: Setting EllLinConfMetric::Bnd::Robin failed")
         end if
         call Ell_SetIntKey(ierr,max_robin_falloff,
     $                      "EllLinConfMetric::Bnd::Robin::falloff")
         if (ierr < 0) then
            call CCTK_WARN(0,"maximal: Setting EllLinConfMetric::Bnd::Robin::falloff failed")
         end if
         call Ell_SetRealKey(ierr,max_robin_inf,
     $                       "EllLinConfMetric::Bnd::Robin::inf")
         if (ierr < 0) then
            call CCTK_WARN(0,"maximal: Setting EllLinConfMetric::Bnd::Robin::inf failed")
         end if
      else
         call Ell_SetStrKey(ierr, "no", "EllLinConfMetric::Bnd::Robin")
      end if

c     Solver.

      call CCTK_SyncGroup(ierr,cctkGH,"einstein::lapse")
      call CCTK_SyncGroup(ierr,cctkGH,"einstein::metric")

      if (CCTK_EQUALS(max_solver,"sor")) then

c        Set any SOR specific keys
         call Ell_SetIntKey(ierr,sor_maxit,"Ell::SORmaxit");
         call Ell_SetStrKey(ierr,"const","Ell::SORaccel");

c        Call SOR solver
         call Ell_LinConfMetricSolver(ierr, cctkGH, 
     $              metpsi_index,field_index,Mlin_index,Nsrc_index,
     $              AbsTol,RelTol,"sor")
      end if

      if (CCTK_EQUALS(max_solver,"bam")) then
         call Ell_LinConfMetricSolver(ierr, cctkGH, 
     $   metpsi_index,field_index,Mlin_index,Nsrc_index,
     $   AbsTol,RelTol,"bam")
      end if

      if (CCTK_EQUALS(max_solver,"petsc")) then
         call Ell_LinConfMetricSolver(ierr, cctkGH, 
     $   metpsi_index, field_index,Mlin_index,Nsrc_index,
     $   AbsTol,RelTol,"petsc")
      end if

c     Synchronization and symmetry boundaries.

      call CCTK_SyncGroup(ierr,cctkGH,"einstein::lapse")
      call CartSymGN(ierr,cctkGH,"einstein::lapse")
      
      return
      end