/*@@
   @file      GRHydro_HLLE.F90
   @date      Sat Jan 26 01:40:14 2002
   @author    Ian Hawke, Pedro Montero, Toni Font
   @desc 
   The HLLE solver. Called from the wrapper function, so works in 
   all directions.
   @enddesc 
 @@*/
   
#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
#include "cctk_Functions.h"

#include "GRHydro_Macros.h"
#include "SpaceMask.h"

 /*@@
   @routine    GRHydro_HLLE
   @date       Sat Jan 26 01:41:02 2002
   @author     Ian Hawke, Pedro Montero, Toni Font
   @desc 
   The HLLE solver. Sufficiently simple that its just one big routine.
   @enddesc 
   @calls     
   @calledby   
   @history 
   Altered from Cactus 3 routines originally written by Toni Font.
   @endhistory 

@@*/

subroutine GRHydro_HLLE(CCTK_ARGUMENTS)
  USE GRHydro_Eigenproblem

  implicit none
  
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS
  
  integer :: i, j, k, m
  CCTK_REAL, dimension(5) :: consp,consm_i,fplus,fminus,lamplus
  CCTK_REAL, dimension(5) :: f1,lamminus
  CCTK_REAL, dimension(5) :: qdiff
  CCTK_REAL ::  charmin, charmax, charpm,avg_alp,avg_det
  CCTK_REAL :: gxxh, gxyh, gxzh, gyyh, gyzh, gzzh, uxxh, uxyh, &
       uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r
  CCTK_REAL :: xtemp
    
  CCTK_INT :: type_bits, trivial, not_trivial

  if (flux_direction == 1) then
    call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemX")
    call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemX", &
         &"trivial")
  else if (flux_direction == 2) then
    call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemY")
    call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemY", &
         &"trivial")
  else if (flux_direction == 3) then
    call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemZ")
    call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemZ", &
         &"trivial")
  else
    call CCTK_WARN(0, "Flux direction not x,y,z")
  end if

  !$OMP PARALLEL DO PRIVATE(j,i,f1,lamminus,lamplus,consp,consm_i, &
  !$OMP                     fplus,fminus,qdiff,avg_beta,avg_alp,   &
  !$OMP                     avg_det,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, &
  !$OMP                     uxxh,uxyh,uxzh,uyyh,uyzh,uzzh,         &
  !$OMP                     usendh, xtemp, charmin, charmax, charpm, &
  !$OMP                     m)
  do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil
    do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil
      do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil
        
        f1 = 0.d0
        lamminus = 0.d0
        lamplus = 0.d0
        consp = 0.d0
        consm_i = 0.d0
        fplus = 0.d0
        fminus = 0.d0
        qdiff = 0.d0
        
!!$        Set the left (p for plus) and right (m_i for minus, i+1) states
        
        consp(1)   = densplus(i,j,k) 
        consp(2)   = sxplus(i,j,k)
        consp(3)   = syplus(i,j,k)
        consp(4)   = szplus(i,j,k)
        consp(5)   = tauplus(i,j,k)
        
        consm_i(1) = densminus(i+xoffset,j+yoffset,k+zoffset)
        consm_i(2) = sxminus(i+xoffset,j+yoffset,k+zoffset)
        consm_i(3) = syminus(i+xoffset,j+yoffset,k+zoffset)
        consm_i(4) = szminus(i+xoffset,j+yoffset,k+zoffset)
        consm_i(5) = tauminus(i+xoffset,j+yoffset,k+zoffset) 
          
!!$        Calculate various metric terms here.
!!$        Note also need the average of the lapse at the 
!!$        left and right points.

        if (shift_state .ne. 0) then
          if (flux_direction == 1) then
            avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
                 betax(i,j,k))
          else if (flux_direction == 2) then
            avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
                 betay(i,j,k))
          else if (flux_direction == 3) then
            avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
                 betaz(i,j,k))
          else
            call CCTK_WARN(0, "Flux direction not x,y,z")
          end if
        else
          avg_beta = 0.d0
        end if

        avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset))
        
        gxxh = 0.5d0 * (gxx(i+xoffset,j+yoffset,k+zoffset) + gxx(i,j,k))
        gxyh = 0.5d0 * (gxy(i+xoffset,j+yoffset,k+zoffset) + gxy(i,j,k))
        gxzh = 0.5d0 * (gxz(i+xoffset,j+yoffset,k+zoffset) + gxz(i,j,k))
        gyyh = 0.5d0 * (gyy(i+xoffset,j+yoffset,k+zoffset) + gyy(i,j,k))
        gyzh = 0.5d0 * (gyz(i+xoffset,j+yoffset,k+zoffset) + gyz(i,j,k))
        gzzh = 0.5d0 * (gzz(i+xoffset,j+yoffset,k+zoffset) + gzz(i,j,k))

        avg_det =  SPATIAL_DETERMINANT(gxxh,gxyh,gxzh,\
             gyyh,gyzh,gzzh)

