file/parameter string replacement from whisky to GRHydro

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@112 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45

1 files changed, 576 insertions, 0 deletions

diff --git a/src/GRHydro_HLLE.F90 b/src/GRHydro_HLLE.F90
new file mode 100644
index 0000000..d32f4a7
--- /dev/null
+++ b/src/GRHydro_HLLE.F90
@@ -0,0 +1,576 @@
+ /*@@
+   @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 "SpaceMask.h"
+
+ /*@@
+   @routine    GRHydro_HLLEGeneral
+   @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.
+   Rewritten for the new EOS interface.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+   Altered from Cactus 3 routines originally written by Toni Font.
+   @endhistory 
+
+@@*/
+
+subroutine GRHydro_HLLEGeneral(CCTK_ARGUMENTS)
+
+  USE GRHydro_Eigenproblem
+  implicit none
+  
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_FUNCTIONS
+  
+  integer :: i, j, k, m
+  CCTK_REAL, dimension(5) :: cons_p,cons_m,fplus,fminus,lamplus
+  CCTK_REAL, dimension(5) :: f1,lamminus
+  CCTK_REAL, dimension(5) :: qdiff
+  CCTK_REAL, dimension(6) :: prim_p, prim_m
+  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, psi4h, &
+       cs2_p, cs2_m, dpdeps_p, dpdeps_m
+    
+  CCTK_INT :: type_bits, trivial, not_trivial
+
+  integer tadmor
+
+  if(CCTK_EQUALS(HLLE_type,"Tadmor")) then
+     tadmor = 1
+  else
+     tadmor = 0
+  endif
+
+
+  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
+        cons_p = 0.d0
+        cons_m = 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
+        
+        cons_p(1)   = densplus(i,j,k) 
+        cons_p(2)   = sxplus(i,j,k)
+        cons_p(3)   = syplus(i,j,k)
+        cons_p(4)   = szplus(i,j,k)
+        cons_p(5)   = tauplus(i,j,k)
+        
+        cons_m(1) = densminus(i+xoffset,j+yoffset,k+zoffset)
+        cons_m(2) = sxminus(i+xoffset,j+yoffset,k+zoffset)
+        cons_m(3) = syminus(i+xoffset,j+yoffset,k+zoffset)
+        cons_m(4) = szminus(i+xoffset,j+yoffset,k+zoffset)
+        cons_m(5) = tauminus(i+xoffset,j+yoffset,k+zoffset) 
+        
+        prim_p(1)   = rhoplus(i,j,k) 
+        prim_p(2)   = velxplus(i,j,k)
+        prim_p(3)   = velyplus(i,j,k) 
+        prim_p(4)   = velzplus(i,j,k)
+        prim_p(5)   = epsplus(i,j,k)
+        
+        prim_m(1) = rhominus(i+xoffset,j+yoffset,k+zoffset)
+        prim_m(2) = velxminus(i+xoffset,j+yoffset,k+zoffset)
+        prim_m(3) = velyminus(i+xoffset,j+yoffset,k+zoffset)
+        prim_m(4) = velzminus(i+xoffset,j+yoffset,k+zoffset)
+        prim_m(5) = epsminus(i+xoffset,j+yoffset,k+zoffset) 
+        
+        prim_p(6) = pressplus(i,j,k)
+        cs2_p     = eos_cs2_p(i,j,k)
+        dpdeps_p  = eos_dpdeps_p(i,j,k)
+
+        prim_m(6) = pressminus(i+xoffset,j+yoffset,k+zoffset)
+        cs2_m     = eos_cs2_m(i+xoffset,j+yoffset,k+zoffset)
+        dpdeps_m  = eos_dpdeps_m(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))
+
+        if (conformal_state .eq. 0) then
+          psi4h = 1.d0
+        else
+          psi4h = (0.5d0*(psi(i,j,k)+psi(i+xoffset,j+yoffset,k+zoffset)))**4
+        end if
+          
+        call SpatialDeterminant(psi4h*gxxh,psi4h*gxyh,psi4h*gxzh,&
+             psi4h*gyyh,psi4h*gyzh,psi4h*gzzh,avg_det)
+
+!!$ 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(cons_m(1),cons_m(2),cons_m(3),cons_m(4),cons_m(5),&
+                 f1(1),f1(2),f1(3),f1(4),f1(5),&
+                 prim_m(2),prim_m(6), &
+                 avg_det,avg_alp,avg_beta)
+          else if (flux_direction == 2) then
+            call num_x_flux(cons_m(1),cons_m(3),cons_m(4),cons_m(2),cons_m(5),&
+                 f1(1),f1(3),f1(4),f1(2),f1(5),&
+                 prim_m(3),prim_m(6), &
+                 avg_det,avg_alp,avg_beta)
+          else if (flux_direction == 3) then
+            call num_x_flux(cons_m(1),cons_m(4),cons_m(2),cons_m(3),cons_m(5),&
+                 f1(1),f1(4),f1(2),f1(3),f1(5), &
+                 prim_m(4),prim_m(6), &
+                 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,psi4h*gxxh, psi4h*gxyh, psi4h*gxzh, &
+               psi4h*gyyh, psi4h*gyzh, psi4h*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) = cons_m(1) - cons_p(1)
+          qdiff(2) = cons_m(2) - cons_p(2)
+          qdiff(3) = cons_m(3) - cons_p(3)
+          qdiff(4) = cons_m(4) - cons_p(4)
+          qdiff(5) = cons_m(5) - cons_p(5)
+          
+!!$        Eigenvalues and fluxes either side of the cell interface
+          
+          if (flux_direction == 1) then
+            call eigenvalues_general(&
+                 prim_m(2),prim_m(3),prim_m(4),cs2_m, &
+                 lamminus,&
+                 psi4h*gxxh,psi4h*gxyh,psi4h*gxzh,&
+                 psi4h*gyyh,psi4h*gyzh,psi4h*gzzh,&
+                 usendh,avg_alp,avg_beta)
+            call eigenvalues_general(&
+                 prim_p(2),prim_p(3),prim_p(4),cs2_p, &
+                 lamplus,&
+                 psi4h*gxxh,psi4h*gxyh,psi4h*gxzh,&
+                 psi4h*gyyh,psi4h*gyzh,psi4h*gzzh,&
+                 usendh,avg_alp,avg_beta)
+            call num_x_flux(cons_p(1),cons_p(2),cons_p(3),cons_p(4),cons_p(5),&
+                 fplus(1),fplus(2),fplus(3),fplus(4),fplus(5),&
+                 prim_p(2),prim_p(6), &
+                 avg_det,avg_alp,avg_beta)
+            call num_x_flux(cons_m(1),cons_m(2),cons_m(3),cons_m(4),cons_m(5),&
+                 fminus(1),fminus(2),fminus(3),fminus(4),fminus(5),&
+                 prim_m(2),prim_m(6), &
+                 avg_det,avg_alp,avg_beta)
+          else if (flux_direction == 2) then
+            call eigenvalues_general(&
+                 prim_m(3),prim_m(4),prim_m(2),cs2_m, &
+                 lamminus,&
+                 psi4h*gyyh,psi4h*gyzh,psi4h*gxyh,&
+                 psi4h*gzzh,psi4h*gxzh,psi4h*gxxh,&
+                 usendh,avg_alp,avg_beta)
+            call eigenvalues_general(&
+                 prim_p(3),prim_p(4),prim_p(2),cs2_p, &
+                 lamplus,&
+                 psi4h*gyyh,psi4h*gyzh,psi4h*gxyh,&
+                 psi4h*gzzh,psi4h*gxzh,psi4h*gxxh,&
+                 usendh,avg_alp,avg_beta)
+            call num_x_flux(cons_p(1),cons_p(3),cons_p(4),cons_p(2),cons_p(5),&
+                 fplus(1),fplus(3),fplus(4),fplus(2),fplus(5),&
+                 prim_p(3),prim_p(6), &
+                 avg_det,avg_alp,avg_beta)
+            call num_x_flux(cons_m(1),cons_m(3),cons_m(4),cons_m(2),cons_m(5),&
+                 fminus(1),fminus(3),fminus(4),fminus(2),fminus(5),&
+                 prim_m(3),prim_m(6), &
+                 avg_det,avg_alp,avg_beta)
+          else if (flux_direction == 3) then
+            call eigenvalues_general(&
+                 prim_m(4),prim_m(2),prim_m(3),cs2_m, &
+                 lamminus,&
+                 psi4h*gzzh,psi4h*gxzh,psi4h*gyzh,&
+                 psi4h*gxxh,psi4h*gxyh,psi4h*gyyh,&
+                 usendh,avg_alp,avg_beta)
+            call eigenvalues_general(&
+                 prim_p(4),prim_p(2),prim_p(3),cs2_p, &
+                 lamplus,&
+                 psi4h*gzzh,psi4h*gxzh,psi4h*gyzh,&
+                 psi4h*gxxh,psi4h*gxyh,psi4h*gyyh,&
+                 usendh,avg_alp,avg_beta)
+            call num_x_flux(cons_p(1),cons_p(4),cons_p(2),cons_p(3),cons_p(5),&
+                 fplus(1),fplus(4),fplus(2),fplus(3),fplus(5), &
+                 prim_p(4),prim_p(6), &
+                 avg_det,avg_alp,avg_beta)
+            call num_x_flux(cons_m(1),cons_m(4),cons_m(2),cons_m(3),cons_m(5),&
+                 fminus(1),fminus(4),fminus(2),fminus(3),fminus(5), &
+                 prim_m(4),prim_m(6), &
+                 avg_det,avg_alp,avg_beta)
+          else
+            call CCTK_WARN(0, "Flux direction not x,y,z")
+          end if
+
+          if(tadmor.eq.0) then
+
+!!$        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) = cons_m(m) - cons_p(m)
+               
+               f1(m) = (charmax * fplus(m) - charmin * fminus(m) + &
+                    charmax * charmin * qdiff(m)) / charpm
+               
+            end do
+            
+         else
+            ! Tadmor's semi-discrete scheme: JcP 160, 241 (2000)
+
+            charmax = max(abs(lamplus(1)), abs(lamplus(2)), abs(lamplus(3)), &
+                 abs(lamplus(4)),abs(lamplus(5)),abs(lamminus(1)),abs(lamminus(2)), &
+                 abs(lamminus(3)),abs(lamminus(4)),abs(lamminus(5)))
+
+            do m = 1,5 
+            
+               qdiff(m) = cons_m(m) - cons_p(m)
+               
+               f1(m) = 0.5d0 * (fplus(m) + fminus(m)) - 0.5d0*charmax* &
+                    qdiff(m)
+               
+            end do
+            
+
+         end if
+
+
+
+        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)
+
+      end do
+    end do
+  end do
+      
+end subroutine GRHydro_HLLEGeneral
+
+
+subroutine GRHydro_HLLE_TracerGeneral(CCTK_ARGUMENTS)
+
+  USE GRHydro_Eigenproblem
+
+  implicit none
+  
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_FUNCTIONS
+  
+  integer :: i, j, k, m
+  CCTK_REAL, dimension(number_of_tracers) :: cons_p,cons_m,fplus,fminus,f1
+  CCTK_REAL, dimension(5) :: lamminus,lamplus
+  CCTK_REAL, dimension(number_of_tracers) :: qdiff
+  CCTK_REAL, dimension(6) :: prim_p, prim_m
+  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, psi4h, &
+       cs2_p, cs2_m, dpdeps_p, dpdeps_m
+
+  integer tadmor
+
+  if(CCTK_EQUALS(HLLE_type,"Tadmor")) then
+     tadmor = 1
+  else
+     tadmor = 0
+  endif
+    
+
+  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
+        cons_p = 0.d0
+        cons_m = 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
+        
+        cons_p(:) = cons_tracerplus(i,j,k,:) 
+        cons_m(:) = cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,:)
+        
+        prim_p(1)   = rhoplus(i,j,k) 
+        prim_p(2)   = velxplus(i,j,k)
+        prim_p(3)   = velyplus(i,j,k) 
+        prim_p(4)   = velzplus(i,j,k)
+        prim_p(5)   = epsplus(i,j,k)
+        
+        prim_m(1) = rhominus(i+xoffset,j+yoffset,k+zoffset)
+        prim_m(2) = velxminus(i+xoffset,j+yoffset,k+zoffset)
+        prim_m(3) = velyminus(i+xoffset,j+yoffset,k+zoffset)
+        prim_m(4) = velzminus(i+xoffset,j+yoffset,k+zoffset)
+        prim_m(5) = epsminus(i+xoffset,j+yoffset,k+zoffset) 
+        
+        prim_p(6) = pressplus(i,j,k)
+        cs2_p     = eos_cs2_p(i,j,k)
+        dpdeps_p  = eos_dpdeps_p(i,j,k)
+
+        prim_m(6) = pressminus(i+xoffset,j+yoffset,k+zoffset)
+        cs2_m     = eos_cs2_m(i+xoffset,j+yoffset,k+zoffset)
+        dpdeps_m  = eos_dpdeps_m(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))
+
+        if (conformal_state .eq. 0) then
+          psi4h = 1.d0
+        else
+          psi4h = (0.5d0*(psi(i,j,k)+psi(i+xoffset,j+yoffset,k+zoffset)))**4
+        end if
+          
+        call SpatialDeterminant(psi4h*gxxh,psi4h*gxyh,psi4h*gxzh,&
+             psi4h*gyyh,psi4h*gyzh,psi4h*gzzh,avg_det)
+
+!!$ If the Riemann problem is trivial, just calculate the fluxes from the 
+!!$ left state and skip to the next cell
+          
+        call UpperMetric(uxxh, uxyh, uxzh, uyyh, uyzh, uzzh, &
+             avg_det,psi4h*gxxh, psi4h*gxyh, psi4h*gxzh, &
+             psi4h*gyyh, psi4h*gyzh, psi4h*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
+          
+!!$        Eigenvalues and fluxes either side of the cell interface
+          
+          if (flux_direction == 1) then
+            call eigenvalues_general(&
+                 prim_m(2),prim_m(3),prim_m(4),cs2_m, &
+                 lamminus,&
+                 psi4h*gxxh,psi4h*gxyh,psi4h*gxzh,&
+                 psi4h*gyyh,psi4h*gyzh,psi4h*gzzh,&
+                 usendh,avg_alp,avg_beta)
+            call eigenvalues_general(&
+                 prim_p(2),prim_p(3),prim_p(4),cs2_p, &
+                 lamplus,&
+                 psi4h*gxxh,psi4h*gxyh,psi4h*gxzh,&
+                 psi4h*gyyh,psi4h*gyzh,psi4h*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_general(&
+                 prim_m(3),prim_m(4),prim_m(2),cs2_m, &
+                 lamminus,&
+                 psi4h*gyyh,psi4h*gyzh,psi4h*gxyh,&
+                 psi4h*gzzh,psi4h*gxzh,psi4h*gxxh,&
+                 usendh,avg_alp,avg_beta)
+            call eigenvalues_general(&
+                 prim_p(3),prim_p(4),prim_p(2),cs2_p, &
+                 lamplus,&
+                 psi4h*gyyh,psi4h*gyzh,psi4h*gxyh,&
+                 psi4h*gzzh,psi4h*gxzh,psi4h*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_general(&
+                 prim_m(4),prim_m(2),prim_m(3),cs2_m, &
+                 lamminus,&
+                 psi4h*gzzh,psi4h*gxzh,psi4h*gyzh,&
+                 psi4h*gxxh,psi4h*gxyh,psi4h*gyyh,&
+                 usendh,avg_alp,avg_beta)
+            call eigenvalues_general(&
+                 prim_p(4),prim_p(2),prim_p(3),cs2_p, &
+                 lamplus,&
+                 psi4h*gzzh,psi4h*gxzh,psi4h*gyzh,&
+                 psi4h*gxxh,psi4h*gxyh,psi4h*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
+
+          if(tadmor.eq.0) then
+          
+!!$        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) = cons_m(m) - cons_p(m)
+            
+                f1(m) = (charmax * fplus(m) - charmin * fminus(m) + &
+                     charmax * charmin * qdiff(m)) / charpm
+             end do
+
+         else
+            ! Tadmor's semi-descrite scheme: JcP 160, 241 (2000)
+
+            charmax = max(abs(lamplus(1)), abs(lamplus(2)), abs(lamplus(3)), &
+                 abs(lamplus(4)),abs(lamplus(5)),abs(lamminus(1)),abs(lamminus(2)), &
+                 abs(lamminus(3)),abs(lamminus(4)),abs(lamminus(5)))
+
+            do m = 1,number_of_tracers 
+            
+               qdiff(m) = cons_m(m) - cons_p(m)
+               
+               f1(m) = 0.5d0 * (fplus(m) + fminus(m)) - 0.5d0*charmax* &
+                    qdiff(m)
+               
+            end do
+            
+
+         end if
+
+            
+         cons_tracerflux(i,j,k,:) = f1(:)
+
+      end do
+    end do
+  end do
+
+end subroutine GRHydro_HLLE_TracerGeneral
+