Current RIT GRMHD code contributions:

Add the magnetized counterparts for several GRHydro routines.

Adjust interface.ccl, param.ccl and schedule.ccl appropriately.



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

1 files changed, 642 insertions, 0 deletions

diff --git a/src/GRHydro_HLLEM.F90 b/src/GRHydro_HLLEM.F90
new file mode 100644
index 0000000..98c4874
--- /dev/null
+++ b/src/GRHydro_HLLEM.F90
@@ -0,0 +1,642 @@
+ /*@@
+   @file      GRHydro_HLLEM.F90
+   @date      Aug 30, 2010
+   @author    Joshua Faber, Scott Noble, Bruno Mundim, 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_HLLEGeneralM
+   @date       Aug 30, 2010
+   @author     Joshua Faber, Scott Noble, Bruno Mundim, 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_HLLEGeneralM(CCTK_ARGUMENTS)
+  USE GRHydro_EigenproblemM
+
+  implicit none
+  
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_FUNCTIONS
+  
+  integer :: i, j, k, m
+  CCTK_REAL, dimension(8) :: cons_p,cons_m,fplus,fminus,f1,qdiff
+  CCTK_REAL, dimension(6) :: prim_p, prim_m
+  CCTK_REAL, dimension(5) :: lamminus,lamplus
+  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
+  CCTK_REAL :: rhoenth_p, rhoenth_m, avg_betax, avg_betay, avg_betaz
+  CCTK_REAL :: vxtp,vytp,vztp,vxtm,vytm,vztm,ab0p,ab0m,b2p,b2m,bdotvp,bdotvm
+  CCTK_REAL :: wp,wm,v2p,v2m,bxlowp,bxlowm,bylowp,bylowm,bzlowp,bzlowm,vA2m,vA2p
+  CCTK_REAL :: Bvecxlowp,Bvecxlowm,Bvecylowp,Bvecylowm,Bveczlowp,Bveczlowm
+  CCTK_REAL :: pressstarp,pressstarm,velxlowp,velxlowm,velylowp,velylowm,velzlowp,velzlowm
+    
+  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_p(6)   = Bvecxplus(i,j,k)
+        cons_p(7)   = Bvecyplus(i,j,k)
+        cons_p(8)   = Bveczplus(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) 
+        cons_m(6) = Bvecxminus(i+xoffset,j+yoffset,k+zoffset) 
+        cons_m(7) = Bvecyminus(i+xoffset,j+yoffset,k+zoffset) 
+        cons_m(8) = Bveczminus(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)
+        rhoenth_p    = prim_p(1)*(1.0d0+prim_p(5))+prim_p(6)
+
+        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)
+        rhoenth_m    = prim_m(1)*(1.0d0+prim_m(5))+prim_m(6)
+          
+!!$        Calculate various metric terms here.
+!!$        Note also need the average of the lapse at the 
+!!$        left and right points.
+!!$        
+!!$        In MHD, we need all three shift components regardless of the flux direction
+
+        if (shift_state .ne. 0) then
+           avg_betax = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
+                betax(i,j,k))
+           avg_betay = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
+                betay(i,j,k))
+           avg_betaz = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
+                betaz(i,j,k))
+           if (flux_direction == 1) then
+              avg_beta=avg_betax
+           else if (flux_direction == 2) then
+              avg_beta=avg_betay
+           else if (flux_direction == 3) then
+              avg_beta=avg_betaz
+           else
+              call CCTK_WARN(0, "Flux direction not x,y,z")
+           end if
+        else
+           avg_beta = 0.d0
+           avg_betax = 0.d0
+           avg_betay = 0.d0
+           avg_betaz = 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)
+
+        vxtp = prim_p(2)-avg_betax/avg_alp
+        vytp = prim_p(3)-avg_betay/avg_alp
+        vztp = prim_p(4)-avg_betaz/avg_alp
+        vxtm = prim_m(2)-avg_betax/avg_alp
+        vytm = prim_m(3)-avg_betay/avg_alp
+        vztm = prim_m(4)-avg_betaz/avg_alp
+
+        call calc_vlow_blow(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, &
+             prim_p(2),prim_p(3),prim_p(4),cons_p(6),cons_p(7),cons_p(8), &
+             velxlowp,velylowp,velzlowp,Bvecxlowp,Bvecylowp,Bveczlowp, &
+             bdotvp,b2p,v2p,wp,bxlowp,bylowp,bzlowp)
+        call calc_vlow_blow(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, &
+             prim_m(2),prim_m(3),prim_m(4),cons_m(6),cons_m(7),cons_m(8), &
+             velxlowm,velylowm,velzlowm,Bvecxlowm,Bvecylowm,Bveczlowm, &
+             bdotvm,b2m,v2m,wm,bxlowm,bylowm,bzlowm)
+        
+        ab0p = wp*bdotvp
+        ab0m = wm*bdotvm
+
+        vA2p = b2p/(rhoenth_p+b2p)
+        vA2m = b2m/(rhoenth_m+b2m)
+
+!!$ p^*  = p+b^2/2  in Anton et al.
+        pressstarp = prim_p(6)+0.5d0*b2p
+        pressstarm = prim_m(6)+0.5d0*b2m        
+
+
+!!$ 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
+          
+!!$ we now pass in the B-field as conserved and flux, the vtilde's instead of v's,
+!!$ pressstar instead of P, b_i, alp b^0, w, metric determinant, 
+!!$ alp, and beta in the flux dir
+
+          if (flux_direction == 1) then
+            call num_x_fluxM(cons_m(1),cons_m(2),cons_m(3),cons_m(4),cons_m(5),&
+                 cons_m(6),cons_m(7),cons_m(8),&
+                 f1(1),f1(2),f1(3),f1(4),f1(5),f1(6),f1(7),f1(8),&
+                 vxtm,vytm,vztm,pressstarm,bxlowm,bylowm,bzlowm,ab0m,wm, &
+                 avg_det,avg_alp,avg_beta)
+          else if (flux_direction == 2) then
+            call num_x_fluxM(cons_m(1),cons_m(3),cons_m(4),cons_m(2),cons_m(5),&
+                 cons_m(7),cons_m(8),cons_m(6),&
+                 f1(1),f1(3),f1(4),f1(2),f1(5),f1(7),f1(8),f1(6),&
+                 vytm,vztm,vxtm,pressstarm,bylowm,bzlowm,bxlowm,ab0m,wm, &
+                 avg_det,avg_alp,avg_beta)
+          else if (flux_direction == 3) then
+            call num_x_fluxM(cons_m(1),cons_m(4),cons_m(2),cons_m(3),cons_m(5),&
+                 cons_m(8),cons_m(6),cons_m(7),&
+                 f1(1),f1(4),f1(2),f1(3),f1(5),f1(8),f1(6),f1(7), &
+                 vztm,vxtm,vytm,pressstarm,bzlowm,bxlowm,bylowm,ab0m,wm, &
+                 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) = 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(6) = cons_m(6) - cons_p(6)
+          qdiff(7) = cons_m(7) - cons_p(7)
+          qdiff(8) = cons_m(8) - cons_p(8)
+          
+!!$        Eigenvalues and fluxes either side of the cell interface
+          
+          if (flux_direction == 1) then
+            call eigenvalues_generalM(&
+                 prim_m(2),prim_m(3),prim_m(4),cs2_m,vA2m, &
+                 lamminus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,&
+                 usendh,avg_alp,avg_beta)
+            call eigenvalues_generalM(&
+                 prim_p(2),prim_p(3),prim_p(4),cs2_p,vA2p, &
+                 lamplus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,&
+                 usendh,avg_alp,avg_beta)
+            call num_x_fluxM(cons_p(1),cons_p(2),cons_p(3),cons_p(4),cons_p(5),&
+                 cons_p(6),cons_p(7),cons_p(8),&
+                 fplus(1),fplus(2),fplus(3),fplus(4),fplus(5),fplus(6),fplus(7),fplus(8),&
+                 vxtp,vytp,vztp,pressstarp,bxlowp,bylowp,bzlowp,ab0p,wp, &
+                 avg_det,avg_alp,avg_beta)
+            call num_x_fluxM(cons_m(1),cons_m(2),cons_m(3),cons_m(4),cons_m(5),&
+                 cons_m(6),cons_m(7),cons_m(8),&
+                 fminus(1),fminus(2),fminus(3),fminus(4),fminus(5),&
+                 fminus(6),fminus(7),fminus(8),&
+                 vxtm,vytm,vztm,pressstarm,bxlowm,bylowm,bzlowm,ab0m,wm, &
+                 avg_det,avg_alp,avg_beta)
+          else if (flux_direction == 2) then
+            call eigenvalues_generalM(&
+                 prim_m(3),prim_m(4),prim_m(2),cs2_m,vA2m, &
+                 lamminus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,&
+                 usendh,avg_alp,avg_beta)
+            call eigenvalues_generalM(&
+                 prim_p(3),prim_p(4),prim_p(2),cs2_p,vA2p, &
+                 lamplus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,&
+                 usendh,avg_alp,avg_beta)
+            call num_x_fluxM(cons_p(1),cons_p(3),cons_p(4),cons_p(2),cons_p(5),&
+                 cons_p(7),cons_p(8),cons_p(6),&
+                 fplus(1),fplus(3),fplus(4),fplus(2),fplus(5),fplus(7),fplus(8),fplus(6),&
+                 vytp,vztp,vxtp,pressstarp,bylowp,bzlowp,bxlowp,ab0p,wp, &
+                 avg_det,avg_alp,avg_beta)
+            call num_x_fluxM(cons_m(1),cons_m(3),cons_m(4),cons_m(2),cons_m(5),&
+                 cons_m(7),cons_m(8),cons_m(6),&
+                 fminus(1),fminus(3),fminus(4),fminus(2),fminus(5),&
+                 fminus(7),fminus(8),fminus(6),&
+                 vytm,vztm,vxtm,pressstarm,bylowm,bzlowm,bxlowm,ab0m,wm, &
+                 avg_det,avg_alp,avg_beta)
+          else if (flux_direction == 3) then
+            call eigenvalues_generalM(&
+                 prim_m(4),prim_m(2),prim_m(3),cs2_m,vA2m, &
+                 lamminus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,&
+                 usendh,avg_alp,avg_beta)
+            call eigenvalues_generalM(&
+                 prim_p(4),prim_p(2),prim_p(3),cs2_p,vA2p, &
+                 lamplus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,&
+                 usendh,avg_alp,avg_beta)
+            call num_x_fluxM(cons_p(1),cons_p(4),cons_p(2),cons_p(3),cons_p(5),&
+                 cons_p(8),cons_p(6),cons_p(7),&
+                 fplus(1),fplus(4),fplus(2),fplus(3),fplus(5),fplus(8),fplus(6),fplus(7), &
+                 vztp,vxtp,vytp,pressstarp,bzlowp,bxlowp,bylowp,ab0p,wp, &
+                 avg_det,avg_alp,avg_beta)
+            call num_x_fluxM(cons_m(1),cons_m(4),cons_m(2),cons_m(3),cons_m(5),&
+                 cons_m(8),cons_m(6),cons_m(7),&
+                 fminus(1),fminus(4),fminus(2),fminus(3),fminus(5), &
+                 fminus(8),fminus(6),fminus(7), &
+                 vztm,vxtm,vytm,pressstarm,bzlowm,bxlowm,bylowm,ab0m,wm, &
+                 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,8 
+            
+               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,8 
+            
+               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)
+        Bvecxflux(i, j, k) = f1(6)
+        Bvecyflux(i, j, k) = f1(7)
+        Bveczflux(i, j, k) = f1(8)
+
+      end do
+    end do
+  end do
+      
+end subroutine GRHydro_HLLEGeneralM
+
+
+subroutine GRHydro_HLLE_TracerGeneralM(CCTK_ARGUMENTS)
+
+  USE GRHydro_EigenproblemM
+
+  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, dimension(3) :: mag_p, mag_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
+  CCTK_REAL :: b2p,b2m,vA2m,vA2p
+
+  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
+        mag_p = 0.d0
+        mag_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,:)
+        
+        mag_p(1) = Bvecxplus(i,j,k)
+        mag_p(2) = Bvecyplus(i,j,k)
+        mag_p(3) = Bveczplus(i,j,k)
+
+        mag_m(1) = Bvecxminus(i+xoffset,j+yoffset,k+zoffset)
+        mag_m(2) = Bvecyminus(i+xoffset,j+yoffset,k+zoffset)
+        mag_m(3) = Bveczminus(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))
+
+        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
+        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
+
+!!$        b^2 = (1-v^2)B^2+(B dot v)^2 
+        b2p=(1.d0-DOTP2(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,prim_p(2),prim_p(3),prim_p(4)))* &
+             DOTP2(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,mag_p(1),mag_p(2),mag_p(3)) + &
+             (DOTP(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,prim_p(2),prim_p(3),prim_p(4),mag_p(1),mag_p(2),mag_p(3)))**2
+        b2m=(1.d0-DOTP2(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,prim_m(2),prim_m(3),prim_m(4)))* &
+             DOTP2(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,mag_m(1),mag_m(2),mag_m(3)) + &
+             (DOTP(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,prim_m(2),prim_m(3),prim_m(4),mag_m(1),mag_m(2),mag_m(3)))**2
+        
+        vA2p = b2p/(prim_p(1)*(1.0d0+prim_p(5))+prim_p(6)+b2p)
+        vA2m = b2m/(prim_m(1)*(1.0d0+prim_m(5))+prim_m(6)+b2m)
+        
+!!$        Eigenvalues and fluxes either side of the cell interface
+        
+        if (flux_direction == 1) then
+           call eigenvalues_generalM(&
+                prim_m(2),prim_m(3),prim_m(4),cs2_m,vA2m, &
+                lamminus,&
+                gxxh,gxyh,gxzh,&
+                gyyh,gyzh,gzzh,&
+                usendh,avg_alp,avg_beta)
+           call eigenvalues_generalM(&
+                prim_p(2),prim_p(3),prim_p(4),cs2_p,vA2p, &
+                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_generalM(&
+                prim_m(3),prim_m(4),prim_m(2),cs2_m,vA2m, &
+                lamminus,&
+                gyyh,gyzh,gxyh,&
+                gzzh,gxzh,gxxh,&
+                usendh,avg_alp,avg_beta)
+           call eigenvalues_generalM(&
+                prim_p(3),prim_p(4),prim_p(2),cs2_p,vA2p, &
+                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_generalM(&
+                prim_m(4),prim_m(2),prim_m(3),cs2_m,vA2m, &
+                lamminus,&
+                gzzh,gxzh,gyzh,&
+                gxxh,gxyh,gyyh,&
+                usendh,avg_alp,avg_beta)
+           call eigenvalues_generalM(&
+                prim_p(4),prim_p(2),prim_p(3),cs2_p,vA2p, &
+                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
+           
+!!$        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_TracerGeneralM
+