!!$ If the Riemann problem is trivial, just calculate the fluxes from the 
!!$ left state and skip to the next cell
          
        if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, trivial)) then

          if (flux_direction == 1) then
            call num_x_flux(consp(1),consp(2),consp(3),consp(4),consp(5),&
                 f1(1),f1(2),f1(3),&
                 f1(4),f1(5),&
                 velxplus(i,j,k),pressplus(i,j,k),&
                 avg_det,avg_alp,avg_beta)
          else if (flux_direction == 2) then
            call num_x_flux(consp(1),consp(3),consp(4),consp(2),consp(5),&
                 f1(1),f1(3),f1(4),&
                 f1(2),f1(5),&
                 velyplus(i,j,k),pressplus(i,j,k),&
                 avg_det,avg_alp,avg_beta)
          else if (flux_direction == 3) then
            call num_x_flux(consp(1),consp(4),consp(2),consp(3),consp(5),&
                 f1(1),f1(4),f1(2),&
                 f1(3),f1(5),&
                 velzplus(i,j,k),pressplus(i,j,k),&
                 avg_det,avg_alp,avg_beta)
          else
            call CCTK_WARN(0, "Flux direction not x,y,z")
          end if
          
        else !!! The end of this branch is right at the bottom of the routine
            
          call UpperMetric(uxxh, uxyh, uxzh, uyyh, uyzh, uzzh, &
               avg_det,gxxh, gxyh, gxzh, &
               gyyh, gyzh, gzzh)
          
          if (flux_direction == 1) then
            usendh = uxxh
          else if (flux_direction == 2) then
            usendh = uyyh
          else if (flux_direction == 3) then
            usendh = uzzh
          else
            call CCTK_WARN(0, "Flux direction not x,y,z")
          end if
          
!!$        Calculate the jumps in the conserved variables
          
          qdiff(1) = consm_i(1) - consp(1)
          qdiff(2) = consm_i(2) - consp(2)
          qdiff(3) = consm_i(3) - consp(3)
          qdiff(4) = consm_i(4) - consp(4)
          qdiff(5) = consm_i(5) - consp(5)
          
