RIT GRMHD dev: split reconstruct routine.

Comments:

1) The reconstruct routine is the most expensive routine
in GRHydro and may require some serious thought in the
future regarding optimization.

2) This routine was merged by the end of the last year.
With too many loops, the compiler wasn't able to fully
optimize it (and it would warn about that).

3) With this split, the code became easier to read,
the compiler was able to fully optimize it, but we only
have a marginal gain in the routine performance (3.2%). 
Further optimization and clean-up is still desirable.



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

5 files changed, 1132 insertions, 784 deletions

diff --git a/src/GRHydro_ENOReconstruct_drv.F90 b/src/GRHydro_ENOReconstruct_drv.F90
new file mode 100644
index 0000000..656fd6e
--- /dev/null
+++ b/src/GRHydro_ENOReconstruct_drv.F90
@@ -0,0 +1,403 @@
+ /*@@
+   @file      GRHydro_ENOReconstruct_drv.F90
+   @date      Tue Jul 19 13:22:03 EDT 2011
+   @author    Bruno C. Mundim, Joshua Faber
+   @desc 
+   Driver routine to perform the ENO reconstruction.
+   @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)
+#define Bvecx(i,j,k) Bvec(i,j,k,1)
+#define Bvecy(i,j,k) Bvec(i,j,k,2)
+#define Bvecz(i,j,k) Bvec(i,j,k,3)
+#define Bconsx(i,j,k) Bcons(i,j,k,1)
+#define Bconsy(i,j,k) Bcons(i,j,k,2)
+#define Bconsz(i,j,k) Bcons(i,j,k,3)
+
+
+ /*@@
+   @routine    GRHydro_ENOReconstruct_drv
+   @date       Tue Jul 19 13:24:34 EDT 2011
+   @author     Luca Baiotti, Ian Hawke, Bruno C. Mundim, Joshua Faber
+   @desc 
+   A driver routine to do ENO reconstruction. 
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+subroutine GRHydro_ENOReconstruct_drv(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, dum, dump, dumm
+!!$  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
+  
+!!$ The dum variables are used as dummies if MHD on but divergence cleaning off
+  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)),&
+       dum(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),&
+       dump(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),&
+       dumm(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)
+
+  !$OMP PARALLEL
+  !$OMP WORKSHARE
+  trivial_rp = .false.
+  !$OMP END WORKSHARE NOWAIT
+
+!!$  Currently only option is reconstruction on primitive variables.
+!!$  Should change this.
+
+  !$OMP WORKSHARE
+  psi4 = 1.d0
+  !$OMP END WORKSHARE NOWAIT
+  if (shift_state .ne. 0) then
+    !$OMP WORKSHARE
+    lbetax = betax
+    lbetay = betay
+    lbetaz = betaz
+    !$OMP END WORKSHARE NOWAIT
+  else
+    !$OMP WORKSHARE
+    lbetax = 0.d0
+    lbetay = 0.d0
+    lbetaz = 0.d0
+    !$OMP END WORKSHARE NOWAIT
+  end if
+  !$OMP END PARALLEL
+
+!!$ ENO starts:
+     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 (evolve_Y_e .ne. 0) then
+        call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
+             Y_e(:,:,:), Y_e_plus(:,:,:), &
+             Y_e_minus(:,:,:), &
+             trivial_rp, hydro_excision_mask)
+     endif
+
+    if (flux_direction == 1) then
+      !$OMP PARALLEL DO PRIVATE(i, j, k)
+      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))
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
+                 eps(:,j,k),epsminus(:,j,k),epsplus(:,j,k),&
+                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
+            if(evolve_mhd.ne.0) then
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
+                    Bvecx(:,j,k),Bvecxminus(:,j,k),Bvecxplus(:,j,k),&
+                    trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
+                    Bvecy(:,j,k),Bvecyminus(:,j,k),Bvecyplus(:,j,k),&
+                    trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
+                    Bvecz(:,j,k),Bveczminus(:,j,k),Bveczplus(:,j,k),&
+                    trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
+            endif
+          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))
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
+                 tau(:,j,k),tauminus(:,j,k),tauplus(:,j,k),&
+                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
+            if(evolve_mhd.ne.0) then
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
+                    Bconsx(:,j,k),Bconsxminus(:,j,k),Bconsxplus(:,j,k),&
+                    trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
+                    Bconsy(:,j,k),Bconsyminus(:,j,k),Bconsyplus(:,j,k),&
+                    trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
+                    Bconsz(:,j,k),Bconszminus(:,j,k),Bconszplus(:,j,k),&
+                    trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
+            endif
+          else
+            !$OMP CRITICAL
+            call CCTK_WARN(0, "Variable type to reconstruct not recognized.")
+            !$OMP END CRITICAL
+          end if
+
+          if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then
+             call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
+                  psidc(:,j,k),psidcminus(:,j,k),psidcplus(:,j,k),&
+                  trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
+          endif
+
+          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, k)
+      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))
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
+                 eps(j,:,k),epsminus(j,:,k),epsplus(j,:,k),&
+                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
+            if(evolve_mhd.ne.0) then
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
+                    Bvecx(j,:,k),Bvecxminus(j,:,k),Bvecxplus(j,:,k),&
+                    trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
+                    Bvecy(j,:,k),Bvecyminus(j,:,k),Bvecyplus(j,:,k),&
+                    trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
+                    Bvecz(j,:,k),Bveczminus(j,:,k),Bveczplus(j,:,k),&
+                    trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
+            endif
+         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))
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
+                 tau(j,:,k),tauminus(j,:,k),tauplus(j,:,k),&
+                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
+            if(evolve_mhd.ne.0) then
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
+                    Bconsx(j,:,k),Bconsxminus(j,:,k),Bconsxplus(j,:,k),&
+                    trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
+                    Bconsy(j,:,k),Bconsyminus(j,:,k),Bconsyplus(j,:,k),&
+                    trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
+                    Bconsz(j,:,k),Bconszminus(j,:,k),Bconszplus(j,:,k),&
+                    trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
+            endif
+          else
+            !$OMP CRITICAL
+            call CCTK_WARN(0, "Variable type to reconstruct not recognized.")
+            !$OMP END CRITICAL
+          end if
+
+          if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
+                 psidc(j,:,k),psidcminus(j,:,k),psidcplus(j,:,k),&
+                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
+          endif
+
+          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, k)
+      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,:))
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
+                 eps(j,k,:),epsminus(j,k,:),epsplus(j,k,:),&
+                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
+            if(evolve_mhd.ne.0) then
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
+                    Bvecx(j,k,:),Bvecxminus(j,k,:),Bvecxplus(j,k,:),&
+                    trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
+                    Bvecy(j,k,:),Bvecyminus(j,k,:),Bvecyplus(j,k,:),&
+                    trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
+                    Bvecz(j,k,:),Bveczminus(j,k,:),Bveczplus(j,k,:),&
+                    trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
+            endif
+          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,:))
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
+                 tau(j,k,:),tauminus(j,k,:),tauplus(j,k,:),&
+                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
+            if(evolve_mhd.ne.0) then
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
+                    Bconsx(j,k,:),Bconsxminus(j,k,:),Bconsxplus(j,k,:),&
+                    trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
+                    Bconsy(j,k,:),Bconsyminus(j,k,:),Bconsyplus(j,k,:),&
+                    trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
+                    Bconsz(j,k,:),Bconszminus(j,k,:),Bconszplus(j,k,:),&
+                    trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
+            endif
+          else
+            !$OMP CRITICAL
+            call CCTK_WARN(0, "Variable type to reconstruct not recognized.")
+            !$OMP END CRITICAL
+          end if
+
+          if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then
+             call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
+                  psidc(j,k,:),psidcminus(j,k,:),psidcplus(j,k,:),&
+                  trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
+          endif
+
+          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
+!!$ ENO ends.
+
+  deallocate(trivial_rp)
+  deallocate(psi4, lbetax, lbetay, lbetaz)
+  deallocate(dum,dump,dumm)
+
+end subroutine GRHydro_ENOReconstruct_drv
+
diff --git a/src/GRHydro_PPMReconstruct_drv.F90 b/src/GRHydro_PPMReconstruct_drv.F90
new file mode 100644
index 0000000..ddd68cd
--- /dev/null
+++ b/src/GRHydro_PPMReconstruct_drv.F90
@@ -0,0 +1,482 @@
+ /*@@
+   @file      GRHydro_PPMReconstruct_drv.F90
+   @date      Tue Jul 19 13:22:03 EDT 2011
+   @author    Bruno C. Mundim, Joshua Faber
+   @desc 
+   Driver routine to perform the PPM reconstructions.
+   @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)
+#define Bvecx(i,j,k) Bvec(i,j,k,1)
+#define Bvecy(i,j,k) Bvec(i,j,k,2)
+#define Bvecz(i,j,k) Bvec(i,j,k,3)
+#define Bconsx(i,j,k) Bcons(i,j,k,1)
+#define Bconsy(i,j,k) Bcons(i,j,k,2)
+#define Bconsz(i,j,k) Bcons(i,j,k,3)
+
+
+ /*@@
+   @routine    GRHydro_PPMReconstruct_drv
+   @date       Tue Jul 19 13:24:34 EDT 2011
+   @author     Luca Baiotti, Ian Hawke, Bruno C. Mundim, Joshua Faber
+   @desc 
+   A driver routine to do PPM reconstructions. 
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+subroutine GRHydro_PPMReconstruct_drv(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, dum, dump, dumm
+!!$  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
+  
+!!$ The dum variables are used as dummies if MHD on but divergence cleaning off
+  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)),&
+       dum(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),&
+       dump(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),&
+       dumm(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)
+
+  !$OMP PARALLEL
+  !$OMP WORKSHARE
+  trivial_rp = .false.
+  !$OMP END WORKSHARE NOWAIT
+
+!!$  Currently only option is reconstruction on primitive variables.
+!!$  Should change this.
+
+  !$OMP WORKSHARE
+  psi4 = 1.d0
+  !$OMP END WORKSHARE NOWAIT
+  if (shift_state .ne. 0) then
+    !$OMP WORKSHARE
+    lbetax = betax
+    lbetay = betay
+    lbetaz = betaz
+    !$OMP END WORKSHARE NOWAIT
+  else
+    !$OMP WORKSHARE
+    lbetax = 0.d0
+    lbetay = 0.d0
+    lbetaz = 0.d0
+    !$OMP END WORKSHARE NOWAIT
+  end if
+  !$OMP END PARALLEL
+
+
+!!$ PPM starts:
+    if (flux_direction == 1) then
+       if(evolve_mhd.ne.0) then
+          if(clean_divergence.ne.0) then
+            !$OMP PARALLEL DO PRIVATE(i, j, k)
+            do k = GRHydro_stencil, nz - GRHydro_stencil + 1
+              do j = GRHydro_stencil, ny - GRHydro_stencil + 1
+                 call SimplePPM_1dM(GRHydro_eos_handle,0,nx,CCTK_DELTA_SPACE(1),&
+                      rho(:,j,k),velx(:,j,k),vely(:,j,k),velz(:,j,k),&
+                      Bvecx(:,j,k),Bvecy(:,j,k),Bvecz(:,j,k),psidc(:,j,k),eps(:,j,k),press(:,j,k),&
+                      rhominus(:,j,k),velxminus(:,j,k),velyminus(:,j,k),velzminus(:,j,k),&
+                      Bvecxminus(:,j,k),Bvecyminus(:,j,k),Bveczminus(:,j,k),psidcminus(:,j,k),epsminus(:,j,k),&
+                      rhoplus(:,j,k),velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),&
+                      Bvecxplus(:,j,k),Bvecyplus(:,j,k),Bveczplus(:,j,k),psidcplus(:,j,k),epsplus(:,j,k),&
+                      1,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
+          else  !clean_divergence
+            !$OMP PARALLEL DO PRIVATE(i, j, k)
+            do k = GRHydro_stencil, nz - GRHydro_stencil + 1
+              do j = GRHydro_stencil, ny - GRHydro_stencil + 1
+                 call SimplePPM_1dM(GRHydro_eos_handle,0,nx,CCTK_DELTA_SPACE(1),&
+                      rho(:,j,k),velx(:,j,k),vely(:,j,k),velz(:,j,k),&
+                      Bvecx(:,j,k),Bvecy(:,j,k),Bvecz(:,j,k),dum(:,j,k),eps(:,j,k),press(:,j,k),&
+                      rhominus(:,j,k),velxminus(:,j,k),velyminus(:,j,k),velzminus(:,j,k),&
+                      Bvecxminus(:,j,k),Bvecyminus(:,j,k),Bveczminus(:,j,k),dumm(:,j,k),epsminus(:,j,k),&
+                      rhoplus(:,j,k),velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),&
+                      Bvecxplus(:,j,k),Bvecyplus(:,j,k),Bveczplus(:,j,k),dump(:,j,k),epsplus(:,j,k),&
+                      0,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
+          endif  !clean_divergence
+       else  !evolve_mhd
+         !$OMP PARALLEL DO PRIVATE(i, j, k)
+         do k = GRHydro_stencil, nz - GRHydro_stencil + 1
+           do j = GRHydro_stencil, ny - GRHydro_stencil + 1
+           call SimplePPM_1d(GRHydro_eos_handle,0,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
+       endif  !evolve_mhd
+      if(evolve_tracer.ne.0) then
+         !$OMP PARALLEL DO PRIVATE(j, k)
+         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
+      if(evolve_Y_e.ne.0) then
+         !$OMP PARALLEL DO PRIVATE(j, k)
+         do k = GRHydro_stencil, nz - GRHydro_stencil + 1
+            do j = GRHydro_stencil, ny - GRHydro_stencil + 1
+               call SimplePPM_ye_1d(nx,CCTK_DELTA_SPACE(1),rho(:,j,k), &
+                    velx(:,j,k),vely(:,j,k),velz(:,j,k), &
+                    Y_e(:,j,k),Y_e_minus(:,j,k),Y_e_plus(:,j,k), &
+                    press(:,j,k))
+            end do
+         end do
+         !$OMP END PARALLEL DO
+      end if
+
+    else if (flux_direction == 2) then
+      if(evolve_mhd.ne.0) then
+        if(clean_divergence.ne.0) then
+          !$OMP PARALLEL DO PRIVATE(i, j, k)
+          do k = GRHydro_stencil, nz - GRHydro_stencil + 1
+            do j = GRHydro_stencil, nx - GRHydro_stencil + 1
+              call SimplePPM_1dM(GRHydro_eos_handle,0,ny,CCTK_DELTA_SPACE(2),&
+                   rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),&
+                   Bvecy(j,:,k),Bvecz(j,:,k),Bvecx(j,:,k),psidc(j,:,k),eps(j,:,k),press(j,:,k),&
+                   rhominus(j,:,k),velyminus(j,:,k),velzminus(j,:,k),velxminus(j,:,k),&
+                   Bvecyminus(j,:,k),Bveczminus(j,:,k),Bvecxminus(j,:,k),psidcminus(j,:,k),epsminus(j,:,k),&
+                   rhoplus(j,:,k),velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),&
+                   Bvecyplus(j,:,k),Bveczplus(j,:,k),Bvecxplus(j,:,k),psidcplus(j,:,k),epsplus(j,:,k),&
+                   1,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
+        else  !clean_divergence
+          !$OMP PARALLEL DO PRIVATE(i, j, k)
+          do k = GRHydro_stencil, nz - GRHydro_stencil + 1
+            do j = GRHydro_stencil, nx - GRHydro_stencil + 1
+              call SimplePPM_1dM(GRHydro_eos_handle,0,ny,CCTK_DELTA_SPACE(2),&
+                   rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),&
+                   Bvecy(j,:,k),Bvecz(j,:,k),Bvecx(j,:,k),dum(j,:,k),eps(j,:,k),press(j,:,k),&
+                   rhominus(j,:,k),velyminus(j,:,k),velzminus(j,:,k),velxminus(j,:,k),&
+                   Bvecyminus(j,:,k),Bveczminus(j,:,k),Bvecxminus(j,:,k),dumm(j,:,k),epsminus(j,:,k),&
+                   rhoplus(j,:,k),velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),&
+                   Bvecyplus(j,:,k),Bveczplus(j,:,k),Bvecxplus(j,:,k),dump(j,:,k),epsplus(j,:,k),&
+                   0,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
+        endif  !clean_divergence
+      else  !evolve_mhd
+        !$OMP PARALLEL DO PRIVATE(i, j, k)
+        do k = GRHydro_stencil, nz - GRHydro_stencil + 1
+          do j = GRHydro_stencil, nx - GRHydro_stencil + 1
+            call SimplePPM_1d(GRHydro_eos_handle,0,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
+      endif  !evolve_mhd
+      if(evolve_tracer.ne.0) then
+         !$OMP PARALLEL DO PRIVATE(j, k)
+         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
+      if(evolve_Y_e.ne.0) then
+         !$OMP PARALLEL DO PRIVATE(j, k)
+         do k = GRHydro_stencil, nz - GRHydro_stencil + 1
+            do j = GRHydro_stencil, nx - GRHydro_stencil + 1
+               call SimplePPM_ye_1d(ny,CCTK_DELTA_SPACE(2),rho(j,:,k), &
+                    velx(j,:,k),vely(j,:,k),velz(j,:,k), &
+                    Y_e(j,:,k),Y_e_minus(j,:,k),Y_e_plus(j,:,k), &
+                    press(j,:,k))
+            end do
+         end do
+         !$OMP END PARALLEL DO
+      end if
+
+    else if (flux_direction == 3) then
+      if(evolve_mhd.ne.0) then
+        if(clean_divergence.ne.0) then
+          !$OMP PARALLEL DO PRIVATE(i, j, k)
+          do k = GRHydro_stencil, ny - GRHydro_stencil + 1
+            do j = GRHydro_stencil, nx - GRHydro_stencil + 1
+              call SimplePPM_1dM(GRHydro_eos_handle,0,nz,CCTK_DELTA_SPACE(3),&
+                   rho(j,k,:),velz(j,k,:),velx(j,k,:),vely(j,k,:),&
+                   Bvecz(j,k,:),Bvecx(j,k,:),Bvecy(j,k,:),psidc(j,k,:),eps(j,k,:),press(j,k,:),&
+                   rhominus(j,k,:),velzminus(j,k,:),velxminus(j,k,:),velyminus(j,k,:),&
+                   Bveczminus(j,k,:),Bvecxminus(j,k,:),Bvecyminus(j,k,:),psidcminus(j,k,:),epsminus(j,k,:),&
+                   rhoplus(j,k,:),velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),&
+                   Bveczplus(j,k,:),Bvecxplus(j,k,:),Bvecyplus(j,k,:),psidcplus(j,k,:),epsplus(j,k,:),&
+                   1,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
+        else  !clean_divergence
+          !$OMP PARALLEL DO PRIVATE(i, j, k)
+          do k = GRHydro_stencil, ny - GRHydro_stencil + 1
+            do j = GRHydro_stencil, nx - GRHydro_stencil + 1
+              call SimplePPM_1dM(GRHydro_eos_handle,0,nz,CCTK_DELTA_SPACE(3),&
+                   rho(j,k,:),velz(j,k,:),velx(j,k,:),vely(j,k,:),&
+                   Bvecz(j,k,:),Bvecx(j,k,:),Bvecy(j,k,:),dum(j,k,:),eps(j,k,:),press(j,k,:),&
+                   rhominus(j,k,:),velzminus(j,k,:),velxminus(j,k,:),velyminus(j,k,:),&
+                   Bveczminus(j,k,:),Bvecxminus(j,k,:),Bvecyminus(j,k,:),dumm(j,k,:),epsminus(j,k,:),&
+                   rhoplus(j,k,:),velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),&
+                   Bveczplus(j,k,:),Bvecxplus(j,k,:),Bvecyplus(j,k,:),dump(j,k,:),epsplus(j,k,:),&
+                   0,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
+        endif  !clean_divergence
+      else  !evolve_mhd
+        !$OMP PARALLEL DO PRIVATE(i, j, k)
+        do k = GRHydro_stencil, ny - GRHydro_stencil + 1
+          do j = GRHydro_stencil, nx - GRHydro_stencil + 1
+            call SimplePPM_1d(GRHydro_eos_handle,0,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
+      endif  !evolve_mhd
+      if(evolve_tracer.ne.0) then
+         !$OMP PARALLEL DO PRIVATE(j, k)
+         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
+      if(evolve_Y_e.ne.0) then
+         !$OMP PARALLEL DO PRIVATE(j, k)
+         do k = GRHydro_stencil, ny - GRHydro_stencil + 1
+            do j = GRHydro_stencil, nx - GRHydro_stencil + 1
+               call SimplePPM_ye_1d(nz,CCTK_DELTA_SPACE(3),rho(j,k,:), &
+                    velx(j,k,:),vely(j,k,:),velz(j,k,:), &
+                    Y_e(j,k,:),Y_e_minus(j,k,:),Y_e_plus(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
+!!$ PPM ends.
+
+  deallocate(trivial_rp)
+  deallocate(psi4, lbetax, lbetay, lbetaz)
+  deallocate(dum,dump,dumm)
+
+end subroutine GRHydro_PPMReconstruct_drv
+
diff --git a/src/GRHydro_Reconstruct.F90 b/src/GRHydro_Reconstruct.F90
index 5b447c0..a8ac2ec 100644
--- a/src/GRHydro_Reconstruct.F90
+++ b/src/GRHydro_Reconstruct.F90
@@ -14,20 +14,6 @@
 
 #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)
-#define Bvecx(i,j,k) Bvec(i,j,k,1)
-#define Bvecy(i,j,k) Bvec(i,j,k,2)
-#define Bvecz(i,j,k) Bvec(i,j,k,3)
-#define Bconsx(i,j,k) Bcons(i,j,k,1)
-#define Bconsy(i,j,k) Bcons(i,j,k,2)
-#define Bconsz(i,j,k) Bcons(i,j,k,3)
-
-
  /*@@
    @routine    Reconstruction
    @date       Sat Jan 26 02:13:47 2002
@@ -52,89 +38,11 @@ subroutine Reconstruction(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, dum, dump, dumm
-!!$  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
-  
-!!$ The dum variables are used as dummies if MHD on but divergence cleaning off
-  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)),&
-       dum(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),&
-       dump(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),&
-       dumm(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)
-
-  !$OMP PARALLEL
-  !$OMP WORKSHARE
-  trivial_rp = .false.
-  !$OMP END WORKSHARE NOWAIT
-
-!!$  Currently only option is reconstruction on primitive variables.
-!!$  Should change this.
-
-  !$OMP WORKSHARE
-  psi4 = 1.d0
-  !$OMP END WORKSHARE NOWAIT
-  if (shift_state .ne. 0) then
-    !$OMP WORKSHARE
-    lbetax = betax
-    lbetay = betay
-    lbetaz = betaz
-    !$OMP END WORKSHARE NOWAIT
-  else
-    !$OMP WORKSHARE
-    lbetax = 0.d0
-    lbetay = 0.d0
-    lbetaz = 0.d0
-    !$OMP END WORKSHARE NOWAIT
-  end if
-
 !!$ Initialize variables that store reconstructed quantities
 
+  !$OMP PARALLEL
   !$OMP WORKSHARE
   rhoplus = 0.0d0
   rhominus = 0.0d0
@@ -183,706 +91,20 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
 
   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 (evolve_Y_e .ne. 0) then
-        call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-             Y_e(:,:,:), Y_e_plus(:,:,:), &
-             Y_e_minus(:,:,:), &
-             trivial_rp, hydro_excision_mask)
-     endif
-
-    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)
-      call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-           eps, epsplus, epsminus, trivial_rp, hydro_excision_mask)
-      if(evolve_mhd.ne.0) then
-         call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-              Bvec(:,:,:,1), Bvecxplus, Bvecxminus, trivial_rp, hydro_excision_mask)
-         call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-              Bvec(:,:,:,2), Bvecyplus, Bvecyminus, trivial_rp, hydro_excision_mask)
-         call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-              Bvec(:,:,:,3), Bveczplus, Bveczminus, trivial_rp, hydro_excision_mask)
-      endif
-
-    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)
-      call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-           tau, tauplus, tauminus, trivial_rp, hydro_excision_mask)
-      if(evolve_mhd.ne.0) then
-         call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-              Bcons(:,:,:,1), Bconsxplus, Bconsxminus, trivial_rp, hydro_excision_mask)
-         call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-              Bcons(:,:,:,2), Bconsyplus, Bconsyminus, trivial_rp, hydro_excision_mask)
-         call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-              Bcons(:,:,:,3), Bconszplus, Bconszminus, trivial_rp, hydro_excision_mask)
-      endif
-
-    else
-      call CCTK_WARN(0, "Variable type to reconstruct not recognized.")
-    end if
-
-      if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then
-       call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-            psidc, psidcplus, psidcminus, trivial_rp, hydro_excision_mask)
-    endif
-
-    !$OMP PARALLEL DO PRIVATE(i, j, k)
-    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
+    call GRHydro_TVDReconstruct_drv(CCTK_PASS_FTOF)
 
   else if (CCTK_EQUALS(recon_method,"ppm")) then
 
-    if (flux_direction == 1) then
-       if(evolve_mhd.ne.0) then
-          if(clean_divergence.ne.0) then
-            !$OMP PARALLEL DO PRIVATE(i, j, k)
-            do k = GRHydro_stencil, nz - GRHydro_stencil + 1
-              do j = GRHydro_stencil, ny - GRHydro_stencil + 1
-                 call SimplePPM_1dM(GRHydro_eos_handle,0,nx,CCTK_DELTA_SPACE(1),&
-                      rho(:,j,k),velx(:,j,k),vely(:,j,k),velz(:,j,k),&
-                      Bvecx(:,j,k),Bvecy(:,j,k),Bvecz(:,j,k),psidc(:,j,k),eps(:,j,k),press(:,j,k),&
-                      rhominus(:,j,k),velxminus(:,j,k),velyminus(:,j,k),velzminus(:,j,k),&
-                      Bvecxminus(:,j,k),Bvecyminus(:,j,k),Bveczminus(:,j,k),psidcminus(:,j,k),epsminus(:,j,k),&
-                      rhoplus(:,j,k),velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),&
-                      Bvecxplus(:,j,k),Bvecyplus(:,j,k),Bveczplus(:,j,k),psidcplus(:,j,k),epsplus(:,j,k),&
-                      1,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
-          else  !clean_divergence
-            !$OMP PARALLEL DO PRIVATE(i, j, k)
-            do k = GRHydro_stencil, nz - GRHydro_stencil + 1
-              do j = GRHydro_stencil, ny - GRHydro_stencil + 1
-                 call SimplePPM_1dM(GRHydro_eos_handle,0,nx,CCTK_DELTA_SPACE(1),&
-                      rho(:,j,k),velx(:,j,k),vely(:,j,k),velz(:,j,k),&
-                      Bvecx(:,j,k),Bvecy(:,j,k),Bvecz(:,j,k),dum(:,j,k),eps(:,j,k),press(:,j,k),&
-                      rhominus(:,j,k),velxminus(:,j,k),velyminus(:,j,k),velzminus(:,j,k),&
-                      Bvecxminus(:,j,k),Bvecyminus(:,j,k),Bveczminus(:,j,k),dumm(:,j,k),epsminus(:,j,k),&
-                      rhoplus(:,j,k),velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),&
-                      Bvecxplus(:,j,k),Bvecyplus(:,j,k),Bveczplus(:,j,k),dump(:,j,k),epsplus(:,j,k),&
-                      0,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
-          endif  !clean_divergence
-       else  !evolve_mhd
-         !$OMP PARALLEL DO PRIVATE(i, j, k)
-         do k = GRHydro_stencil, nz - GRHydro_stencil + 1
-           do j = GRHydro_stencil, ny - GRHydro_stencil + 1
-           call SimplePPM_1d(GRHydro_eos_handle,0,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
-       endif  !evolve_mhd
-      if(evolve_tracer.ne.0) then
-         !$OMP PARALLEL DO PRIVATE(j, k)
-         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
-      if(evolve_Y_e.ne.0) then
-         !$OMP PARALLEL DO PRIVATE(j, k)
-         do k = GRHydro_stencil, nz - GRHydro_stencil + 1
-            do j = GRHydro_stencil, ny - GRHydro_stencil + 1
-               call SimplePPM_ye_1d(nx,CCTK_DELTA_SPACE(1),rho(:,j,k), &
-                    velx(:,j,k),vely(:,j,k),velz(:,j,k), &
-                    Y_e(:,j,k),Y_e_minus(:,j,k),Y_e_plus(:,j,k), &
-                    press(:,j,k))
-            end do
-         end do
-         !$OMP END PARALLEL DO
-      end if
-
-    else if (flux_direction == 2) then
-      if(evolve_mhd.ne.0) then
-        if(clean_divergence.ne.0) then
-          !$OMP PARALLEL DO PRIVATE(i, j, k)
-          do k = GRHydro_stencil, nz - GRHydro_stencil + 1
-            do j = GRHydro_stencil, nx - GRHydro_stencil + 1
-              call SimplePPM_1dM(GRHydro_eos_handle,0,ny,CCTK_DELTA_SPACE(2),&
-                   rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),&
-                   Bvecy(j,:,k),Bvecz(j,:,k),Bvecx(j,:,k),psidc(j,:,k),eps(j,:,k),press(j,:,k),&
-                   rhominus(j,:,k),velyminus(j,:,k),velzminus(j,:,k),velxminus(j,:,k),&
-                   Bvecyminus(j,:,k),Bveczminus(j,:,k),Bvecxminus(j,:,k),psidcminus(j,:,k),epsminus(j,:,k),&
-                   rhoplus(j,:,k),velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),&
-                   Bvecyplus(j,:,k),Bveczplus(j,:,k),Bvecxplus(j,:,k),psidcplus(j,:,k),epsplus(j,:,k),&
-                   1,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
-        else  !clean_divergence
-          !$OMP PARALLEL DO PRIVATE(i, j, k)
-          do k = GRHydro_stencil, nz - GRHydro_stencil + 1
-            do j = GRHydro_stencil, nx - GRHydro_stencil + 1
-              call SimplePPM_1dM(GRHydro_eos_handle,0,ny,CCTK_DELTA_SPACE(2),&
-                   rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),&
-                   Bvecy(j,:,k),Bvecz(j,:,k),Bvecx(j,:,k),dum(j,:,k),eps(j,:,k),press(j,:,k),&
-                   rhominus(j,:,k),velyminus(j,:,k),velzminus(j,:,k),velxminus(j,:,k),&
-                   Bvecyminus(j,:,k),Bveczminus(j,:,k),Bvecxminus(j,:,k),dumm(j,:,k),epsminus(j,:,k),&
-                   rhoplus(j,:,k),velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),&
-                   Bvecyplus(j,:,k),Bveczplus(j,:,k),Bvecxplus(j,:,k),dump(j,:,k),epsplus(j,:,k),&
-                   0,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
-        endif  !clean_divergence
-      else  !evolve_mhd
-        !$OMP PARALLEL DO PRIVATE(i, j, k)
-        do k = GRHydro_stencil, nz - GRHydro_stencil + 1
-          do j = GRHydro_stencil, nx - GRHydro_stencil + 1
-            call SimplePPM_1d(GRHydro_eos_handle,0,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
-      endif  !evolve_mhd
-      if(evolve_tracer.ne.0) then
-         !$OMP PARALLEL DO PRIVATE(j, k)
-         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
-      if(evolve_Y_e.ne.0) then
-         !$OMP PARALLEL DO PRIVATE(j, k)
-         do k = GRHydro_stencil, nz - GRHydro_stencil + 1
-            do j = GRHydro_stencil, nx - GRHydro_stencil + 1
-               call SimplePPM_ye_1d(ny,CCTK_DELTA_SPACE(2),rho(j,:,k), &
-                    velx(j,:,k),vely(j,:,k),velz(j,:,k), &
-                    Y_e(j,:,k),Y_e_minus(j,:,k),Y_e_plus(j,:,k), &
-                    press(j,:,k))
-            end do
-         end do
-         !$OMP END PARALLEL DO
-      end if
+    call GRHydro_PPMReconstruct_drv(CCTK_PASS_FTOF)
 
-    else if (flux_direction == 3) then
-      if(evolve_mhd.ne.0) then
-        if(clean_divergence.ne.0) then
-          !$OMP PARALLEL DO PRIVATE(i, j, k)
-          do k = GRHydro_stencil, ny - GRHydro_stencil + 1
-            do j = GRHydro_stencil, nx - GRHydro_stencil + 1
-              call SimplePPM_1dM(GRHydro_eos_handle,0,nz,CCTK_DELTA_SPACE(3),&
-                   rho(j,k,:),velz(j,k,:),velx(j,k,:),vely(j,k,:),&
-                   Bvecz(j,k,:),Bvecx(j,k,:),Bvecy(j,k,:),psidc(j,k,:),eps(j,k,:),press(j,k,:),&
-                   rhominus(j,k,:),velzminus(j,k,:),velxminus(j,k,:),velyminus(j,k,:),&
-                   Bveczminus(j,k,:),Bvecxminus(j,k,:),Bvecyminus(j,k,:),psidcminus(j,k,:),epsminus(j,k,:),&
-                   rhoplus(j,k,:),velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),&
-                   Bveczplus(j,k,:),Bvecxplus(j,k,:),Bvecyplus(j,k,:),psidcplus(j,k,:),epsplus(j,k,:),&
-                   1,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
-        else  !clean_divergence
-          !$OMP PARALLEL DO PRIVATE(i, j, k)
-          do k = GRHydro_stencil, ny - GRHydro_stencil + 1
-            do j = GRHydro_stencil, nx - GRHydro_stencil + 1
-              call SimplePPM_1dM(GRHydro_eos_handle,0,nz,CCTK_DELTA_SPACE(3),&
-                   rho(j,k,:),velz(j,k,:),velx(j,k,:),vely(j,k,:),&
-                   Bvecz(j,k,:),Bvecx(j,k,:),Bvecy(j,k,:),dum(j,k,:),eps(j,k,:),press(j,k,:),&
-                   rhominus(j,k,:),velzminus(j,k,:),velxminus(j,k,:),velyminus(j,k,:),&
-                   Bveczminus(j,k,:),Bvecxminus(j,k,:),Bvecyminus(j,k,:),dumm(j,k,:),epsminus(j,k,:),&
-                   rhoplus(j,k,:),velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),&
-                   Bveczplus(j,k,:),Bvecxplus(j,k,:),Bvecyplus(j,k,:),dump(j,k,:),epsplus(j,k,:),&
-                   0,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
-        endif  !clean_divergence
-      else  !evolve_mhd
-        !$OMP PARALLEL DO PRIVATE(i, j, k)
-        do k = GRHydro_stencil, ny - GRHydro_stencil + 1
-          do j = GRHydro_stencil, nx - GRHydro_stencil + 1
-            call SimplePPM_1d(GRHydro_eos_handle,0,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
-      endif  !evolve_mhd
-      if(evolve_tracer.ne.0) then
-         !$OMP PARALLEL DO PRIVATE(j, k)
-         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
-      if(evolve_Y_e.ne.0) then
-         !$OMP PARALLEL DO PRIVATE(j, k)
-         do k = GRHydro_stencil, ny - GRHydro_stencil + 1
-            do j = GRHydro_stencil, nx - GRHydro_stencil + 1
-               call SimplePPM_ye_1d(nz,CCTK_DELTA_SPACE(3),rho(j,k,:), &
-                    velx(j,k,:),vely(j,k,:),velz(j,k,:), &
-                    Y_e(j,k,:),Y_e_minus(j,k,:),Y_e_plus(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 (evolve_Y_e .ne. 0) then
-        call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-             Y_e(:,:,:), Y_e_plus(:,:,:), &
-             Y_e_minus(:,:,:), &
-             trivial_rp, hydro_excision_mask)
-     endif
-
-    if (flux_direction == 1) then
-      !$OMP PARALLEL DO PRIVATE(i, j, k)
-      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))
-            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
-                 eps(:,j,k),epsminus(:,j,k),epsplus(:,j,k),&
-                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
-            if(evolve_mhd.ne.0) then
-               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
-                    Bvecx(:,j,k),Bvecxminus(:,j,k),Bvecxplus(:,j,k),&
-                    trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
-               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
-                    Bvecy(:,j,k),Bvecyminus(:,j,k),Bvecyplus(:,j,k),&
-                    trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
-               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
-                    Bvecz(:,j,k),Bveczminus(:,j,k),Bveczplus(:,j,k),&
-                    trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
-            endif
-          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))
-            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
-                 tau(:,j,k),tauminus(:,j,k),tauplus(:,j,k),&
-                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
-            if(evolve_mhd.ne.0) then
-               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
-                    Bconsx(:,j,k),Bconsxminus(:,j,k),Bconsxplus(:,j,k),&
-                    trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
-               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
-                    Bconsy(:,j,k),Bconsyminus(:,j,k),Bconsyplus(:,j,k),&
-                    trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
-               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
-                    Bconsz(:,j,k),Bconszminus(:,j,k),Bconszplus(:,j,k),&
-                    trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
-            endif
-          else
-            !$OMP CRITICAL
-            call CCTK_WARN(0, "Variable type to reconstruct not recognized.")
-            !$OMP END CRITICAL
-          end if
-
-          if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then
-             call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
-                  psidc(:,j,k),psidcminus(:,j,k),psidcplus(:,j,k),&
-                  trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
-          endif
-
-          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, k)
-      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))
-            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
-                 eps(j,:,k),epsminus(j,:,k),epsplus(j,:,k),&
-                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
-            if(evolve_mhd.ne.0) then
-               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
-                    Bvecx(j,:,k),Bvecxminus(j,:,k),Bvecxplus(j,:,k),&
-                    trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
-               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
-                    Bvecy(j,:,k),Bvecyminus(j,:,k),Bvecyplus(j,:,k),&
-                    trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
-               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
-                    Bvecz(j,:,k),Bveczminus(j,:,k),Bveczplus(j,:,k),&
-                    trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
-            endif
-         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))
-            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
-                 tau(j,:,k),tauminus(j,:,k),tauplus(j,:,k),&
-                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
-            if(evolve_mhd.ne.0) then
-               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
-                    Bconsx(j,:,k),Bconsxminus(j,:,k),Bconsxplus(j,:,k),&
-                    trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
-               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
-                    Bconsy(j,:,k),Bconsyminus(j,:,k),Bconsyplus(j,:,k),&
-                    trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
-               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
-                    Bconsz(j,:,k),Bconszminus(j,:,k),Bconszplus(j,:,k),&
-                    trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
-            endif
-          else
-            !$OMP CRITICAL
-            call CCTK_WARN(0, "Variable type to reconstruct not recognized.")
-            !$OMP END CRITICAL
-          end if
+    call GRHydro_ENOReconstruct_drv(CCTK_PASS_FTOF)
 
-          if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then
-            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
-                 psidc(j,:,k),psidcminus(j,:,k),psidcplus(j,:,k),&
-                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
-          endif
-
-          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, k)
-      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,:))
-            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
-                 eps(j,k,:),epsminus(j,k,:),epsplus(j,k,:),&
-                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
-            if(evolve_mhd.ne.0) then
-               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
-                    Bvecx(j,k,:),Bvecxminus(j,k,:),Bvecxplus(j,k,:),&
-                    trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
-               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
-                    Bvecy(j,k,:),Bvecyminus(j,k,:),Bvecyplus(j,k,:),&
-                    trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
-               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
-                    Bvecz(j,k,:),Bveczminus(j,k,:),Bveczplus(j,k,:),&
-                    trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
-            endif
-          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,:))
-            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
-                 tau(j,k,:),tauminus(j,k,:),tauplus(j,k,:),&
-                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
-            if(evolve_mhd.ne.0) then
-               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
-                    Bconsx(j,k,:),Bconsxminus(j,k,:),Bconsxplus(j,k,:),&
-                    trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
-               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
-                    Bconsy(j,k,:),Bconsyminus(j,k,:),Bconsyplus(j,k,:),&
-                    trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
-               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
-                    Bconsz(j,k,:),Bconszminus(j,k,:),Bconszplus(j,k,:),&
-                    trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
-            endif
-          else
-            !$OMP CRITICAL
-            call CCTK_WARN(0, "Variable type to reconstruct not recognized.")
-            !$OMP END CRITICAL
-          end if
-
-          if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then
-             call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
-                  psidc(j,k,:),psidcminus(j,k,:),psidcplus(j,k,:),&
-                  trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
-          endif
-
-          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)
-  deallocate(dum,dump,dumm)
-
   !$OMP PARALLEL WORKSHARE
   where ( (rhoplus < GRHydro_rho_min).or.(rhominus < GRHydro_rho_min) )
     rhoplus = rho
diff --git a/src/GRHydro_TVDReconstruct_drv.F90 b/src/GRHydro_TVDReconstruct_drv.F90
new file mode 100644
index 0000000..6fc2356
--- /dev/null
+++ b/src/GRHydro_TVDReconstruct_drv.F90
@@ -0,0 +1,238 @@
+ /*@@
+   @file      GRHydro_TVDReconstruct_drv.F90
+   @date      Tue Jul 19 13:22:03 EDT 2011
+   @author    Bruno C. Mundim, Joshua Faber
+   @desc 
+   Driver routine to perform the TVD reconstruction.
+   @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)
+#define Bvecx(i,j,k) Bvec(i,j,k,1)
+#define Bvecy(i,j,k) Bvec(i,j,k,2)
+#define Bvecz(i,j,k) Bvec(i,j,k,3)
+#define Bconsx(i,j,k) Bcons(i,j,k,1)
+#define Bconsy(i,j,k) Bcons(i,j,k,2)
+#define Bconsz(i,j,k) Bcons(i,j,k,3)
+
+
+ /*@@
+   @routine    GRHydro_TVDReconstruct_drv
+   @date       Tue Jul 19 13:24:34 EDT 2011
+   @author     Luca Baiotti, Ian Hawke, Bruno C. Mundim, Joshua Faber
+   @desc 
+   A driver routine to do TVD reconstruction. Currently just does
+   TVD on the primitive variables.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+subroutine GRHydro_TVDReconstruct_drv(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, dum, dump, dumm
+!!$  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
+  
+!!$ The dum variables are used as dummies if MHD on but divergence cleaning off
+  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)),&
+       dum(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),&
+       dump(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),&
+       dumm(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)
+
+  !$OMP PARALLEL
+  !$OMP WORKSHARE
+  trivial_rp = .false.
+  !$OMP END WORKSHARE NOWAIT
+
+!!$  Currently only option is reconstruction on primitive variables.
+!!$  Should change this.
+
+  !$OMP WORKSHARE
+  psi4 = 1.d0
+  !$OMP END WORKSHARE NOWAIT
+  if (shift_state .ne. 0) then
+    !$OMP WORKSHARE
+    lbetax = betax
+    lbetay = betay
+    lbetaz = betaz
+    !$OMP END WORKSHARE NOWAIT
+  else
+    !$OMP WORKSHARE
+    lbetax = 0.d0
+    lbetay = 0.d0
+    lbetaz = 0.d0
+    !$OMP END WORKSHARE NOWAIT
+  end if
+  !$OMP END PARALLEL
+
+
+!!$ TVD starts:
+     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 (evolve_Y_e .ne. 0) then
+        call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
+             Y_e(:,:,:), Y_e_plus(:,:,:), &
+             Y_e_minus(:,:,:), &
+             trivial_rp, hydro_excision_mask)
+     endif
+
+    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)
+      call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
+           eps, epsplus, epsminus, trivial_rp, hydro_excision_mask)
+      if(evolve_mhd.ne.0) then
+         call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
+              Bvec(:,:,:,1), Bvecxplus, Bvecxminus, trivial_rp, hydro_excision_mask)
+         call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
+              Bvec(:,:,:,2), Bvecyplus, Bvecyminus, trivial_rp, hydro_excision_mask)
+         call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
+              Bvec(:,:,:,3), Bveczplus, Bveczminus, trivial_rp, hydro_excision_mask)
+      endif
+
+    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)
+      call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
+           tau, tauplus, tauminus, trivial_rp, hydro_excision_mask)
+      if(evolve_mhd.ne.0) then
+         call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
+              Bcons(:,:,:,1), Bconsxplus, Bconsxminus, trivial_rp, hydro_excision_mask)
+         call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
+              Bcons(:,:,:,2), Bconsyplus, Bconsyminus, trivial_rp, hydro_excision_mask)
+         call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
+              Bcons(:,:,:,3), Bconszplus, Bconszminus, trivial_rp, hydro_excision_mask)
+      endif
+
+    else
+      call CCTK_WARN(0, "Variable type to reconstruct not recognized.")
+    end if
+
+      if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then
+       call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
+            psidc, psidcplus, psidcminus, trivial_rp, hydro_excision_mask)
+    endif
+
+    !$OMP PARALLEL DO PRIVATE(i, j, k)
+    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
+!!$ TVD ends.
+
+  deallocate(trivial_rp)
+  deallocate(psi4, lbetax, lbetay, lbetaz)
+  deallocate(dum,dump,dumm)
+
+end subroutine GRHydro_TVDReconstruct_drv
+
+
diff --git a/src/make.code.defn b/src/make.code.defn
index 27c2ecb..93c6b5d 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -7,9 +7,10 @@ SRCS = 	Utils.F90 \
 	GRHydro_Boundaries.F90 \
 	GRHydro_CalcUpdate.F90 \
 	GRHydro_Con2Prim.F90 \
-        GRHydro_DivergenceClean.F90 \
+   GRHydro_DivergenceClean.F90 \
 	GRHydro_Eigenproblem.F90 \
 	GRHydro_Eigenproblem_Marquina.F90 \
+   GRHydro_ENOReconstruct_drv.F90 \
 	GRHydro_ENOReconstruct.F90 \
 	GRHydro_ENOScalars.F90 \
 	GRHydro_EOS.c \
@@ -19,6 +20,7 @@ SRCS = 	Utils.F90 \
 	GRHydro_Loop.F90 \
 	GRHydro_Minima.F90 \
 	GRHydro_Minima.cc \
+   GRHydro_PPMReconstruct_drv.F90 \
 	GRHydro_PPM.F90 \
 	GRHydro_ParamCheck.F90 \
 	GRHydro_Particle.F90 \
@@ -36,6 +38,7 @@ SRCS = 	Utils.F90 \
 	GRHydro_SlopeLimiter.F90 \
 	GRHydro_Source.F90 \
 	GRHydro_Startup.F90 \
+   GRHydro_TVDReconstruct_drv.F90 \
 	GRHydro_TVDReconstruct.F90 \
 	GRHydro_Marquina.F90 \
 	GRHydro_UpdateMask.F90 \
@@ -49,7 +52,7 @@ SRCS = 	Utils.F90 \
 	GRHydro_EigenproblemM.F90 \
 	GRHydro_FluxM.F90 \
 	GRHydro_HLLEM.F90 \
-        GRHydro_PPMM.F90 \
+   GRHydro_PPMM.F90 \
 	GRHydro_Prim2ConM.F90 \
 	GRHydro_RiemannSolveM.F90 \
 	GRHydro_SourceM.F90 \