* rename HLLE routine source code files to reflect

  true content



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

1 files changed, 357 insertions, 312 deletions

diff --git a/src/GRHydro_HLLE.F90 b/src/GRHydro_HLLE.F90
index 7a3f5f8..bd976e1 100644
--- a/src/GRHydro_HLLE.F90
+++ b/src/GRHydro_HLLE.F90
@@ -12,16 +12,16 @@
 #include "cctk_Parameters.h"
 #include "cctk_Arguments.h"
 #include "cctk_Functions.h"
+
 #include "GRHydro_Macros.h"
 #include "SpaceMask.h"
 
  /*@@
-   @routine    GRHydro_HLLEGeneral
+   @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.
-   Rewritten for the new EOS interface.
    @enddesc 
    @calls     
    @calledby   
@@ -31,36 +31,25 @@
 
 @@*/
 
-subroutine GRHydro_HLLEGeneral(CCTK_ARGUMENTS)
-
+subroutine GRHydro_HLLE(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) :: consp,consm_i,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, &
-       cs2_p, cs2_m, dpdeps_p, dpdeps_m
+       uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r
+  CCTK_REAL :: xtemp
     
   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", &
@@ -84,45 +73,25 @@ subroutine GRHydro_HLLEGeneral(CCTK_ARGUMENTS)
         f1 = 0.d0
         lamminus = 0.d0
         lamplus = 0.d0
-        cons_p = 0.d0
-        cons_m = 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
         
-        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)
+        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)
         
-        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)
+        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 
@@ -154,27 +123,31 @@ subroutine GRHydro_HLLEGeneral(CCTK_ARGUMENTS)
         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)
+        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(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), &
+            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(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), &
+            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(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), &
+            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")
@@ -183,8 +156,8 @@ subroutine GRHydro_HLLEGeneral(CCTK_ARGUMENTS)
         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)
+               avg_det,gxxh, gxyh, gxzh, &
+               gyyh, gyzh, gzzh)
           
           if (flux_direction == 1) then
             usendh = uxxh
@@ -198,127 +171,201 @@ subroutine GRHydro_HLLEGeneral(CCTK_ARGUMENTS)
           
 !!$        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)
+          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
-            call eigenvalues_general(&
-                 prim_m(2),prim_m(3),prim_m(4),cs2_m, &
-                 lamminus,&
-                 gxxh,gxyh,gxzh,&
-                 gyyh,gyzh,gzzh,&
-                 usendh,avg_alp,avg_beta)
-            call eigenvalues_general(&
-                 prim_p(2),prim_p(3),prim_p(4),cs2_p, &
-                 lamplus,&
-                 gxxh,gxyh,gxzh,&
-                 gyyh,gyzh,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), &
+             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(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), &
+            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
-            call eigenvalues_general(&
-                 prim_m(3),prim_m(4),prim_m(2),cs2_m, &
-                 lamminus,&
-                 gyyh,gyzh,gxyh,&
-                 gzzh,gxzh,gxxh,&
-                 usendh,avg_alp,avg_beta)
-            call eigenvalues_general(&
-                 prim_p(3),prim_p(4),prim_p(2),cs2_p, &
-                 lamplus,&
-                 gyyh,gyzh,gxyh,&
-                 gzzh,gxzh,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), &
+             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(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), &
+            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
-            call eigenvalues_general(&
-                 prim_m(4),prim_m(2),prim_m(3),cs2_m, &
-                 lamminus,&
-                 gzzh,gxzh,gyzh,&
-                 gxxh,gxyh,gyyh,&
-                 usendh,avg_alp,avg_beta)
-            call eigenvalues_general(&
-                 prim_p(4),prim_p(2),prim_p(3),cs2_p, &
-                 lamplus,&
-                 gzzh,gxzh,gyzh,&
-                 gxxh,gxyh,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), &
+             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
-
-          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))  
+          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))
+          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
+          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
+          do m = 1,5 
             
-         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) = consm_i(m) - consp(m)
             