!!$        Eigenvalues and fluxes either side of the cell interface
          
          if (flux_direction == 1) then
             if(evolve_temper.ne.1) then
                call eigenvalues(GRHydro_eos_handle,&
                     rhominus(i+xoffset,j+yoffset,k+zoffset),&
                     velxminus(i+xoffset,j+yoffset,k+zoffset),&
                     velyminus(i+xoffset,j+yoffset,k+zoffset),&
                     velzminus(i+xoffset,j+yoffset,k+zoffset),&
                     epsminus(i+xoffset,j+yoffset,k+zoffset),&
                     w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),&
                     lamminus,gxxh,gxyh,gxzh,gyyh,&
                     gyzh,gzzh,&
                     usendh,avg_alp,avg_beta)
                call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),&
                     velxplus(i,j,k),velyplus(i,j,k),&
                     velzplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),&
                     lamplus,gxxh,gxyh,gxzh,gyyh,&
                     gyzh,gzzh,&
                     usendh,avg_alp,avg_beta)
             else
                xtemp = temperature(i,j,k)
                call eigenvalues_hot(GRHydro_eos_handle,i,j,k,&
                     rhominus(i+xoffset,j+yoffset,k+zoffset),&
                     velxminus(i+xoffset,j+yoffset,k+zoffset),&
                     velyminus(i+xoffset,j+yoffset,k+zoffset),&
                     velzminus(i+xoffset,j+yoffset,k+zoffset),&
                     epsminus(i+xoffset,j+yoffset,k+zoffset),&
                     xtemp,&
                     y_e_minus(i+xoffset,j+yoffset,k+zoffset),&
                     w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),&
                     lamminus,gxxh,gxyh,gxzh,gyyh,&
                     gyzh,gzzh,&
                     usendh,avg_alp,avg_beta)
                xtemp = temperature(i,j,k)
                call eigenvalues_hot(GRHydro_eos_handle,i,j,k,&
                     rhoplus(i,j,k),&
                     velxplus(i,j,k),velyplus(i,j,k),&
                     velzplus(i,j,k),epsplus(i,j,k), &
                     xtemp,y_e_plus(i,j,k),&
                     w_lorentzplus(i,j,k),&
                     lamplus,gxxh,gxyh,gxzh,gyyh,&
                     gyzh,gzzh,&
                     usendh,avg_alp,avg_beta)
             endif
            call num_x_flux(consp(1),consp(2),consp(3),consp(4),consp(5),&
                 fplus(1),fplus(2),fplus(3),fplus(4),&
                 fplus(5),velxplus(i,j,k),pressplus(i,j,k),&
                 avg_det,avg_alp,avg_beta)
            call num_x_flux(consm_i(1),consm_i(2),consm_i(3),&
                 consm_i(4),consm_i(5),fminus(1),fminus(2),fminus(3),&
                 fminus(4), fminus(5),&
                 velxminus(i+xoffset,j+yoffset,k+zoffset),&
                 pressminus(i+xoffset,j+yoffset,k+zoffset),&
                 avg_det,avg_alp,avg_beta)
          else if (flux_direction == 2) then
             if(evolve_temper.ne.1) then
                call eigenvalues(GRHydro_eos_handle,&
                     rhominus(i+xoffset,j+yoffset,k+zoffset),&
                     velyminus(i+xoffset,j+yoffset,k+zoffset),&
                     velzminus(i+xoffset,j+yoffset,k+zoffset),&
                     velxminus(i+xoffset,j+yoffset,k+zoffset),&
                     epsminus(i+xoffset,j+yoffset,k+zoffset),&
                     w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),&
                     lamminus,gyyh,gyzh,gxyh,gzzh,&
                     gxzh,gxxh,&
                     usendh,avg_alp,avg_beta)
                call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),&
                     velyplus(i,j,k),velzplus(i,j,k),&
                     velxplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),&
                     lamplus,gyyh,gyzh,gxyh,gzzh,&
                     gxzh,gxxh,&
                     usendh,avg_alp,avg_beta)
             else
                xtemp = temperature(i,j,k)
                call eigenvalues_hot(GRHydro_eos_handle,i,j,k,&
                     rhominus(i+xoffset,j+yoffset,k+zoffset),&
                     velyminus(i+xoffset,j+yoffset,k+zoffset),&
                     velzminus(i+xoffset,j+yoffset,k+zoffset),&
                     velxminus(i+xoffset,j+yoffset,k+zoffset),&
                     epsminus(i+xoffset,j+yoffset,k+zoffset),&
                     xtemp,&
                     y_e_minus(i+xoffset,j+yoffset,k+zoffset),&
                     w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),&
                     lamminus,gyyh,gyzh,gxyh,gzzh,&
                     gxzh,gxxh,&
                     usendh,avg_alp,avg_beta)
                xtemp = temperature(i,j,k)
                call eigenvalues_hot(GRHydro_eos_handle,i,j,k,&
                     rhoplus(i,j,k),&
                     velyplus(i,j,k),velzplus(i,j,k),&
                     velxplus(i,j,k),epsplus(i,j,k),&
                     xtemp,y_e_plus(i,j,k),&
                     w_lorentzplus(i,j,k),&
                     lamplus,gyyh,gyzh,gxyh,gzzh,&
                     gxzh,gxxh,&
                     usendh,avg_alp,avg_beta)
             endif
            call num_x_flux(consp(1),consp(3),consp(4),consp(2),consp(5),&
                 fplus(1),fplus(3),fplus(4),fplus(2),&
                 fplus(5),velyplus(i,j,k),pressplus(i,j,k),&
                 avg_det,avg_alp,avg_beta)
            call num_x_flux(consm_i(1),consm_i(3),consm_i(4),&
                 consm_i(2),consm_i(5),fminus(1),fminus(3),fminus(4),&
                 fminus(2), fminus(5),&
                 velyminus(i+xoffset,j+yoffset,k+zoffset),&
                 pressminus(i+xoffset,j+yoffset,k+zoffset),&
                 avg_det,avg_alp,avg_beta)
          else if (flux_direction == 3) then
             if(evolve_temper.ne.1) then
                call eigenvalues(GRHydro_eos_handle,&
                     rhominus(i+xoffset,j+yoffset,k+zoffset),&
                     velzminus(i+xoffset,j+yoffset,k+zoffset),&
                     velxminus(i+xoffset,j+yoffset,k+zoffset),&
                     velyminus(i+xoffset,j+yoffset,k+zoffset),&
                     epsminus(i+xoffset,j+yoffset,k+zoffset),&
                     w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),&
                     lamminus,gzzh,gxzh,gyzh,&
                     gxxh,gxyh,gyyh,&
                     usendh,avg_alp,avg_beta)
                call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),&
                     velzplus(i,j,k),velxplus(i,j,k),&
                     velyplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),&
                     lamplus,gzzh,gxzh,gyzh,&
                     gxxh,gxyh,gyyh,&
                     usendh,avg_alp,avg_beta)
             else
                xtemp = temperature(i,j,k)
                call eigenvalues_hot(GRHydro_eos_handle,i,j,k,&
                     rhominus(i+xoffset,j+yoffset,k+zoffset),&
                     velzminus(i+xoffset,j+yoffset,k+zoffset),&
                     velxminus(i+xoffset,j+yoffset,k+zoffset),&
                     velyminus(i+xoffset,j+yoffset,k+zoffset),&
                     epsminus(i+xoffset,j+yoffset,k+zoffset),&
                     xtemp,y_e_minus(i+xoffset,j+yoffset,k+zoffset),&
                     w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),&
                     lamminus,gzzh,gxzh,gyzh,&
                     gxxh,gxyh,gyyh,&
                     usendh,avg_alp,avg_beta)
                xtemp = temperature(i,j,k)
                call eigenvalues_hot(GRHydro_eos_handle,i,j,k,&
                     rhoplus(i,j,k),&
                     velzplus(i,j,k),velxplus(i,j,k),&
                     velyplus(i,j,k),epsplus(i,j,k),&
                     xtemp,y_e_plus(i,j,k),&
                     w_lorentzplus(i,j,k),&
                     lamplus,gzzh,gxzh,gyzh,&
                     gxxh,gxyh,gyyh,&
                     usendh,avg_alp,avg_beta)
             endif
            call num_x_flux(consp(1),consp(4),consp(2),consp(3),consp(5),&
                 fplus(1),fplus(4),fplus(2),fplus(3),&
                 fplus(5),velzplus(i,j,k),pressplus(i,j,k),avg_det,&
                 avg_alp,avg_beta)
            call num_x_flux(consm_i(1),consm_i(4),consm_i(2),&
                 consm_i(3),consm_i(5),fminus(1),fminus(4),fminus(2),&
                 fminus(3), fminus(5),&
                 velzminus(i+xoffset,j+yoffset,k+zoffset),&
                 pressminus(i+xoffset,j+yoffset,k+zoffset),&
                 avg_det,avg_alp,avg_beta)
          else
            call CCTK_WARN(0, "Flux direction not x,y,z")
          end if
          
