/*@@
   @file      GRHydro_BoundariesM.F90
   @date      Aug 30, 2010
   @author    
   @desc 
   The two routines for dealing with boundary conditions.
   @enddesc 
 @@*/

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

#include "util_Table.h"

#define velx(i,j,k) vel(i,j,k,1)
#define vely(i,j,k) vel(i,j,k,2)
#define velz(i,j,k) vel(i,j,k,3)
#define sx(i,j,k) scon(i,j,k,1)
#define sy(i,j,k) scon(i,j,k,2)
#define sz(i,j,k) scon(i,j,k,3)
#define Bvecx(i,j,k) Bvec(i,j,k,1)
#define Bvecy(i,j,k) Bvec(i,j,k,2)
#define Bvecz(i,j,k) Bvec(i,j,k,3)

 /*@@
   @routine    GRHydro_InitSymBoundM
   @date       Aug 30, 2010
   @author     Joshua Faber, Ian Hawke
   @desc 
   Sets up the symmetries at the boundaries of the hydrodynamical variables.
   @enddesc 
   @calls     
   @calledby   
   @history 
   Direct translation of routines from GR3D, GRAstro_Hydro, 
   written by Mark Miller, or WaveToy routines, or...
   @endhistory 

@@*/

subroutine GRHydro_InitSymBoundM(CCTK_ARGUMENTS)

  implicit none
  
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS
  DECLARE_CCTK_FUNCTIONS

  integer, dimension(3) :: sym
  integer :: ierr
  integer :: itracer

  character(len=100) tracername
  character(len=100) tracerindex

  sym(1) = 1
  sym(2) = 1
  sym(3) = 1
  
  call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::rho")
  call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::press")
  call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::dens")
  call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::tau")
  call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::w_lorentz")
  call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::eps")
  call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::GRHydro_C2P_failed")
  if(clean_divergence.ne.0) then
     call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::psidc")
  endif

!!$ handle multiple tracer variables
  if(evolve_tracer.eq.1) then
     call SetCartSymGN(ierr, cctkGH, sym, "GRHydro::GRHydro_tracers")
     call SetCartSymGN(ierr, cctkGH, sym, "GRHydro::GRHydro_cons_tracers")
  endif 

  if(CCTK_EQUALS(Y_e_evolution_method,"GRHydro")) then
     call SetCartSymGN(ierr, cctkGH, sym, "HydroBase::Y_e")
     call SetCartSymGN(ierr, cctkGH, sym, "GRHydro::Y_e_con")
  endif
  
  sym(1) = -1
  sym(2) = 1
  sym(3) = 1
  
  call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::vel[0]")
  call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::scon[0]")
  call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::Bvec[0]")

  
  sym(1) = 1
  sym(2) = -1
  sym(3) = 1
  
  call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::vel[1]")
  call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::scon[1]")
  call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::Bvec[1]")
  
  sym(1) = 1
  sym(2) = 1
  sym(3) = -1
  
  call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::vel[2]")
  call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::scon[2]")
  call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::Bvec[2]")
  
end subroutine GRHydro_InitSymBoundM

 /*@@
   @routine    GRHydro_BoundariesM
   @date       Aug 30, 2010
   @author     
   @desc 
   Calls the appropriate boundary routines
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/

subroutine GRHydro_BoundariesM(CCTK_ARGUMENTS)
  
  implicit none
  
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS
  DECLARE_CCTK_FUNCTIONS
  
  integer, dimension(3) :: sw
  integer :: ierr
  CCTK_REAL :: pi = 3.141569d0
  integer :: i,j,k

  CCTK_INT, parameter :: faces=CCTK_ALL_FACES
  CCTK_INT, parameter :: ione=1
  
  sw = GRHydro_stencil

!!$Flat boundaries if required  

  if (CCTK_EQUALS(bound,"flat")) then
    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
         "GRHydro::dens", "Flat")
    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
         "GRHydro::tau", "Flat")
    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
         "GRHydro::scon", "Flat")
    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
         "HydroBase::w_lorentz", "Flat")
    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
         "HydroBase::rho", "Flat")
    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
         "HydroBase::press", "Flat")
    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
         "HydroBase::eps", "Flat")
    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
         "HydroBase::vel", "Flat")
    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
         "HydroBase::Bvec", "Flat")
    if(clean_divergence.ne.0) then
       ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
            "GRHydro::psidc", "Flat")
    endif

    if(evolve_tracer .ne. 0) then 
      ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
            "GRHydro::GRHydro_tracers", "Flat")
      ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
            "GRHydro::GRHydro_cons_tracers", "Flat")
    endif

    if(CCTK_EQUALS(Y_e_evolution_method,"GRHydro")) then
      ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
            "HydroBase::Y_e", "Flat")
      ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
            "GRHydro::Y_e_con", "Flat")
    endif

  endif

  if (CCTK_EQUALS(bound,"none")) then
    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
         "GRHydro::dens", "None")
    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
         "GRHydro::tau", "None")
    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
         "GRHydro::scon", "None")
    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
         "HydroBase::w_lorentz", "None")
    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
         "HydroBase::rho", "None")
    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
         "HydroBase::press", "None")
    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
         "HydroBase::eps", "None")
    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
         "HydroBase::vel", "None")
    ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
         "HydroBase::Bvec", "None")
    if(clean_divergence.ne.0) then
       ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
            "GRHydro::psidc", "None")
    endif

    if(evolve_tracer .ne. 0) then 
      ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
            "GRHydro::GRHydro_tracers", "None")
      ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
            "GRHydro::GRHydro_cons_tracers", "None")
    endif

    if(CCTK_EQUALS(Y_e_evolution_method,"GRHydro")) then
      ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
            "HydroBase::Y_e", "None")
      ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
            "GRHydro::Y_e_con", "None")
    endif

  end if

  if (CCTK_EQUALS(bound,"scalar")) then
    call CCTK_WARN(0, "Until somebody uses this I see no reason to support it")
  end if

  if (ierr < 0) call CCTK_WARN(0, "problems with applying the chosen boundary condition")

end subroutine GRHydro_BoundariesM