-               qdiff(m) = cons_m(m) - cons_p(m)
-               
-               f1(m) = 0.5d0 * (fplus(m) + fminus(m)) - 0.5d0*charmax* &
-                    qdiff(m)
-               
-            end do
+            f1(m) = (charmax * fplus(m) - charmin * fminus(m) + &
+                 charmax * charmin * qdiff(m)) / charpm
             
-         end if
-
+          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)
@@ -330,22 +377,36 @@ subroutine GRHydro_HLLEGeneral(CCTK_ARGUMENTS)
               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
+           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
   end do
+
       
+end subroutine GRHydro_HLLE
 
-end subroutine GRHydro_HLLEGeneral
+ /*@@
+   @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_TracerGeneral(CCTK_ARGUMENTS)
+subroutine GRHydro_HLLE_Tracer(CCTK_ARGUMENTS)
 
   USE GRHydro_Eigenproblem
 
@@ -353,26 +414,32 @@ subroutine GRHydro_HLLE_TracerGeneral(CCTK_ARGUMENTS)
   
   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
+  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, 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, &
-       cs2_p, cs2_m, dpdeps_p, dpdeps_m
-
-  integer tadmor
+       uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r
+    
+  CCTK_INT :: type_bits, trivial, not_trivial
 
-  if(CCTK_EQUALS(HLLE_type,"Tadmor")) then
-     tadmor = 1
+  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
-     tadmor = 0
-  endif
-    
+    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
@@ -381,37 +448,18 @@ subroutine GRHydro_HLLE_TracerGeneral(CCTK_ARGUMENTS)
         f1 = 0.d0
         lamminus = 0.d0
         lamplus = 0.d0
-        cons_p = 0.d0
-        cons_m = 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
         
-        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)
-          
+        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.
@@ -442,135 +490,132 @@ subroutine GRHydro_HLLE_TracerGeneral(CCTK_ARGUMENTS)
         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)
+        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
-          
         call UpperMetric(uxxh, uxyh, uxzh, uyyh, uyzh, uzzh, &
              avg_det,gxxh, gxyh, gxzh, &
              gyyh, gyzh, gzzh)
-          
+        
         if (flux_direction == 1) then
-           usendh = uxxh
+          usendh = uxxh
         else if (flux_direction == 2) then
-           usendh = uyyh
+          usendh = uyyh
         else if (flux_direction == 3) then
-           usendh = uzzh
+          usendh = uzzh
         else
-           call CCTK_WARN(0, "Flux direction not x,y,z")
+          call CCTK_WARN(0, "Flux direction not x,y,z")
         end if
           
-!!$        Eigenvalues and fluxes either side of the cell interface
+!!$        Calculate the jumps in the conserved variables
           
-          if (flux_direction == 1) then
-            call eigenvalues_general(&
-                 prim_m(2),prim_m(3),prim_m(4),cs2_m, &
-                 lamminus,&
-                 gxxh,gxyh,gxzh,&
-                 gyyh,gyzh,gzzh,&
-                 usendh,avg_alp,avg_beta)
-            call eigenvalues_general(&
-                 prim_p(2),prim_p(3),prim_p(4),cs2_p, &
-                 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_general(&
-                 prim_m(3),prim_m(4),prim_m(2),cs2_m, &
-                 lamminus,&
-                 gyyh,gyzh,gxyh,&
-                 gzzh,gxzh,gxxh,&
-                 usendh,avg_alp,avg_beta)
-            call eigenvalues_general(&
-                 prim_p(3),prim_p(4),prim_p(2),cs2_p, &
-                 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_general(&
-                 prim_m(4),prim_m(2),prim_m(3),cs2_m, &
-                 lamminus,&
-                 gzzh,gxzh,gyzh,&
-                 gxxh,gxyh,gyyh,&
-                 usendh,avg_alp,avg_beta)
-            call eigenvalues_general(&
-                 prim_p(4),prim_p(2),prim_p(3),cs2_p, &
-                 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
-
-          if(tadmor.eq.0) then
+        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))  
+        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
+        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
+        do m = 1,number_of_tracers
             
-
-         end if
-
+          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(:)
-
+        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_TracerGeneral
-
-
-
+      
+end subroutine GRHydro_HLLE_Tracer