!!$        Find minimum and maximum wavespeeds
      
          charmin = min(0.d0, lamplus(1), lamplus(2), lamplus(3), &
               lamplus(4),lamplus(5),  lamminus(1),lamminus(2),lamminus(3),&
               lamminus(4),lamminus(5))  
          
          charmax = max(0.d0, lamplus(1), lamplus(2), lamplus(3), &
               lamplus(4),lamplus(5),  lamminus(1),lamminus(2),lamminus(3),&
               lamminus(4),lamminus(5))
          
          charpm = charmax - charmin
          
!!$        Calculate flux by standard formula
          
          do m = 1,5 
            
            qdiff(m) = consm_i(m) - consp(m)
            
            f1(m) = (charmax * fplus(m) - charmin * fminus(m) + &
                 charmax * charmin * qdiff(m)) / charpm
            
          end do

        end if !!! The end of the SpaceMask check for a trivial RP.
            
        densflux(i, j, k) = f1(1)
        sxflux(i, j, k) = f1(2)
        syflux(i, j, k) = f1(3)
        szflux(i, j, k) = f1(4)
        tauflux(i, j, k) = f1(5)

        if(evolve_Y_e.ne.0) then
           if (densflux(i, j, k) > 0.d0) then
              Y_e_con_flux(i, j, k) = &
                   Y_e_plus(i, j, k) * &
                   densflux(i, j, k)
           else
              Y_e_con_flux(i, j, k) = &
                   Y_e_minus(i + xoffset, j + yoffset, k + zoffset) * &
                   densflux(i, j, k)
           endif
        endif

      end do
    end do
  end do
  !$OMP END PARALLEL DO
      
