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

diff --git a/src/GRHydro_ReconstructPoly.F90 b/src/GRHydro_ReconstructPoly.F90
new file mode 100644
index 0000000..0467c27
--- /dev/null
+++ b/src/GRHydro_ReconstructPoly.F90
@@ -0,0 +1,530 @@
+ /*@@
+   @file      GRHydro_ReconstructPoly.F90
+   @date      Sat Jan 26 02:13:25 2002
+   @author    
+   @desc 
+   Wrapper routine to perform the reconstruction for polytropes.
+   @enddesc 
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+#include "cctk_Functions.h"
+
+#include "SpaceMask.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    ReconstructionPolytype
+   @date       Tue Mar 19 11:36:55 2002
+   @author     Ian Hawke  
+   @desc 
+   If using a polytropic type EOS, we do not have to reconstruct all the 
+   variables.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+subroutine ReconstructionPolytype(CCTK_ARGUMENTS)
+  
+  implicit none
+  
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_FUNCTIONS
+
+  integer :: nx, ny, nz, i, j, k, itracer
+
+  logical, dimension(:,:,:), allocatable :: trivial_rp
+!!$  logical, dimension(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) :: trivial_rp
+
+  CCTK_INT :: type_bitsx, trivialx, not_trivialx, &
+       &type_bitsy, trivialy, not_trivialy, &
+       &type_bitsz, trivialz, not_trivialz
+
+  CCTK_REAL, dimension(:,:,:),allocatable :: &
+       &psi4, lbetax, lbetay, lbetaz
+!!$  CCTK_REAL, dimension(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) :: &
+!!$       &psi4, lbetax, lbetay, lbetaz
+
+  CCTK_INT :: ierr
+
+  CCTK_REAL :: local_min_tracer
+
+  allocate(trivial_rp(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),STAT=ierr)
+  if (ierr .ne. 0) then
+    call CCTK_WARN(0, "Allocation problems with trivial_rp")
+  end if
+  
+  allocate(psi4(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),&
+       lbetax(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),&
+       lbetay(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),&
+       lbetaz(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),STAT=ierr)
+  if (ierr .ne. 0) then
+    call CCTK_WARN(0, "Allocation problems with lbeta")
+  end if
+  
+  call SpaceMask_GetTypeBits(type_bitsx, "Hydro_RiemannProblemX")
+  call SpaceMask_GetStateBits(trivialx, "Hydro_RiemannProblemX", &
+       &"trivial")
+  call SpaceMask_GetStateBits(not_trivialx, "Hydro_RiemannProblemX", &
+       &"not_trivial")
+  call SpaceMask_GetTypeBits(type_bitsy, "Hydro_RiemannProblemY")
+  call SpaceMask_GetStateBits(trivialy, "Hydro_RiemannProblemY", &
+       &"trivial")
+  call SpaceMask_GetStateBits(not_trivialy, "Hydro_RiemannProblemY", &
+       &"not_trivial")
+  call SpaceMask_GetTypeBits(type_bitsz, "Hydro_RiemannProblemZ")
+  call SpaceMask_GetStateBits(trivialz, "Hydro_RiemannProblemZ", &
+       &"trivial")
+  call SpaceMask_GetStateBits(not_trivialz, "Hydro_RiemannProblemZ", &
+       &"not_trivial")
+
+  nx = cctk_lsh(1)
+  ny = cctk_lsh(2)
+  nz = cctk_lsh(3)
+
+  trivial_rp = .false.
+
+!!$  Currently only option is reconstruction on primitive variables.
+!!$  Should change this.
+
+  if (conformal_state .eq. 0) then
+    psi4 = 1.d0
+  else
+    psi4 = psi**4
+  end if
+  if (shift_state .ne. 0) then
+    lbetax = betax
+    lbetay = betay
+    lbetaz = betaz
+  else
+    lbetax = 0.d0
+    lbetay = 0.d0
+    lbetaz = 0.d0
+  end if
+
+!!$ Initialize variables that store reconstructed quantities
+
+  rhoplus = 0.0d0
+  rhominus = 0.0d0
+  epsplus = 0.0d0
+  epsminus = 0.0d0
+  velxplus = 0.0d0
+  velxminus = 0.0d0
+  velyplus = 0.0d0
+  velyminus = 0.0d0
+  velzplus = 0.0d0
+  velzminus = 0.0d0
+
+  if (evolve_tracer .ne. 0) then
+     tracerplus = 0.0d0
+     tracerminus = 0.0d0
+  endif
+
+  if (CCTK_EQUALS(recon_method,"tvd")) then
+
+     if (evolve_tracer .ne. 0) then
+        do itracer=1,number_of_tracers
+           call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
+                tracer(:,:,:,itracer), tracerplus(:,:,:,itracer), &
+                tracerminus(:,:,:,itracer), &
+                trivial_rp, hydro_excision_mask)
+        enddo
+     end if
+     
+
+    if (CCTK_EQUALS(recon_vars,"primitive")) then
+      call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
+           rho, rhoplus, rhominus, trivial_rp, hydro_excision_mask)
+      call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
+           vel(:,:,:,1), velxplus, velxminus, trivial_rp, hydro_excision_mask)
+      call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
+           vel(:,:,:,2), velyplus, velyminus, trivial_rp, hydro_excision_mask)
+      call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
+           vel(:,:,:,3), velzplus, velzminus, trivial_rp, hydro_excision_mask)
+
+    else if (CCTK_EQUALS(recon_vars,"conservative")) then
+      call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
+           dens, densplus, densminus, trivial_rp, hydro_excision_mask)
+      call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
+           scon(:,:,:,1), sxplus, sxminus, trivial_rp, hydro_excision_mask)
+      call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
+           scon(:,:,:,2), syplus, syminus, trivial_rp, hydro_excision_mask)
+      call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
+           scon(:,:,:,3), szplus, szminus, trivial_rp, hydro_excision_mask)
+
+    else
+      call CCTK_WARN(0, "Variable type to reconstruct not recognized.")
+    end if
+
+    !$OMP PARALLEL DO PRIVATE(i, j)
+    do k = 1, nz
+      do j = 1, ny
+        do i = 1, nx
+          if (trivial_rp(i,j,k)) then
+            if (flux_direction == 1) then
+              SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, trivialx)
+            else if (flux_direction == 2) then
+              SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsy, trivialy)
+            else if (flux_direction == 3) then
+              SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsz, trivialz)
+            end if
+          else
+            if (flux_direction == 1) then
+              SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, not_trivialx)
+            else if (flux_direction == 2) then
+              SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsy, not_trivialy)
+            else if (flux_direction == 3) then
+              SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsz, not_trivialz)
+            end if
+          end if
+        end do
+      end do
+    end do
+    !$OMP END PARALLEL DO
+
+  else if (CCTK_EQUALS(recon_method,"ppm")) then
+
+    if (flux_direction == 1) then
+      !$OMP PARALLEL DO PRIVATE(i, j)
+      do k = GRHydro_stencil, nz - GRHydro_stencil + 1
+        do j = GRHydro_stencil, ny - GRHydro_stencil + 1
+          call SimplePPM_1d(GRHydro_eos_handle,1,nx,CCTK_DELTA_SPACE(1),&
+               rho(:,j,k),velx(:,j,k),vely(:,j,k),velz(:,j,k),eps(:,j,k),&
+               press(:,j,k),rhominus(:,j,k),velxminus(:,j,k),velyminus(:,j,k),&
+               velzminus(:,j,k),epsminus(:,j,k),rhoplus(:,j,k),&
+               velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),epsplus(:,j,k),&
+               trivial_rp(:,j,k), hydro_excision_mask(:,j,k),&
+               gxx(:,j,k), gxy(:,j,k), gxz(:,j,k), gyy(:,j,k), gyz(:,j,k), &
+               gzz(:,j,k), psi4(:,j,k), lbetax(:,j,k), alp(:,j,k),&
+               w_lorentz(:,j,k), &
+               1, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_x_left, &
+                                    GRHydro_mppm_eigenvalue_x_right, &
+                                    GRHydro_mppm_xwind)
+          do i = 1, nx
+            if (trivial_rp(i,j,k)) then
+              SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, trivialx)
+            else
+              SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, not_trivialx)
+            end if
+          end do
+        end do
+      end do
+      !$OMP END PARALLEL DO
+      if(evolve_tracer.ne.0) then
+         !$OMP PARALLEL DO PRIVATE(j)
+         do k = GRHydro_stencil, nz - GRHydro_stencil + 1
+            do j = GRHydro_stencil, ny - GRHydro_stencil + 1
+               call SimplePPM_tracer_1d(nx,CCTK_DELTA_SPACE(1),rho(:,j,k), &
+                    velx(:,j,k),vely(:,j,k),velz(:,j,k), &
+                    tracer(:,j,k,:),tracerminus(:,j,k,:),tracerplus(:,j,k,:), &
+                    press(:,j,k))
+            end do
+         end do
+         !$OMP END PARALLEL DO
+      end if
+
+    else if (flux_direction == 2) then
+      !$OMP PARALLEL DO PRIVATE(i, j)
+      do k = GRHydro_stencil, nz - GRHydro_stencil + 1
+        do j = GRHydro_stencil, nx - GRHydro_stencil + 1
+          call SimplePPM_1d(GRHydro_eos_handle,1,ny,CCTK_DELTA_SPACE(2),&
+               rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),eps(j,:,k),&
+               press(j,:,k),rhominus(j,:,k),velyminus(j,:,k),velzminus(j,:,k),&
+               velxminus(j,:,k),epsminus(j,:,k),rhoplus(j,:,k),&
+               velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),epsplus(j,:,k),&
+               trivial_rp(j,:,k), hydro_excision_mask(j,:,k),&
+               gyy(j,:,k), gyz(j,:,k), gxy(j,:,k), gzz(j,:,k), gxz(j,:,k), &
+               gxx(j,:,k), psi4(j,:,k), lbetay(j,:,k), alp(j,:,k),&
+               w_lorentz(j,:,k), &
+               2, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_y_left, &
+                                    GRHydro_mppm_eigenvalue_y_right, &
+                                    GRHydro_mppm_xwind)
+          do i = 1, ny
+            if (trivial_rp(j,i,k)) then
+              SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, trivialy)
+            else
+              SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, not_trivialy)
+            end if
+          end do
+        end do
+      end do
+      !$OMP END PARALLEL DO
+      if(evolve_tracer.ne.0) then
+         !$OMP PARALLEL DO PRIVATE(j)
+         do k = GRHydro_stencil, nz - GRHydro_stencil + 1
+            do j = GRHydro_stencil, nx - GRHydro_stencil + 1
+              call SimplePPM_tracer_1d(ny,CCTK_DELTA_SPACE(2),rho(j,:,k), &
+                   vely(j,:,k),velz(j,:,k),velx(j,:,k), &
+                   tracer(j,:,k,:),tracerminus(j,:,k,:),tracerplus(j,:,k,:), &
+                   press(j,:,k))
+            end do
+         end do
+         !$OMP END PARALLEL DO
+      end if
+
+    else if (flux_direction == 3) then
+      !$OMP PARALLEL DO PRIVATE(i, j)
+      do k = GRHydro_stencil, ny - GRHydro_stencil + 1
+        do j = GRHydro_stencil, nx - GRHydro_stencil + 1
+          call SimplePPM_1d(GRHydro_eos_handle,1,nz,CCTK_DELTA_SPACE(3),&
+               rho(j,k,:),velz(j,k,:),velx(j,k,:),vely(j,k,:),eps(j,k,:),&
+               press(j,k,:),rhominus(j,k,:),velzminus(j,k,:),velxminus(j,k,:),&
+               velyminus(j,k,:),epsminus(j,k,:),rhoplus(j,k,:),&
+               velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),epsplus(j,k,:),&
+               trivial_rp(j,k,:), hydro_excision_mask(j,k,:),&
+               gzz(j,k,:), gxz(j,k,:), gyz(j,k,:), gxx(j,k,:), gxy(j,k,:), &
+               gyy(j,k,:), psi4(j,k,:), lbetaz(j,k,:), alp(j,k,:),&
+               w_lorentz(j,k,:), &
+               3, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_z_left, &
+                                    GRHydro_mppm_eigenvalue_z_right, &
+                                    GRHydro_mppm_xwind)
+          do i = 1, nz
+            if (trivial_rp(j,k,i)) then
+              SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, trivialz)
+            else
+              SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, not_trivialz)
+            end if
+          end do
+        end do
+      end do
+      !$OMP END PARALLEL DO
+      if(evolve_tracer.ne.0) then
+         !$OMP PARALLEL DO PRIVATE(j)
+         do k = GRHydro_stencil, ny - GRHydro_stencil + 1
+            do j = GRHydro_stencil, nx - GRHydro_stencil + 1
+
+               call SimplePPM_tracer_1d(nz,CCTK_DELTA_SPACE(3),rho(j,k,:), &
+                    velz(j,k,:),velx(j,k,:),vely(j,k,:), &
+                    tracer(j,k,:,:),tracerminus(j,k,:,:),tracerplus(j,k,:,:), &
+                    press(j,k,:))
+            end do
+         end do
+         !$OMP END PARALLEL DO
+      end if
+
+    else
+      call CCTK_WARN(0, "Flux direction not x,y,z")
+    end if
+  else if (CCTK_EQUALS(recon_method,"eno")) then
+
+     if (evolve_tracer .ne. 0) then
+        do itracer=1,number_of_tracers
+           call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
+                tracer(:,:,:,itracer), tracerplus(:,:,:,itracer), &
+                tracerminus(:,:,:,itracer), trivial_rp, &
+                hydro_excision_mask)
+        enddo
+     end if
+
+    if (flux_direction == 1) then
+      !$OMP PARALLEL DO PRIVATE(i, j)
+      do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + 1
+        do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + 1
+          if (CCTK_EQUALS(recon_vars,"primitive")) then
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
+                 rho(:,j,k),rhominus(:,j,k),rhoplus(:,j,k),&
+                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
+                 velx(:,j,k),velxminus(:,j,k),velxplus(:,j,k),&
+                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
+                 vely(:,j,k),velyminus(:,j,k),velyplus(:,j,k),&
+                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
+                 velz(:,j,k),velzminus(:,j,k),velzplus(:,j,k),&
+                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
+          else if (CCTK_EQUALS(recon_vars,"conservative")) then
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
+                 dens(:,j,k),densminus(:,j,k),densplus(:,j,k),&
+                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
+                 sx(:,j,k),sxminus(:,j,k),sxplus(:,j,k),&
+                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
+                 sy(:,j,k),syminus(:,j,k),syplus(:,j,k),&
+                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
+                 sz(:,j,k),szminus(:,j,k),szplus(:,j,k),&
+                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
+          else
+            !$OMP CRITICAL
+            call CCTK_WARN(0, "Variable type to reconstruct not recognized.")
+            !$OMP END CRITICAL
+          end if
+          do i = 1, cctk_lsh(1)
+            if (trivial_rp(i,j,k)) then
+              SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, trivialx)
+            else
+              SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, not_trivialx)
+            end if
+          end do
+        end do
+      end do
+      !$OMP END PARALLEL DO
+    else if (flux_direction == 2) then
+      !$OMP PARALLEL DO PRIVATE(i, j)
+      do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + 1
+        do j = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + 1
+          if (CCTK_EQUALS(recon_vars,"primitive")) then
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
+                 rho(j,:,k),rhominus(j,:,k),rhoplus(j,:,k),&
+                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
+                 velx(j,:,k),velxminus(j,:,k),velxplus(j,:,k),&
+                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
+                 vely(j,:,k),velyminus(j,:,k),velyplus(j,:,k),&
+                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
+                 velz(j,:,k),velzminus(j,:,k),velzplus(j,:,k),&
+                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
+          else if (CCTK_EQUALS(recon_vars,"conservative")) then
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
+                 dens(j,:,k),densminus(j,:,k),densplus(j,:,k),&
+                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
+                 sx(j,:,k),sxminus(j,:,k),sxplus(j,:,k),&
+                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
+                 sy(j,:,k),syminus(j,:,k),syplus(j,:,k),&
+                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
+                 sz(j,:,k),szminus(j,:,k),szplus(j,:,k),&
+                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
+          else
+            !$OMP CRITICAL
+            call CCTK_WARN(0, "Variable type to reconstruct not recognized.")
+            !$OMP END CRITICAL
+          end if
+          do i = 1, cctk_lsh(2)
+            if (trivial_rp(j,i,k)) then
+              SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, trivialy)
+            else
+              SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, not_trivialy)
+            end if
+          end do
+        end do
+      end do
+      !$OMP END PARALLEL DO
+    else if (flux_direction == 3) then
+      !$OMP PARALLEL DO PRIVATE(i, j)
+      do k = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + 1
+        do j = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + 1
+          if (CCTK_EQUALS(recon_vars,"primitive")) then
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
+                 rho(j,k,:),rhominus(j,k,:),rhoplus(j,k,:),&
+                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
+                 velx(j,k,:),velxminus(j,k,:),velxplus(j,k,:),&
+                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
+                 vely(j,k,:),velyminus(j,k,:),velyplus(j,k,:),&
+                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
+                 velz(j,k,:),velzminus(j,k,:),velzplus(j,k,:),&
+                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
+          else if (CCTK_EQUALS(recon_vars,"conservative")) then
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
+                 dens(j,k,:),densminus(j,k,:),densplus(j,k,:),&
+                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
+                 sx(j,k,:),sxminus(j,k,:),sxplus(j,k,:),&
+                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
+                 sy(j,k,:),syminus(j,k,:),syplus(j,k,:),&
+                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
+                 sz(j,k,:),szminus(j,k,:),szplus(j,k,:),&
+                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
+          else
+            !$OMP CRITICAL
+            call CCTK_WARN(0, "Variable type to reconstruct not recognized.")
+            !$OMP END CRITICAL
+          end if
+          do i = 1, cctk_lsh(3)
+            if (trivial_rp(j,k,i)) then
+              SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, trivialz)
+            else
+              SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, not_trivialz)
+            end if
+          end do
+        end do
+      end do
+      !$OMP END PARALLEL DO
+    else
+      call CCTK_WARN(0, "Flux direction not x,y,z")
+    end if
+  else
+    call CCTK_WARN(0, "Reconstruction method not recognized!")
+  end if
+
+  deallocate(trivial_rp)
+  deallocate(psi4, lbetax, lbetay, lbetaz)
+
+  !$OMP WORKSHARE
+  where ( (rhoplus < GRHydro_rho_min).or.(rhominus < GRHydro_rho_min).or.&
+          (epsplus < 0.d0).or.(epsminus < 0.d0) )
+    rhoplus = rho
+    rhominus = rho
+    velxplus = vel(:,:,:,1)
+    velxminus = vel(:,:,:,1)
+    velyplus = vel(:,:,:,2)
+    velyminus = vel(:,:,:,2)
+    velzplus = vel(:,:,:,3)
+    velzminus = vel(:,:,:,3)
+    epsplus = eps
+    epsminus = eps
+  end where
+  !$OMP END WORKSHARE
+
+  if (evolve_tracer .ne. 0) then
+    if (use_min_tracer .ne. 0) then
+      local_min_tracer = min_tracer
+    else
+      local_min_tracer = 0d0
+    end if
+   
+    !$OMP WORKSHARE
+    where( (tracerplus  .le. local_min_tracer).or.&
+           (tracerminus .le. local_min_tracer) )
+      tracerplus = tracer
+      tracerminus = tracer
+    end where
+    !$OMP END WORKSHARE
+    ! Call the conserved tracer routine in any case because (accord. to
+    ! Christian Ott) this is the only way this works
+    call Prim2ConservativeTracer(CCTK_PASS_FTOF)
+  endif
+
+  if (CCTK_EQUALS(recon_vars,"primitive").or.&
+       CCTK_EQUALS(recon_method,"ppm")) then
+    if (use_eosgeneral == 0) then
+      call Prim2ConservativePolytype(CCTK_PASS_FTOF)
+    else
+      call primitive2conservativegeneral(CCTK_PASS_FTOF)
+    end if
+  else if (CCTK_EQUALS(recon_vars,"conservative")) then
+    call Con2PrimBoundsPolytype(CCTK_PASS_FTOF)
+  else
+    call CCTK_WARN(0,"Variable type to reconstruct not recognized.")
+  end if
+
+  return
+
+end subroutine ReconstructionPolytype
+