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

diff --git a/src/GRHydro_Prim2Con.F90 b/src/GRHydro_Prim2Con.F90
new file mode 100644
index 0000000..9710a60
--- /dev/null
+++ b/src/GRHydro_Prim2Con.F90
@@ -0,0 +1,705 @@
+  /*@@
+   @file      primitive2conservative
+   @date      Thu Jan  11 11:03:32 2002
+   @author    Pedro Montero, Ian Hawke
+   @desc 
+   Primitive to conservative routine
+   @enddesc 
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+#include "cctk_Functions.h"
+
+#define velx(i,j,k) vel(i,j,k,1)
+#define vely(i,j,k) vel(i,j,k,2)
+#define velz(i,j,k) vel(i,j,k,3)
+#define sx(i,j,k) scon(i,j,k,1)
+#define sy(i,j,k) scon(i,j,k,2)
+#define sz(i,j,k) scon(i,j,k,3)
+
+
+ /*@@
+   @routine   primitive2conservative.f90 
+   @date       Thu Jan 11 11:03:32 2002
+   @author     Pedro Montero, Ian Hawke
+   @desc 
+   Converts primitive to conserved variables for the boundary extended data.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+subroutine primitive2conservative(CCTK_ARGUMENTS)
+
+  implicit none
+  
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  
+  integer :: i, j, k
+  CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,&
+       gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr,psi4l,psi4r
+  
+  do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1
+    do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1
+      do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1
+        
+        gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset))
+        gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset))
+        gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset))
+        gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset))
+        gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset))
+        gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset))
+        gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset))
+        gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset))
+        gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset))
+        gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset))
+        gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset))
+        gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset))
+
+        if (.not.(conformal_state .eq. 0)) then
+          psi4r = (0.5d0*(psi(i,j,k)+psi(i+xoffset,j+yoffset,k+zoffset)))**4
+          psi4l = (0.5d0*(psi(i,j,k)+psi(i-xoffset,j-yoffset,k-zoffset)))**4
+          gxxr = gxxr * psi4r
+          gxyr = gxyr * psi4r
+          gxzr = gxzr * psi4r
+          gyyr = gyyr * psi4r
+          gyzr = gyzr * psi4r
+          gzzr = gzzr * psi4r
+          gxxl = gxxl * psi4l
+          gxyl = gxyl * psi4l
+          gxzl = gxzl * psi4l
+          gyyl = gyyl * psi4l
+          gyzl = gyzl * psi4l
+          gzzl = gzzl * psi4l
+        end if
+
+        call SpatialDeterminant(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl,avg_detl)
+        call SpatialDeterminant(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr,avg_detr)
+
+        call prim2con(GRHydro_eos_handle, gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, &
+             avg_detl,densminus(i,j,k),sxminus(i,j,k),&
+             syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),rhominus(i,j,k), &
+             velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),&
+             epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i, j, k))
+ 
+        call prim2con(GRHydro_eos_handle, gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, &
+             avg_detr, densplus(i,j,k),sxplus(i,j,k),&
+             syplus(i,j,k),szplus(i,j ,k),tauplus(i,j,k),&
+             rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),&
+             velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),&
+             w_lorentzplus(i,j,k)) 
+
+      end do
+    end do
+  end do
+
+end subroutine primitive2conservative
+
+ /*@@
+   @routine    prim2con
+   @date       Sat Jan 26 01:52:18 2002
+   @author     Pedro Montero, Ian Hawke
+   @desc 
+   Converts from primitive to conservative at a single point
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+subroutine prim2con(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, &
+     dsx, dsy, dsz, dtau , drho, dvelx, dvely, dvelz, deps, dpress, w) 
+  
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det
+  CCTK_REAL :: ddens, dsx, dsy, dsz, dtau, drho, dvelx, dvely, dvelz,&
+       deps, dpress, w, vlowx, vlowy, vlowz   
+  CCTK_INT :: handle
+  character(len=256) NaN_WarnLine
+
+#include "EOS_Base.inc"
+  
+  w = 1.d0 / sqrt(1.d0 - (gxx*dvelx*dvelx + gyy*dvely*dvely + gzz &
+       *dvelz*dvelz + 2*gxy*dvelx*dvely + 2*gxz*dvelx *dvelz + 2*gyz&
+       *dvely*dvelz))  
+
+!!$ BEGIN: Check for NaN value
+  if (w .ne. w) then
+    !$OMP CRITICAL
+    write(NaN_WarnLine,'(a100,3g15.6)') 'NaN produced in sqrt(): (dvelx,dvely,dvelz)', dvelx, dvely, dvelz
+    call CCTK_WARN(GRHydro_NaN_verbose, NaN_WarnLine)
+    !$OMP END CRITICAL
+  endif
+!!$ END: Check for NaN value
+
+  dpress = EOS_Pressure(handle, drho, deps)
+
+  vlowx = gxx*dvelx + gxy*dvely + gxz*dvelz
+  vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz
+  vlowz = gxz*dvelx + gyz*dvely + gzz*dvelz
+
+  ddens = sqrt(det) * drho * w 
+  dsx = sqrt(det) * (drho*(1+deps)+dpress)*w*w * vlowx
+  dsy = sqrt(det) * (drho*(1+deps)+dpress)*w*w * vlowy
+  dsz = sqrt(det) * (drho*(1+deps)+dpress)*w*w * vlowz
+  dtau = sqrt(det) * ((drho*(1+deps)+dpress)*w*w - dpress) - ddens 
+
+end subroutine prim2con
+
+
+ /*@@
+   @routine    Primitive2ConservativeCells
+   @date       Sun Mar 10 21:16:20 2002
+   @author     
+   @desc 
+   Wrapper function that converts primitive to conservative at the 
+     cell centres. 
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+
+subroutine Primitive2ConservativeCells(CCTK_ARGUMENTS)
+
+  implicit none
+  
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  
+  CCTK_INT :: i, j, k
+  CCTK_REAL :: det,psi4pt
+  
+  do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1
+    do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1
+      do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1
+        
+        if (conformal_state .eq. 0) then
+          psi4pt = 1d0
+          call SpatialDeterminant(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),&
+               gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),det)
+        else
+          psi4pt = psi(i,j,k)**4
+          call SpatialDeterminant(psi4pt*gxx(i,j,k),psi4pt*gxy(i,j,k),&
+               psi4pt*gxz(i,j,k),psi4pt*gyy(i,j,k),psi4pt*gyz(i,j,k),&
+               psi4pt*gzz(i,j,k),det)
+        end if
+        call prim2con(GRHydro_eos_handle,psi4pt*gxx(i,j,k),&
+             psi4pt*gxy(i,j,k),psi4pt*gxz(i,j,k),&
+             psi4pt*gyy(i,j,k),psi4pt*gyz(i,j,k),psi4pt*gzz(i,j,k),&
+             det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
+             tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),&
+             eps(i,j,k),press(i,j,k),w_lorentz(i,j,k))
+
+      end do
+    end do
+  end do
+
+end subroutine Primitive2ConservativeCells
+
+
+ /*@@
+   @routine    Prim2ConservativePolytype
+   @date       Tue Mar 19 22:52:21 2002
+   @author     Ian Hawke
+   @desc 
+   Same as first routine, only for polytropes.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+
+subroutine Prim2ConservativePolytype(CCTK_ARGUMENTS)
+
+  implicit none
+  
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  
+  integer :: i, j, k
+  CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,&
+       gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr,psi4r,psi4l
+  
+  !$OMP PARALLEL DO PRIVATE(i, j, gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,&
+  !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr,psi4r,psi4l)
+  do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1
+    do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1
+      do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1
+        
+        gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset))
+        gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset))
+        gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset))
+        gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset))
+        gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset))
+        gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset))
+        gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset))
+        gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset))
+        gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset))
+        gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset))
+        gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset))
+        gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset))
+
+        if (.not.(conformal_state .eq. 0)) then
+          psi4r = (0.5d0*(psi(i,j,k)+psi(i+xoffset,j+yoffset,k+zoffset)))**4
+          psi4l = (0.5d0*(psi(i,j,k)+psi(i-xoffset,j-yoffset,k-zoffset)))**4
+          gxxr = gxxr * psi4r
+          gxyr = gxyr * psi4r
+          gxzr = gxzr * psi4r
+          gyyr = gyyr * psi4r
+          gyzr = gyzr * psi4r
+          gzzr = gzzr * psi4r
+          gxxl = gxxl * psi4l
+          gxyl = gxyl * psi4l
+          gxzl = gxzl * psi4l
+          gyyl = gyyl * psi4l
+          gyzl = gyzl * psi4l
+          gzzl = gzzl * psi4l
+        end if
+
+        call SpatialDeterminant(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl,avg_detl)
+        call SpatialDeterminant(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr,avg_detr)
+
+        call prim2conpolytype(GRHydro_eos_handle, gxxl,gxyl,gxzl,&
+             gyyl,gyzl,gzzl, &
+             avg_detl,densminus(i,j,k),sxminus(i,j,k),&
+             syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),rhominus(i,j,k), &
+             velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),&
+             epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i, j, k))
+ 
+        call prim2conpolytype(GRHydro_eos_handle, gxxr,gxyr,gxzr,&
+             gyyr,gyzr,gzzr, &
+             avg_detr, densplus(i,j,k),sxplus(i,j,k),&
+             syplus(i,j,k),szplus(i,j ,k),tauplus(i,j,k),&
+             rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),&
+             velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),&
+             w_lorentzplus(i,j,k)) 
+        if (densminus(i,j,k) .ne. densminus(i,j,k)) then
+        !$OMP CRITICAL
+        call CCTK_WARN(1, "NaN in densminus(i,j,k) (Prim2Con)")
+        !$OMP END CRITICAL
+        endif
+      end do
+    end do
+  end do
+  !$OMP END PARALLEL DO
+end subroutine Prim2ConservativePolytype
+
+ /*@@
+   @routine    prim2conpolytype
+   @date       Sat Jan 26 01:52:18 2002
+   @author     Pedro Montero, Ian Hawke
+   @desc 
+   Converts from primitive to conservative at a single point
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+subroutine prim2conpolytype(handle, gxx, gxy, gxz, gyy, gyz, &
+     gzz, det, ddens, &
+     dsx, dsy, dsz, dtau , drho, dvelx, dvely, dvelz, deps, dpress, w) 
+  
+  implicit none
+  
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det
+  CCTK_REAL :: ddens, dsx, dsy, dsz, dtau, drho, dvelx, dvely, dvelz,&
+       deps, dpress, w_tmp, w, vlowx, vlowy, vlowz, sqrtdet   
+  CCTK_INT :: handle
+  character(len=256) NaN_WarnLine
+
+#ifdef _EOS_BASE_INC_
+#undef _EOS_BASE_INC_
+#endif
+#include "EOS_Base.inc"
+  
+  w_tmp = gxx*dvelx*dvelx + gyy*dvely*dvely + gzz *dvelz*dvelz + &
+          2*gxy*dvelx*dvely + 2*gxz*dvelx*dvelz + 2*gyz*dvely*dvelz
+  if (w_tmp .ge. 1.d0) then
+    ! In theory this should not happen, and even when accepting the fact
+    ! that numerically it can, one might be tempted to set w to some large
+    ! value in that case. However, this would lead to completely bogus
+    ! and hard to trace wrong values below. There is no good value to
+    ! choose in this case, but something small is probably the best of
+    ! all bad choices.
+    !$OMP CRITICAL
+    write(NaN_WarnLine,'(a80,2g15.6)') 'Infinite Lorentz factor reset. rho, w_tmp: ', drho, w_tmp
+    call CCTK_WARN(GRHydro_NaN_verbose, NaN_WarnLine)
+    !$OMP END CRITICAL
+    w = 1.d-20
+  else
+    w = 1.d0 / sqrt(1.d0 - w_tmp)
+  endif
+  
+  if (handle .ge. 0) then
+    dpress = EOS_Pressure(handle, drho, deps)
+    deps = EOS_SpecificIntEnergy(handle, drho, dpress)
+  end if
+
+  vlowx = gxx*dvelx + gxy*dvely + gxz*dvelz
+  vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz
+  vlowz = gxz*dvelx + gyz*dvely + gzz*dvelz
+
+  sqrtdet = sqrt(det)
+  ddens = sqrtdet * drho * w 
+  dsx = sqrtdet * (drho*(1+deps)+dpress)*w*w * vlowx
+  dsy = sqrtdet * (drho*(1+deps)+dpress)*w*w * vlowy
+  dsz = sqrtdet * (drho*(1+deps)+dpress)*w*w * vlowz
+  dtau = sqrtdet * ((drho*(1+deps)+dpress)*w*w - dpress) - ddens 
+
+end subroutine prim2conpolytype
+
+
+ /*@@
+   @routine    Primitive2ConservativePolyCells
+   @date       Sun Mar 10 21:16:20 2002
+   @author     
+   @desc 
+   Wrapper function that converts primitive to conservative at the 
+     cell centres. 
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+
+subroutine Primitive2ConservativePolyCells(CCTK_ARGUMENTS)
+
+  implicit none
+  
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  
+  CCTK_INT :: i, j, k
+  CCTK_REAL :: det
+  
+  do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1
+    do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1
+      do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1
+        
+        call SpatialDeterminant(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),&
+             gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),det)
+        call prim2conpolytype(GRHydro_eos_handle,gxx(i,j,k),gxy(i,j,k),&
+             gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),&
+             det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
+             tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),&
+             eps(i,j,k),press(i,j,k),w_lorentz(i,j,k))
+
+      end do
+    end do
+  end do
+
+end subroutine Primitive2ConservativePolyCells
+
+ /*@@
+   @routine    Prim2ConservativeTracer
+   @date       Mon Mar  8 13:32:32 2004
+   @author     Ian Hawke
+   @desc 
+   Gets the conserved tracer variable from the primitive.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+subroutine Prim2ConservativeTracer(CCTK_ARGUMENTS)
+
+  implicit none
+  
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  
+  integer :: i, j, k
+  CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,&
+       gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr,psi4r,psi4l
+  
+  do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1
+    do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1
+      do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1
+        
+        gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset))
+        gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset))
+        gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset))
+        gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset))
+        gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset))
+        gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset))
+        gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset))
+        gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset))
+        gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset))
+        gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset))
+        gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset))
+        gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset))
+
+        if (.not.(conformal_state .eq. 0)) then
+          psi4r = (0.5d0*(psi(i,j,k)+psi(i+xoffset,j+yoffset,k+zoffset)))**4
+          psi4l = (0.5d0*(psi(i,j,k)+psi(i-xoffset,j-yoffset,k-zoffset)))**4
+          gxxr = gxxr * psi4r
+          gxyr = gxyr * psi4r
+          gxzr = gxzr * psi4r
+          gyyr = gyyr * psi4r
+          gyzr = gyzr * psi4r
+          gzzr = gzzr * psi4r
+          gxxl = gxxl * psi4l
+          gxyl = gxyl * psi4l
+          gxzl = gxzl * psi4l
+          gyyl = gyyl * psi4l
+          gyzl = gyzl * psi4l
+          gzzl = gzzl * psi4l
+        end if
+
+        call SpatialDeterminant(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl,avg_detl)
+        call SpatialDeterminant(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr,avg_detr)
+
+        cons_tracerplus(i,j,k,:) = tracerplus(i,j,k,:) * &
+             sqrt(avg_detr) * rhoplus(i,j,k) / &
+             sqrt(1.d0 - &
+                   (gxxr * velxplus(i,j,k)**2 + &
+                    gyyr * velyplus(i,j,k)**2 + &
+                    gzzr * velzplus(i,j,k)**2 + &
+                    2.d0 * (gxyr * velxplus(i,j,k) * velyplus(i,j,k) + &
+                            gxzr * velxplus(i,j,k) * velzplus(i,j,k) + &
+                            gyzr * velyplus(i,j,k) * velzplus(i,j,k) ) ) )
+        cons_tracerminus(i,j,k,:) = tracerminus(i,j,k,:) * &
+             sqrt(avg_detl) * rhominus(i,j,k) / &
+             sqrt(1.d0 - &
+                   (gxxl * velxminus(i,j,k)**2 + &
+                    gyyl * velyminus(i,j,k)**2 + &
+                    gzzl * velzminus(i,j,k)**2 + &
+                    2.d0 * (gxyl * velxminus(i,j,k) * velyminus(i,j,k) + &
+                            gxzl * velxminus(i,j,k) * velzminus(i,j,k) + &
+                            gyzl * velyminus(i,j,k) * velzminus(i,j,k) ) ) )
+
+      end do
+    end do
+  end do
+
+end subroutine Prim2ConservativeTracer
+
+ /*@@
+   @routine   primitive2conservativegeneral
+   @date       Thu Jan 11 11:03:32 2002
+   @author     Pedro Montero, Ian Hawke
+   @desc 
+   Converts primitive to conserved variables for the boundary extended data.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+subroutine primitive2conservativegeneral(CCTK_ARGUMENTS)
+
+  USE GRHydro_Scalars
+
+  implicit none
+  
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_FUNCTIONS
+  
+  CCTK_INT :: i, j, k, ierr
+  CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl, &
+       gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr,psi4l,psi4r, &
+       vlowx, vlowy, vlowz
+  
+  ierr = EOS_SetGFs(cctkGH, EOS_Prim2ConBndCallPlus)
+  ierr = EOS_SetGFs(cctkGH, EOS_Prim2ConBndCallMinus)
+
+  !$OMP PARALLEL DO PRIVATE(i, j, gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,&
+  !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr,psi4r,psi4l)
+  do k = GRHydro_stencil , cctk_lsh(3) - GRHydro_stencil + 1
+    do j = GRHydro_stencil , cctk_lsh(2) - GRHydro_stencil + 1
+      do i = GRHydro_stencil , cctk_lsh(1) - GRHydro_stencil + 1
+        
+        gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset))
+        gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset))
+        gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset))
+        gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset))
+        gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset))
+        gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset))
+        gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset))
+        gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset))
+        gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset))
+        gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset))
+        gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset))
+        gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset))
+
+        if (.not.(conformal_state .eq. 0)) then
+          psi4r = (0.5d0*(psi(i,j,k)+psi(i+xoffset,j+yoffset,k+zoffset)))**4
+          psi4l = (0.5d0*(psi(i,j,k)+psi(i-xoffset,j-yoffset,k-zoffset)))**4
+          gxxr = gxxr * psi4r
+          gxyr = gxyr * psi4r
+          gxzr = gxzr * psi4r
+          gyyr = gyyr * psi4r
+          gyzr = gyzr * psi4r
+          gzzr = gzzr * psi4r
+          gxxl = gxxl * psi4l
+          gxyl = gxyl * psi4l
+          gxzl = gxzl * psi4l
+          gyyl = gyyl * psi4l
+          gyzl = gyzl * psi4l
+          gzzl = gzzl * psi4l
+        end if
+
+        call SpatialDeterminant(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl,avg_detl)
+        call SpatialDeterminant(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr,avg_detr)
+
+        call prim2conpolytype(-1, gxxl,gxyl,gxzl,&
+             gyyl,gyzl,gzzl, &
+             avg_detl,densminus(i,j,k),sxminus(i,j,k),&
+             syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),rhominus(i,j,k), &
+             velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),&
+             epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i, j, k))
+
+        call prim2conpolytype(-1, gxxr,gxyr,gxzr,&
+             gyyr,gyzr,gzzr, &
+             avg_detr, densplus(i,j,k),sxplus(i,j,k),&
+             syplus(i,j,k),szplus(i,j ,k),tauplus(i,j,k),&
+             rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),&
+             velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),&
+             w_lorentzplus(i,j,k))
+        if (densminus(i,j,k) .ne. densminus(i,j,k)) then
+        !$OMP CRITICAL
+          call CCTK_WARN(1, "NaN in densminus(i,j,k) (Prim2Con)")
+        !$OMP END CRITICAL
+        endif
+      end do
+    end do
+  end do
+  !$OMP END PARALLEL DO
+
+end subroutine primitive2conservativegeneral
+
+ /*@@
+   @routine   Primitive2ConservativeCellsGeneral
+   @date       Thu Jan 11 11:03:32 2002
+   @author     Pedro Montero, Ian Hawke
+   @desc 
+   Converts primitive to conserved variables for the boundary extended data.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+subroutine Primitive2ConservativeCellsGeneral(CCTK_ARGUMENTS)
+
+  USE GRHydro_Scalars
+
+  implicit none
+  
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_FUNCTIONS
+  
+  CCTK_INT :: i, j, k, ierr
+  CCTK_REAL :: det,psi4pt,gxxpt,gxypt,gxzpt,gyypt,gyzpt,gzzpt, &
+       vlowx, vlowy, vlowz
+  
+  ierr = EOS_SetGFs(cctkGH, EOS_Prim2ConCellsCall)
+
+  do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + 1
+    do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + 1
+      do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + 1
+        
+        gxxpt = gxx(i,j,k)
+        gxypt = gxy(i,j,k)
+        gxzpt = gxz(i,j,k)
+        gyypt = gyy(i,j,k)
+        gyzpt = gyz(i,j,k)
+        gzzpt = gzz(i,j,k)
+
+        if (.not.(conformal_state .eq. 0)) then
+          psi4pt = psi(i,j,k)**4
+          gxxpt = gxxpt * psi4pt
+          gxypt = gxypt * psi4pt
+          gxzpt = gxzpt * psi4pt
+          gyypt = gyypt * psi4pt
+          gyzpt = gyzpt * psi4pt
+          gzzpt = gzzpt * psi4pt
+        end if
+
+        call SpatialDeterminant(gxxpt,gxypt,gxzpt,gyypt,gyzpt,gzzpt,det)
+
+        w_lorentz(i,j,k) = 1.d0 / sqrt(1.d0 - &
+             (gxxpt * velx(i,j,k)**2 + &
+              gyypt * vely(i,j,k)**2 + &
+              gzzpt * velz(i,j,k)**2 + &
+              2.d0 * gxypt * velx(i,j,k) * vely(i,j,k) + &
+              2.d0 * gxzpt * velx(i,j,k) * velz(i,j,k) + &
+              2.d0 * gyzpt * vely(i,j,k) * velz(i,j,k) ))
+        vlowx = gxxpt * velx(i,j,k) + &
+                gxypt * vely(i,j,k) + &
+                gxzpt * velz(i,j,k)
+        vlowy = gxypt * velx(i,j,k) + &
+                gyypt * vely(i,j,k) + &
+                gyzpt * velz(i,j,k)
+        vlowz = gxzpt * velx(i,j,k) + &
+                gyzpt * vely(i,j,k) + &
+                gzzpt * velz(i,j,k)
+        dens(i,j,k) = sqrt(det) * rho(i,j,k) * &
+             w_lorentz(i,j,k)
+        sx(i,j,k) =  sqrt(det) * &
+             ( rho(i,j,k) * &
+               (1.d0 + eps(i,j,k)) + press(i,j,k) ) * &
+               w_lorentz(i,j,k)**2 * vlowx
+        sy(i,j,k) =  sqrt(det) * &
+             ( rho(i,j,k) * &
+               (1.d0 + eps(i,j,k)) + press(i,j,k) ) * &
+               w_lorentz(i,j,k)**2 * vlowy
+        sz(i,j,k) =  sqrt(det) * &
+             ( rho(i,j,k) * &
+               (1.d0 + eps(i,j,k)) + press(i,j,k) ) * &
+               w_lorentz(i,j,k)**2 * vlowz
+        tau(i,j,k) = sqrt(det) * &
+             ( ( rho(i,j,k) * &
+                 (1.d0 + eps(i,j,k)) + press(i,j,k) ) * &
+                 w_lorentz(i,j,k)**2 - press(i,j,k) ) - &
+             dens(i,j,k)
+
+      end do
+    end do
+  end do
+
+end subroutine Primitive2ConservativeCellsGeneral