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, 587 insertions, 0 deletions

diff --git a/src/GRHydro_RoeSolver.F90 b/src/GRHydro_RoeSolver.F90
new file mode 100644
index 0000000..c661ed1
--- /dev/null
+++ b/src/GRHydro_RoeSolver.F90
@@ -0,0 +1,587 @@
+ /*@@
+   @file      GRHydro_RoeSolver.F90
+   @date      Sat Jan 26 01:55:27 2002
+   @author    
+   @desc 
+   Calculates the Roe flux
+   @enddesc 
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+#include "cctk_Functions.h"
+
+#include "SpaceMask.h"
+
+ /*@@
+   @routine    GRHydro_RoeSolve
+   @date       Sat Jan 26 01:55:55 2002
+   @author     Pedro Montero, Ian Hawke
+   @desc 
+   Wrapper routine to calculate the Roe fluxes and hence the update 
+   terms.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+   
+   @endhistory 
+
+@@*/
+
+subroutine GRHydro_RoeSolve(CCTK_ARGUMENTS)
+  
+  USE GRHydro_Eigenproblem
+  implicit none
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  CCTK_REAL, dimension(5) :: roeflux,roeave,qdiff,consp,consm_i,&
+       fplus,fminus,f_roe,primp,primm_i,consh
+  CCTK_REAL :: avg_alp,avg_beta,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, &
+       avg_det,uxxh,uxyh,uxzh,uyyh,uyzh,uzzh,&
+       rhoave, velxave, velyave, velzave, pressave, epsave, &
+       w_lorentzave, usendh, alp_l, alp_r, psi4h
+  integer :: m
+  integer :: i,j,k
+  
+  CCTK_INT :: type_bits, trivial, not_trivial
+
+  character(len=256) NaN_WarnLine
+
+  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
+  
+  f_roe = 0.d0
+  
+  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
+
+!!$        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) 
+        
+        primp(1)   = rhoplus(i,j,k) 
+        primp(2)   = velxplus(i,j,k)
+        primp(3)   = velyplus(i,j,k)
+        primp(4)   = velzplus(i,j,k)
+        primp(5)   = epsplus(i,j,k)
+        
+        primm_i(1) = rhominus(i+xoffset,j+yoffset,k+zoffset)
+        primm_i(2) = velxminus(i+xoffset,j+yoffset,k+zoffset)
+        primm_i(3) = velyminus(i+xoffset,j+yoffset,k+zoffset)
+        primm_i(4) = velzminus(i+xoffset,j+yoffset,k+zoffset)
+        primm_i(5) = epsminus(i+xoffset,j+yoffset,k+zoffset) 
+        
+        roeflux = 0.d0
+        qdiff = 0.d0
+
+!!$        Calculate jumps in conserved variables
+        
+        do m = 1,5
+          qdiff(m) = consm_i(m) - consp(m)
+        end do
+
+!!$        Set metric terms at interface
+        
+        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 > 0) then
+          
+          psi4h = (0.5d0*(psi(i,j,k)+psi(i+xoffset,j+yoffset,k+zoffset)))**4
+          gxxh = gxxh * psi4h
+          gxyh = gxyh * psi4h
+          gxzh = gxzh * psi4h
+          gyyh = gyyh * psi4h
+          gyzh = gyzh * psi4h
+          gzzh = gzzh * psi4h
+          
+        end if
+
+        call SpatialDeterminant(gxxh,gxyh,gxzh,gyyh,gyzh,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(consp(1),consp(2),consp(3),consp(4),consp(5),&
+                 f_roe(1),f_roe(2),f_roe(3),&
+                 f_roe(4),f_roe(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),&
+                 f_roe(1),f_roe(3),f_roe(4),&
+                 f_roe(2),f_roe(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),&
+                 f_roe(1),f_roe(4),f_roe(2),&
+                 f_roe(3),f_roe(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
+          
+!!$        Set the Roe average of the fluid variables
+        
+          call roeaverage(primp, primm_i, roeave)
+        
+          rhoave = roeave(1)
+          velxave = roeave(2)
+          velyave = roeave(3)
+          velzave = roeave(4)
+          epsave = roeave(5)
+
+!!$        Convert to conserved variables and find the part of the Roe
+!!$        flux that requires the spectral decomposition.
+!!$        The conversion to conserved variables is just to set the 
+!!$        pressure at this point (means this routine doesn''t need
+!!$        the EOS interface).
+
+!!$        The conversion routine is unnecessary (the pressure is set
+!!$        inside the eigenproblem routine) so instead we just have
+!!$        to set the average W.
+
+          w_lorentzave = 1.d0 / &
+               sqrt(1.d0 - &
+                 (gxxh*velxave*velxave + gyyh*velyave*velyave + &
+                  gzzh*velzave*velzave + 2*gxyh*velxave*velyave + &
+                  2*gxzh*velxave *velzave + 2*gyzh*velyave*velzave))  
+
+!!$ BEGIN: Check for NaN value
+  if (w_lorentzave .ne. w_lorentzave) then
+    write(NaN_WarnLine,'(a100,3g15.6)') 'NaN produced in sqrt(): (velxave, velyave, velzave))', velxave, velyave, velzave
+    call CCTK_WARN(GRHydro_NaN_verbose, NaN_WarnLine)
+  endif
+!!$ END: Check for NaN value
+
+          if (flux_direction == 1) then
+!!$            call prim2con(GRHydro_eos_handle,gxxh, gxyh, gxzh, gyyh, &
+!!$                 gyzh, gzzh, avg_det, &
+!!$                 consh(1), consh(2), consh(3), consh(4), consh(5), rhoave, &
+!!$                 velxave, velyave, velzave, epsave, pressave, w_lorentzave)
+            call eigenproblem(GRHydro_eos_handle,rhoave, velxave, &
+                 velyave, velzave, epsave, w_lorentzave, &
+                 gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, &
+                 usendh, avg_alp,avg_beta,qdiff(1),qdiff(2), &
+                 qdiff(3),qdiff(4),qdiff(5),roeflux(1),roeflux(2),&
+                 roeflux(3),roeflux(4),roeflux(5))
+          else if (flux_direction == 2) then
+!!$            call prim2con(GRHydro_eos_handle,gyyh, gyzh, gxyh, gzzh, &
+!!$                 gxzh, gxxh, avg_det, &
+!!$                 consh(1), consh(3), consh(4), consh(2), consh(5), rhoave, &
+!!$                 velyave, velzave, velxave, epsave, pressave, w_lorentzave)
+            call eigenproblem(GRHydro_eos_handle,rhoave, velyave, &
+                 velzave, velxave, epsave, w_lorentzave, &
+                 gyyh,gyzh,gxyh,gzzh,gxzh,gxxh, &
+                 usendh, avg_alp,avg_beta,qdiff(1),qdiff(3), &
+                 qdiff(4),qdiff(2),qdiff(5),roeflux(1),roeflux(3),&
+                 roeflux(4),roeflux(2),roeflux(5))
+          else if (flux_direction == 3) then
+!!$            call prim2con(GRHydro_eos_handle,gzzh, gxzh, gyzh, gxxh, &
+!!$                 gxyh, gyyh, avg_det, &
+!!$                 consh(1), consh(4), consh(2), consh(3), consh(5), rhoave, &
+!!$                 velzave, velxave, velyave, epsave, pressave, w_lorentzave)
+            call eigenproblem(GRHydro_eos_handle,rhoave, velzave, &
+                 velxave, velyave, epsave, w_lorentzave, &
+                 gzzh,gxzh,gyzh,gxxh,gxyh,gyyh, &
+                 usendh, avg_alp,avg_beta,qdiff(1),qdiff(4), &
+                 qdiff(2),qdiff(3),qdiff(5),roeflux(1),roeflux(4),&
+                 roeflux(2),roeflux(3),roeflux(5))
+          else
+            call CCTK_WARN(0, "Flux direction not x,y,z")
+          end if
+          
+          fplus = 0.d0
+          fminus = 0.d0
+        
+!!$Calculate the fluxes of the original reconstructed data
+        
+          if (flux_direction == 1) then
+            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
+            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
+            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
+        
+!!$The combined Roe flux          
+        
+          do m = 1,5
+          
+            f_roe(m) = 0.5d0 * (fplus(m) + fminus(m) - roeflux(m))
+          
+          end do
+
+        end if !!! The end of the SpaceMask check for a trivial RP.
+        
+        densflux(i,j,k) = f_roe(1)
+        sxflux(i,j,k)   = f_roe(2)
+        syflux(i,j,k)   = f_roe(3)
+        szflux(i,j,k)   = f_roe(4)
+        tauflux(i,j,k)  = f_roe(5)
+
+      enddo
+    enddo
+  enddo
+    
+end subroutine GRHydro_RoeSolve
+
+ /*@@
+   @routine    GRHydro_RoeSolveGeneral
+   @date       Sat Jan 26 01:55:55 2002
+   @author     Pedro Montero, Ian Hawke
+   @desc 
+   Wrapper routine to calculate the Roe fluxes and hence the update 
+   terms.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+   
+   @endhistory 
+
+@@*/
+
+subroutine GRHydro_RoeSolveGeneral(CCTK_ARGUMENTS)
+  
+  USE GRHydro_Scalars
+  use GRHydro_Eigenproblem
+
+  implicit none
+
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_FUNCTIONS
+
+  CCTK_REAL, dimension(5) :: roeflux,roeave,qdiff,cons_p,cons_m,&
+       cons_ave
+  CCTK_REAL, dimension(6) :: prim_p, prim_m, prim_ave
+  CCTK_REAL :: avg_alp,avg_beta,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, &
+       avg_det,uxxh,uxyh,uxzh,uyyh,uyzh,uzzh,&
+       usendh,alp_l,alp_r,psi4h,cs2_ave,dpdeps_ave
+  CCTK_INT :: i,j,k,m
+  
+  CCTK_INT :: type_bits, trivial, not_trivial, ierr
+
+  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
+
+!!$        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)
+
+        prim_m(6) = pressminus(i+xoffset,j+yoffset,k+zoffset)
+
+!!$        Set the Roe average of the fluid variables
+        
+        call roeaverage(prim_p(1:5), prim_m(1:5), roeave)
+        
+        rho_ave(i,j,k)  = roeave(1)
+        velx_ave(i,j,k) = roeave(2)
+        vely_ave(i,j,k) = roeave(3)
+        velz_ave(i,j,k) = roeave(4)
+        eps_ave(i,j,k)  = roeave(5)
+
+      end do
+    end do
+  end do
+  
+  ierr = EOS_SetGFs(cctkGH, EOS_RoeAverageCall)
+  
+  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
+
+!!$        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_ave(1)   = rho_ave(i,j,k) 
+        prim_ave(2)   = velx_ave(i,j,k)
+        prim_ave(3)   = vely_ave(i,j,k)
+        prim_ave(4)   = velz_ave(i,j,k)
+        prim_ave(5)   = eps_ave(i,j,k)
+        prim_ave(6)   = press_ave(i,j,k)
+        cs2_ave       = eos_cs2_ave(i,j,k)
+        dpdeps_ave    = eos_dpdeps_ave(i,j,k)
+
+        roeflux = 0.d0
+        qdiff = 0.d0
+
+!!$        Calculate jumps in conserved variables
+        
+        do m = 1,5
+          qdiff(m) = cons_m(m) - cons_p(m)
+        end do
+
+!!$        Set metric terms at interface
+        
+        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 > 0) then
+          
+          psi4h = (0.5d0*(psi(i,j,k)+psi(i+xoffset,j+yoffset,k+zoffset)))**4
+          gxxh = gxxh * psi4h
+          gxyh = gxyh * psi4h
+          gxzh = gxzh * psi4h
+          gyyh = gyyh * psi4h
+          gyzh = gyzh * psi4h
+          gzzh = gzzh * psi4h
+          
+        end if
+
+        call SpatialDeterminant(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,avg_det)
+        
+!!$ If the Riemann problem is trivial the flux is already correct
+          
+        if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, trivial)) then
+            
+        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
+          
+!!$        Convert to conserved variables and find the part of the Roe
+!!$        flux that requires the spectral decomposition.
+!!$        The conversion to conserved variables is just to set the 
+!!$        pressure at this point (means this routine doesn''t need
+!!$        the EOS interface).
+
+          if (flux_direction == 1) then
+            call eigenproblem_general(&
+                 prim_ave(1),prim_ave(2),prim_ave(3),prim_ave(4),prim_ave(5), &
+                 prim_ave(6),cs2_ave,dpdeps_ave, &
+                 gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, &
+                 usendh,avg_alp,avg_beta, &
+                 qdiff(1),qdiff(2),qdiff(3),qdiff(4),qdiff(5), &
+                 roeflux(1),roeflux(2),roeflux(3),roeflux(4),roeflux(5))
+          else if (flux_direction == 2) then
+            call eigenproblem_general(&
+                 prim_ave(1),prim_ave(3),prim_ave(4),prim_ave(2),prim_ave(5), &
+                 prim_ave(6),cs2_ave,dpdeps_ave, &
+                 gyyh,gyzh,gxyh,gzzh,gxzh,gxxh, &
+                 usendh,avg_alp,avg_beta, &
+                 qdiff(1),qdiff(3),qdiff(4),qdiff(2),qdiff(5), &
+                 roeflux(1),roeflux(3),roeflux(4),roeflux(2),roeflux(5))
+          else if (flux_direction == 3) then
+            call eigenproblem_general(&
+                 prim_ave(1),prim_ave(4),prim_ave(2),prim_ave(3),prim_ave(5), &
+                 prim_ave(6),cs2_ave,dpdeps_ave, &
+                 gzzh,gxzh,gyzh,gxxh,gxyh,gyyh, &
+                 usendh,avg_alp,avg_beta, &
+                 qdiff(1),qdiff(4),qdiff(2),qdiff(3),qdiff(5), &
+                 roeflux(1),roeflux(4),roeflux(2),roeflux(3),roeflux(5))
+          else
+            call CCTK_WARN(0, "Flux direction not x,y,z")
+          end if
+                  
+!!$ Roe flux
+          
+          densflux(i,j,k) = densflux(i,j,k) - 0.5d0 * roeflux(1)
+          sxflux(i,j,k)   = sxflux(i,j,k)   - 0.5d0 * roeflux(2)
+          syflux(i,j,k)   = syflux(i,j,k)   - 0.5d0 * roeflux(3)
+          szflux(i,j,k)   = szflux(i,j,k)   - 0.5d0 * roeflux(4)
+          tauflux(i,j,k)  = tauflux(i,j,k)  - 0.5d0 * roeflux(5)
+
+        end if !!! The end of the SpaceMask check for a trivial RP.
+
+      enddo
+    enddo
+  enddo
+    
+end subroutine GRHydro_RoeSolveGeneral