end subroutine GRHydro_HLLE

 /*@@
   @routine    GRHydro_HLLE_Tracer
   @date       Mon Mar  8 13:47:13 2004
   @author     Ian Hawke
   @desc 
   HLLE just for the tracer.
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/

subroutine GRHydro_HLLE_Tracer(CCTK_ARGUMENTS)

  USE GRHydro_Eigenproblem

  implicit none
  
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS
  
  integer :: i, j, k,  m
  CCTK_REAL, dimension(number_of_tracers) :: consp,consm_i,fplus,fminus,f1
  CCTK_REAL, dimension(5) :: lamminus,lamplus
  CCTK_REAL, dimension(number_of_tracers) :: qdiff
  CCTK_REAL ::  charmin, charmax, charpm,avg_alp,avg_det
  CCTK_REAL :: gxxh, gxyh, gxzh, gyyh, gyzh, gzzh, uxxh, uxyh, &
       uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r
    
  CCTK_INT :: type_bits, trivial, not_trivial

  if (flux_direction == 1) then
    call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemX")
    call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemX", &
         &"trivial")
  else if (flux_direction == 2) then
    call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemY")
    call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemY", &
         &"trivial")
  else if (flux_direction == 3) then
    call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemZ")
    call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemZ", &
         &"trivial")
  else
    call CCTK_WARN(0, "Flux direction not x,y,z")
  end if

  do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil
    do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil
      do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil
        
        f1 = 0.d0
        lamminus = 0.d0
        lamplus = 0.d0
        consp = 0.d0
        consm_i = 0.d0
        fplus = 0.d0
        fminus = 0.d0
        qdiff = 0.d0
        
!!$        Set the left (p for plus) and right (m_i for minus, i+1) states
        
        do m=1,number_of_tracers
           consp(m)   = cons_tracerplus(i,j,k,m)         
           consm_i(m) = cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,m)
        enddo
!!$        Calculate various metric terms here.
!!$        Note also need the average of the lapse at the 
!!$        left and right points.

        if (shift_state .ne. 0) then
          if (flux_direction == 1) then
            avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
                 betax(i,j,k))
          else if (flux_direction == 2) then
            avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
                 betay(i,j,k))
          else if (flux_direction == 3) then
            avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
                 betaz(i,j,k))
          else
            call CCTK_WARN(0, "Flux direction not x,y,z")
          end if
        else
          avg_beta = 0.d0
        end if

        avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset))
        
        gxxh = 0.5d0 * (gxx(i+xoffset,j+yoffset,k+zoffset) + gxx(i,j,k))
        gxyh = 0.5d0 * (gxy(i+xoffset,j+yoffset,k+zoffset) + gxy(i,j,k))
        gxzh = 0.5d0 * (gxz(i+xoffset,j+yoffset,k+zoffset) + gxz(i,j,k))
        gyyh = 0.5d0 * (gyy(i+xoffset,j+yoffset,k+zoffset) + gyy(i,j,k))
        gyzh = 0.5d0 * (gyz(i+xoffset,j+yoffset,k+zoffset) + gyz(i,j,k))
        gzzh = 0.5d0 * (gzz(i+xoffset,j+yoffset,k+zoffset) + gzz(i,j,k))

        avg_det =  SPATIAL_DETERMINANT(gxxh,gxyh,gxzh,\
             gyyh,gyzh,gzzh)

        call UpperMetric(uxxh, uxyh, uxzh, uyyh, uyzh, uzzh, &
             avg_det,gxxh, gxyh, gxzh, &
             gyyh, gyzh, gzzh)
        
        if (flux_direction == 1) then
          usendh = uxxh
        else if (flux_direction == 2) then
          usendh = uyyh
        else if (flux_direction == 3) then
          usendh = uzzh
        else
          call CCTK_WARN(0, "Flux direction not x,y,z")
        end if
          
!!$        Calculate the jumps in the conserved variables
          
        qdiff = consm_i - consp
          
!!$        Eigenvalues and fluxes either side of the cell interface
          
        if (flux_direction == 1) then
          call eigenvalues(GRHydro_eos_handle,&
               rhominus(i+xoffset,j+yoffset,k+zoffset),&
               velxminus(i+xoffset,j+yoffset,k+zoffset),&
               velyminus(i+xoffset,j+yoffset,k+zoffset),&
               velzminus(i+xoffset,j+yoffset,k+zoffset),&
               epsminus(i+xoffset,j+yoffset,k+zoffset),&
               w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),&
               lamminus,gxxh,gxyh,gxzh,gyyh,&
               gyzh,gzzh,&
               usendh,avg_alp,avg_beta)
          call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),&
               velxplus(i,j,k),velyplus(i,j,k),&
               velzplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),&
               lamplus,gxxh,gxyh,gxzh,gyyh,&
               gyzh,gzzh,&
               usendh,avg_alp,avg_beta)
          fplus(:)  = (velxplus(i,j,k)  - avg_beta / avg_alp) * &
               cons_tracerplus(i,j,k,:)
          fminus(:) = (velxminus(i+xoffset,j+yoffset,k+zoffset) - avg_beta / avg_alp) * &
               cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,:)
        else if (flux_direction == 2) then
          call eigenvalues(GRHydro_eos_handle,&
               rhominus(i+xoffset,j+yoffset,k+zoffset),&
               velyminus(i+xoffset,j+yoffset,k+zoffset),&
               velzminus(i+xoffset,j+yoffset,k+zoffset),&
               velxminus(i+xoffset,j+yoffset,k+zoffset),&
               epsminus(i+xoffset,j+yoffset,k+zoffset),&
               w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),&
               lamminus,gyyh,gyzh,gxyh,gzzh,&
               gxzh,gxxh,&
               usendh,avg_alp,avg_beta)
          call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),&
               velyplus(i,j,k),velzplus(i,j,k),&
               velxplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),&
               lamplus,gyyh,gyzh,gxyh,gzzh,&
               gxzh,gxxh,&
               usendh,avg_alp,avg_beta)
          fplus(:)  = (velyplus(i,j,k)  - avg_beta / avg_alp) * &
               cons_tracerplus(i,j,k,:)
          fminus(:) = (velyminus(i+xoffset,j+yoffset,k+zoffset) - avg_beta / avg_alp) * &
               cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,:)
        else if (flux_direction == 3) then
          call eigenvalues(GRHydro_eos_handle,&
               rhominus(i+xoffset,j+yoffset,k+zoffset),&
               velzminus(i+xoffset,j+yoffset,k+zoffset),&
               velxminus(i+xoffset,j+yoffset,k+zoffset),&
               velyminus(i+xoffset,j+yoffset,k+zoffset),&
               epsminus(i+xoffset,j+yoffset,k+zoffset),&
               w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),&
               lamminus,gzzh,gxzh,gyzh,&
               gxxh,gxyh,gyyh,&
               usendh,avg_alp,avg_beta)
          call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),&
               velzplus(i,j,k),velxplus(i,j,k),&
               velyplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),&
               lamplus,gzzh,gxzh,gyzh,&
               gxxh,gxyh,gyyh,&
               usendh,avg_alp,avg_beta)
          fplus(:)  = (velzplus(i,j,k) - avg_beta / avg_alp) * &
               cons_tracerplus(i,j,k,:)
          fminus(:) = (velzminus(i+xoffset,j+yoffset,k+zoffset) - avg_beta / avg_alp) * &
               cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,:)
        else
          call CCTK_WARN(0, "Flux direction not x,y,z")
        end if
        
!!$        Find minimum and maximum wavespeeds
      
        charmin = min(0.d0, lamplus(1), lamplus(2), lamplus(3), &
             lamplus(4),lamplus(5),  lamminus(1),lamminus(2),lamminus(3),&
             lamminus(4),lamminus(5))  
          
        charmax = max(0.d0, lamplus(1), lamplus(2), lamplus(3), &
             lamplus(4),lamplus(5),  lamminus(1),lamminus(2),lamminus(3),&
             lamminus(4),lamminus(5))
        
        charpm = charmax - charmin
          
!!$        Calculate flux by standard formula
          
        do m = 1,number_of_tracers
            
          qdiff(m) = consm_i(m) - consp(m)
          
          f1(m) = (charmax * fplus(m) - charmin * fminus(m) + &
               charmax * charmin * qdiff(m)) / charpm
          
        end do
            
        cons_tracerflux(i, j, k,:) = f1(:)
!!$
!!$        if ( ((flux_direction.eq.3).and.(i.eq.4).and.(j.eq.4)).or.&
!!$             ((flux_direction.eq.2).and.(i.eq.4).and.(k.eq.4)).or.&
!!$             ((flux_direction.eq.1).and.(j.eq.4).and.(k.eq.4))&
!!$             ) then
!!$          write(*,*) flux_direction, i, j, k, f1(1), consm_i(1), consp(1)
!!$        end if
        
      end do
    end do
  end do

      
end subroutine GRHydro_HLLE_Tracer