RIT GRMHD dev:

1) Further splitting of PPM Reconstruction routine into magnetic
and non-magnetic part.

2) Merge the divergence cleaning loop into the non-divergence one,
since profiling indicated no substantial difference between the
two cases. Indeed it was a little bit better after the merger
(0.14%), but that it is not significant. We decided then to apply 
Occam's razor and choose the simplest form. Branch prediction 
seems to work fine in this case.

3) Move reconstruction initialization statements to their repective 
drivers.

4) Get rid of WORKSHARE in the reconstruction routines. Profiling 
showed a 4.16% performance improvement for the hydro ppm reconstruction 
routine when using 1 processor and 2 threads only. Expect a 
better improvement for a larger number of threads.

5) Introduce a parameter to control characteristic speeds
for psidc in HLLE.



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

11 files changed, 634 insertions, 365 deletions

diff --git a/interface.ccl b/interface.ccl
index 449a888..aef9828 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -453,6 +453,8 @@ real GRHydro_MHD_psidc_bext type = GF Timelevels = 1 tags='Prolongation="None" c
   psidcplus,psidcminus
 } "Divergence cleaning variable extended to the cell boundaries for diverence cleaning"
 
+int whichpsidcspeed type = SCALAR tags='checkpoint="no"' "Which speed to set for psidc? Set in ParamCheck"
+
 # real fluxweightvolume type = GF Timelevels = 1
 # {
 #   cell_volume
diff --git a/param.ccl b/param.ccl
index 399ff7d..a2226f1 100644
--- a/param.ccl
+++ b/param.ccl
@@ -516,6 +516,23 @@ CCTK_REAL cp_dc "The c_p parameter for divergence cleaning"
   0:* :: "Any value, but one to 12 is preferred"
 } 1.0
 
+KEYWORD psidcspeed "Which speed to set for psidc"
+{
+  "char speed"  :: "Based on the characteristic speeds"
+  "light speed" :: "Set the characteristic speeds to speed of light"
+  "set speed"   :: "Manually set the characteristic speeds: [setcharmin,setcharmax]"
+} "light speed"
+
+CCTK_REAL setcharmax "Maximum characteristic speed for psidc if psidcspeed is set"
+{
+  0:1 :: "Any value smaller than speed of light"
+} 1.0
+
+CCTK_REAL setcharmin "Minimum characteristic speed for psidc if psidcspeed is set"
+{
+  -1:0 :: "Any value smaller than speed of light - sign should be negative"
+} -1.0
+
 boolean track_divB "Track the magnetic field constraint violations"
 {
 } "no"
diff --git a/schedule.ccl b/schedule.ccl
index 163818f..2a5a8be 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -76,6 +76,7 @@ if(CCTK_Equals(Bvec_evolution_method,"GRHydro"))
   if (clean_divergence)
   {
     STORAGE:psidcrhs
+    STORAGE:whichpsidcspeed 
   }
   if (track_divB)
   {
diff --git a/src/GRHydro_ENOReconstruct_drv.F90 b/src/GRHydro_ENOReconstruct_drv.F90
index 656fd6e..65ff19c 100644
--- a/src/GRHydro_ENOReconstruct_drv.F90
+++ b/src/GRHydro_ENOReconstruct_drv.F90
@@ -67,8 +67,6 @@ subroutine GRHydro_ENOReconstruct_drv(CCTK_ARGUMENTS)
 
   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")
@@ -107,31 +105,62 @@ subroutine GRHydro_ENOReconstruct_drv(CCTK_ARGUMENTS)
   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
+!!$ Initialize variables that store reconstructed quantities:
+
+  !$OMP PARALLEL DO PRIVATE(i,j,k)
+  do k=1,cctk_lsh(3)
+     do j=1,cctk_lsh(2)
+        do i=1,cctk_lsh(1)
+           trivial_rp(i,j,k) = .false.
+           psi4(i,j,k) = 1.0d0
+           rhoplus(i,j,k) = 0.0d0
+           rhominus(i,j,k)= 0.0d0
+           epsplus(i,j,k) = 0.0d0
+           epsminus(i,j,k) = 0.0d0
+           velxplus(i,j,k) = 0.0d0
+           velxminus(i,j,k) = 0.0d0
+           velyplus(i,j,k) = 0.0d0
+           velyminus(i,j,k) = 0.0d0
+           velzplus(i,j,k) = 0.0d0
+           velzminus(i,j,k) = 0.0d0
+
+           if(evolve_mhd.ne.0) then
+             Bvecxplus(i,j,k) = 0.0d0
+             Bvecxminus(i,j,k) = 0.0d0
+             Bvecyplus(i,j,k) = 0.0d0
+             Bvecyminus(i,j,k) = 0.0d0
+             Bveczplus(i,j,k) = 0.0d0
+             Bveczminus(i,j,k) = 0.0d0
+             if(clean_divergence.ne.0) then
+                psidcplus(i,j,k) = 0.0d0
+                psidcminus(i,j,k) = 0.0d0
+             endif
+           endif
+
+           if (evolve_tracer .ne. 0) then
+              tracerplus(i,j,k,:) = 0.0d0
+              tracerminus(i,j,k,:) = 0.0d0
+           endif
+
+           if (evolve_Y_e .ne. 0) then
+              Y_e_plus(i,j,k) = 0.0d0
+              Y_e_minus(i,j,k) = 0.0d0
+           endif
+
+           if (shift_state .ne. 0) then
+             lbetax(i,j,k) = betax(i,j,k)
+             lbetay(i,j,k) = betay(i,j,k)
+             lbetaz(i,j,k) = betaz(i,j,k)
+           else
+             lbetax(i,j,k) = 0.d0
+             lbetay(i,j,k) = 0.d0
+             lbetaz(i,j,k) = 0.d0
+           end if
+
+        enddo
+     enddo
+  enddo
+  !$OMP END PARALLEL DO
 
 !!$ ENO starts:
      if (evolve_tracer .ne. 0) then
diff --git a/src/GRHydro_HLLEM.F90 b/src/GRHydro_HLLEM.F90
index e841450..b85efb6 100644
--- a/src/GRHydro_HLLEM.F90
+++ b/src/GRHydro_HLLEM.F90
@@ -426,19 +426,24 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS)
           end do
 
           if(clean_divergence.ne.0) then
-             psidcdiff = psidcm - psidcp
-
-!!$             psidcf = (charmax * psidcfp - charmin * psidcfm + &
-!!$                  charmax * charmin * psidcdiff) / charpm
-
-!!$   Wavespeeds for psidc are +/-c, not Fast Magnetosonic?
-
+           psidcdiff = psidcm - psidcp
+           select case(whichpsidcspeed)
+             case(0)
+             psidcf = (charmax * psidcfp - charmin * psidcfm + &
+                  charmax * charmin * psidcdiff) / charpm
+             case(1)
+             !!$ Wavespeeds for psidc are +/-c, not Fast Magnetosonic?
              psidcf = (1.d0 * psidcfp - (-1.d0) * psidcfm + &
                   1.d0 * (-1.d0) * psidcdiff) / 2.d0
-
-
+             case(2)
+             charmax = setcharmax
+             charmin = setcharmin
+             psidcf = (charmax * psidcfp - charmin * psidcfm + &
+                  charmax * charmin * psidcdiff) / charpm
+           end select
           endif
 
+
         end if !!! The end of the SpaceMask check for a trivial RP.
             
         densflux(i, j, k) = f1(1)
diff --git a/src/GRHydro_PPMMReconstruct_drv.F90 b/src/GRHydro_PPMMReconstruct_drv.F90
new file mode 100644
index 0000000..9486a77
--- /dev/null
+++ b/src/GRHydro_PPMMReconstruct_drv.F90
@@ -0,0 +1,386 @@
+ /*@@
+   @file      GRHydro_PPMMReconstruct_drv.F90
+   @date      Wed Jul 27 15:17:03 EDT 2011
+   @author    Bruno C. Mundim, Joshua Faber
+   @desc 
+   Driver routine to perform the magnetic version of PPM 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_PPMMReconstruct_drv
+   @date       Wed Jul 27 15:17:55 EDT 2011
+   @author     Luca Baiotti, Ian Hawke, Bruno C. Mundim, Joshua Faber
+   @desc 
+   A driver routine to do magnetic version of PPM reconstructions. 
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+subroutine GRHydro_PPMMReconstruct_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
+
+  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)
+
+!!$ Initialize variables that store reconstructed quantities:
+
+  !$OMP PARALLEL DO PRIVATE(i,j,k)
+  do k=1,cctk_lsh(3)
+     do j=1,cctk_lsh(2)
+        do i=1,cctk_lsh(1)
+           trivial_rp(i,j,k) = .false.
+           psi4(i,j,k) = 1.0d0
+           rhoplus(i,j,k) = 0.0d0
+           rhominus(i,j,k)= 0.0d0
+           epsplus(i,j,k) = 0.0d0
+           epsminus(i,j,k) = 0.0d0
+           velxplus(i,j,k) = 0.0d0
+           velxminus(i,j,k) = 0.0d0
+           velyplus(i,j,k) = 0.0d0
+           velyminus(i,j,k) = 0.0d0
+           velzplus(i,j,k) = 0.0d0
+           velzminus(i,j,k) = 0.0d0
+
+           Bvecxplus(i,j,k) = 0.0d0
+           Bvecxminus(i,j,k) = 0.0d0
+           Bvecyplus(i,j,k) = 0.0d0
+           Bvecyminus(i,j,k) = 0.0d0
+           Bveczplus(i,j,k) = 0.0d0
+           Bveczminus(i,j,k) = 0.0d0
+           if(clean_divergence.ne.0) then
+             psidcplus(i,j,k) = 0.0d0
+             psidcminus(i,j,k) = 0.0d0
+           endif
+
+           if (evolve_tracer .ne. 0) then
+              tracerplus(i,j,k,:) = 0.0d0
+              tracerminus(i,j,k,:) = 0.0d0
+           endif
+
+           if (evolve_Y_e .ne. 0) then
+              Y_e_plus(i,j,k) = 0.0d0
+              Y_e_minus(i,j,k) = 0.0d0
+           endif
+
+           if (shift_state .ne. 0) then
+             lbetax(i,j,k) = betax(i,j,k)
+             lbetay(i,j,k) = betay(i,j,k)
+             lbetaz(i,j,k) = betaz(i,j,k)
+           else
+             lbetax(i,j,k) = 0.d0
+             lbetay(i,j,k) = 0.d0
+             lbetaz(i,j,k) = 0.d0
+           end if
+
+        enddo
+     enddo
+  enddo
+  !$OMP END PARALLEL DO
+
+!!$ PPM starts:
+    if (flux_direction == 1) then
+            !$OMP PARALLEL DO PRIVATE(i, j, k)
+            do k = GRHydro_stencil, nz - GRHydro_stencil + 1
+              do j = GRHydro_stencil, ny - GRHydro_stencil + 1
+               if(clean_divergence.ne.0) then
+                 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)
+               else  !clean_divergence
+                 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)
+               endif  !clean_divergence
+                do i = 1, nx
+                   if (trivial_rp(i,j,k)) then
+                     SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, trivialx)
+                   else
+                     SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, not_trivialx)
+                   end if
+                end do
+              end do
+            end do
+            !$OMP END PARALLEL DO
+
+      if(evolve_tracer.ne.0) then
+         !$OMP PARALLEL DO PRIVATE(j, 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
+          !$OMP PARALLEL DO PRIVATE(i, j, k)
+          do k = GRHydro_stencil, nz - GRHydro_stencil + 1
+            do j = GRHydro_stencil, nx - GRHydro_stencil + 1
+             if(clean_divergence.ne.0) then
+              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)
+             else  !clean_divergence
+              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)
+             endif  !clean_divergence
+              do i = 1, ny
+                if (trivial_rp(j,i,k)) then
+                   SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, trivialy)
+                else
+                   SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, not_trivialy)
+                end if
+              end do
+            end do
+          end do
+          !$OMP END PARALLEL DO
+      if(evolve_tracer.ne.0) then
+         !$OMP PARALLEL DO PRIVATE(j, 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
+          !$OMP PARALLEL DO PRIVATE(i, j, k)
+          do k = GRHydro_stencil, ny - GRHydro_stencil + 1
+            do j = GRHydro_stencil, nx - GRHydro_stencil + 1
+             if(clean_divergence.ne.0) then
+              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)
+             else  !clean_divergence
+              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)
+             endif  !clean_divergence
+              do i = 1, nz
+                if (trivial_rp(j,k,i)) then
+                  SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, trivialz)
+                else
+                  SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, not_trivialz)
+                end if
+              end do
+            end do
+          end do
+          !$OMP END PARALLEL DO
+      if(evolve_tracer.ne.0) then
+         !$OMP PARALLEL DO PRIVATE(j, 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_PPMMReconstruct_drv
+
diff --git a/src/GRHydro_PPMReconstruct_drv.F90 b/src/GRHydro_PPMReconstruct_drv.F90
index ddd68cd..854d77e 100644
--- a/src/GRHydro_PPMReconstruct_drv.F90
+++ b/src/GRHydro_PPMReconstruct_drv.F90
@@ -20,12 +20,6 @@
 #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)
 
 
  /*@@
@@ -61,30 +55,25 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS)
        &type_bitsz, trivialz, not_trivialz
 
   CCTK_REAL, dimension(:,:,:),allocatable :: &
-       &psi4, lbetax, lbetay, lbetaz, dum, dump, dumm
+       &psi4, lbetax, lbetay, lbetaz
 !!$  CCTK_REAL, dimension(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) :: &
 !!$       &psi4, lbetax, lbetay, lbetaz
 
   CCTK_INT :: ierr
 
-  CCTK_REAL :: local_min_tracer
-
   allocate(trivial_rp(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),STAT=ierr)
   if (ierr .ne. 0) then
     call CCTK_WARN(0, "Allocation problems with trivial_rp")
   end if
   
-!!$ 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)
+       STAT=ierr)
 
   if (ierr .ne. 0) then
-    call CCTK_WARN(0, "Allocation problems with lbeta")
+    call CCTK_WARN(0, "Allocation problems with psi4 or lbeta")
   end if
   
   call SpaceMask_GetTypeBits(type_bitsx, "Hydro_RiemannProblemX")
@@ -107,94 +96,53 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS)
   ny = cctk_lsh(2)
   nz = cctk_lsh(3)
 
-  !$OMP PARALLEL
-  !$OMP WORKSHARE
-  trivial_rp = .false.
-  !$OMP END WORKSHARE NOWAIT
+!!$ Initialize variables that store reconstructed quantities:
 
-!!$  Currently only option is reconstruction on primitive variables.
-!!$  Should change this.
+  !$OMP PARALLEL DO PRIVATE(i,j,k)
+  do k=1,cctk_lsh(3)
+     do j=1,cctk_lsh(2)
+        do i=1,cctk_lsh(1)
+           trivial_rp(i,j,k) = .false.
+           psi4(i,j,k) = 1.0d0
+           rhoplus(i,j,k) = 0.0d0
+           rhominus(i,j,k)= 0.0d0
+           epsplus(i,j,k) = 0.0d0
+           epsminus(i,j,k) = 0.0d0
+           velxplus(i,j,k) = 0.0d0
+           velxminus(i,j,k) = 0.0d0
+           velyplus(i,j,k) = 0.0d0
+           velyminus(i,j,k) = 0.0d0
+           velzplus(i,j,k) = 0.0d0
+           velzminus(i,j,k) = 0.0d0
 
-  !$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
+           if (evolve_tracer .ne. 0) then
+              tracerplus(i,j,k,:) = 0.0d0
+              tracerminus(i,j,k,:) = 0.0d0
+           endif
+
+           if (evolve_Y_e .ne. 0) then
+              Y_e_plus(i,j,k) = 0.0d0
+              Y_e_minus(i,j,k) = 0.0d0
+           endif
+
+           if (shift_state .ne. 0) then
+             lbetax(i,j,k) = betax(i,j,k)
+             lbetay(i,j,k) = betay(i,j,k)
+             lbetaz(i,j,k) = betaz(i,j,k)
+           else
+             lbetax(i,j,k) = 0.d0
+             lbetay(i,j,k) = 0.d0
+             lbetaz(i,j,k) = 0.d0
+           end if
+
+        enddo
+     enddo
+  enddo
+  !$OMP END PARALLEL DO
 
 
 !!$ 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
@@ -220,7 +168,6 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS)
            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
@@ -247,65 +194,6 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS)
       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
@@ -331,7 +219,6 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS)
           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
@@ -358,65 +245,6 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS)
       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
@@ -442,7 +270,6 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS)
           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
@@ -476,7 +303,6 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS)
 
   deallocate(trivial_rp)
   deallocate(psi4, lbetax, lbetay, lbetaz)
-  deallocate(dum,dump,dumm)
 
 end subroutine GRHydro_PPMReconstruct_drv
 
diff --git a/src/GRHydro_ParamCheck.F90 b/src/GRHydro_ParamCheck.F90
index f6177f7..2c29d13 100644
--- a/src/GRHydro_ParamCheck.F90
+++ b/src/GRHydro_ParamCheck.F90
@@ -112,5 +112,18 @@ subroutine GRHydro_ParamCheck(CCTK_ARGUMENTS)
     call CCTK_PARAMWARN("Marquina solver is not implemented yet for MHD")
   end if
 
+  
+  if(clean_divergence.ne.0) then
+    if (CCTK_EQUALS(psidcspeed,"char speed")) then
+      whichpsidcspeed = 0
+    else if (CCTK_EQUALS(psidcspeed,"light speed")) then
+      whichpsidcspeed = 1
+    else if (CCTK_EQUALS(psidcspeed,"set speed")) then
+      whichpsidcspeed = 2
+    else
+      call CCTK_PARAMWARN("Unknown type of speed setting for psidc (psidcspeed)")
+    end if
+  end if
+
 end subroutine GRHydro_ParamCheck
 
diff --git a/src/GRHydro_Reconstruct.F90 b/src/GRHydro_Reconstruct.F90
index a8ac2ec..ca93a4c 100644
--- a/src/GRHydro_Reconstruct.F90
+++ b/src/GRHydro_Reconstruct.F90
@@ -38,112 +38,71 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
   DECLARE_CCTK_PARAMETERS
   DECLARE_CCTK_FUNCTIONS
 
+  CCTK_INT  :: i,j,k
   CCTK_REAL :: local_min_tracer
 
-!!$ Initialize variables that store reconstructed quantities
-
-  !$OMP PARALLEL
-  !$OMP WORKSHARE
-  rhoplus = 0.0d0
-  rhominus = 0.0d0
-  epsplus = 0.0d0
-  epsminus = 0.0d0
-  velxplus = 0.0d0
-  velxminus = 0.0d0
-  velyplus = 0.0d0
-  velyminus = 0.0d0
-  velzplus = 0.0d0
-  velzminus = 0.0d0
-  !$OMP END WORKSHARE NOWAIT
-
-  if(evolve_mhd.ne.0) then
-     !$OMP WORKSHARE
-     Bvecxplus = 0.0d0
-     Bvecxminus = 0.0d0
-     Bvecyplus = 0.0d0
-     Bvecyminus = 0.0d0
-     Bveczplus = 0.0d0
-     Bveczminus = 0.0d0
-     !$OMP END WORKSHARE NOWAIT
-     if(clean_divergence.ne.0) then
-        !$OMP WORKSHARE
-        psidcplus = 0.0d0
-        psidcminus = 0.0d0
-        !$OMP END WORKSHARE NOWAIT
-     endif
-  endif
-
-  if (evolve_tracer .ne. 0) then
-     !$OMP WORKSHARE
-     tracerplus = 0.0d0
-     tracerminus = 0.0d0
-     !$OMP END WORKSHARE NOWAIT
-  endif
-
-  if (evolve_Y_e .ne. 0) then
-     !$OMP WORKSHARE
-     Y_e_plus = 0.0d0
-     Y_e_minus = 0.0d0
-     !$OMP END WORKSHARE NOWAIT
-  endif
-  !$OMP END PARALLEL
-
 
   if (CCTK_EQUALS(recon_method,"tvd")) then
-
+    ! this handles MHD and non-MHD
     call GRHydro_TVDReconstruct_drv(CCTK_PASS_FTOF)
 
   else if (CCTK_EQUALS(recon_method,"ppm")) then
 
-    call GRHydro_PPMReconstruct_drv(CCTK_PASS_FTOF)
+    if(evolve_mhd.ne.0) then
+      call GRHydro_PPMMReconstruct_drv(CCTK_PASS_FTOF)
+    else
+      call GRHydro_PPMReconstruct_drv(CCTK_PASS_FTOF)
+    end if 
 
   else if (CCTK_EQUALS(recon_method,"eno")) then
-
+    ! this handles MHD and non-MHD
     call GRHydro_ENOReconstruct_drv(CCTK_PASS_FTOF)
 
   else
     call CCTK_WARN(0, "Reconstruction method not recognized!")
   end if
 
-  !$OMP PARALLEL WORKSHARE
-  where ( (rhoplus < GRHydro_rho_min).or.(rhominus < GRHydro_rho_min) )
-    rhoplus = rho
-    rhominus = rho
-    velxplus = vel(:,:,:,1)
-    velxminus = vel(:,:,:,1)
-    velyplus = vel(:,:,:,2)
-    velyminus = vel(:,:,:,2)
-    velzplus = vel(:,:,:,3)
-    velzminus = vel(:,:,:,3)
-    epsplus = eps
-    epsminus = eps
-  end where
-  !$OMP END PARALLEL WORKSHARE
-
   if (evolve_tracer .ne. 0) then
     if (use_min_tracer .ne. 0) then
       local_min_tracer = min_tracer
     else
       local_min_tracer = 0.0d0
     end if
+  end if
 
-    !$OMP PARALLEL WORKSHARE
-    where( (tracerplus  .le. local_min_tracer).or.&
-           (tracerminus .le. local_min_tracer) )
-      tracerplus = tracer
-      tracerminus = tracer
-    end where
-    !$OMP END PARALLEL WORKSHARE
-    ! Call the conserved tracer routine in any case because (accord. to
-    ! Christian Ott) this is the only way this works
-    call Prim2ConservativeTracer(CCTK_PASS_FTOF)
-  endif
+  !$OMP PARALLEL DO PRIVATE(i,j,k)
+     do k=1,cctk_lsh(3)
+        do j=1,cctk_lsh(2)
+           do i=1,cctk_lsh(1)
+              if(rhoplus(i,j,k).lt.GRHydro_rho_min .or. &
+                   rhominus(i,j,k).lt.GRHydro_rho_min) then
+                 rhoplus(i,j,k)  = rho(i,j,k)
+                 rhominus(i,j,k) = rho(i,j,k)
+                 epsplus(i,j,k)  = eps(i,j,k)
+                 epsminus(i,j,k) = eps(i,j,k)
+                 velxplus(i,j,k) = vel(i,j,k,1)
+                 velyplus(i,j,k) = vel(i,j,k,2)
+                 velzplus(i,j,k) = vel(i,j,k,3)
+                 velxminus(i,j,k) = vel(i,j,k,1)
+                 velyminus(i,j,k) = vel(i,j,k,2)
+                 velzminus(i,j,k) = vel(i,j,k,3)
+              end if
+              if(evolve_tracer.ne.0) then
+                 where(tracerplus(i,j,k,:).le.local_min_tracer .or. &
+                      tracerminus(i,j,k,:).le.local_min_tracer)
+                   tracerplus(i,j,k,:) = tracer(i,j,k,:)
+                   tracerminus(i,j,k,:) = tracer(i,j,k,:)
+                end where
+             end if
+          enddo
+       enddo
+    enddo
+    !$OMP END PARALLEL DO
 
   if (CCTK_EQUALS(recon_vars,"primitive").or.&
        CCTK_EQUALS(recon_method,"ppm")) then
      if(evolve_mhd.ne.0) then
         call primitive2conservativeM(CCTK_PASS_FTOF)
-
      else
         call primitive2conservative(CCTK_PASS_FTOF)
      endif
@@ -157,6 +116,12 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
     call CCTK_WARN(0,"Variable type to reconstruct not recognized.")
   end if
 
+ if(evolve_tracer.ne.0) then
+    ! Call the conserved tracer routine in any case because (accord. to
+    ! Christian Ott) this is the only way this works
+    call Prim2ConservativeTracer(CCTK_PASS_FTOF)
+ end if
+
   return
 
 end subroutine Reconstruction
diff --git a/src/GRHydro_TVDReconstruct_drv.F90 b/src/GRHydro_TVDReconstruct_drv.F90
index 6fc2356..f5bf492 100644
--- a/src/GRHydro_TVDReconstruct_drv.F90
+++ b/src/GRHydro_TVDReconstruct_drv.F90
@@ -62,27 +62,22 @@ subroutine GRHydro_TVDReconstruct_drv(CCTK_ARGUMENTS)
        &type_bitsz, trivialz, not_trivialz
 
   CCTK_REAL, dimension(:,:,:),allocatable :: &
-       &psi4, lbetax, lbetay, lbetaz, dum, dump, dumm
+       &psi4, lbetax, lbetay, lbetaz
 !!$  CCTK_REAL, dimension(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) :: &
 !!$       &psi4, lbetax, lbetay, lbetaz
 
   CCTK_INT :: ierr
 
-  CCTK_REAL :: local_min_tracer
-
   allocate(trivial_rp(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),STAT=ierr)
   if (ierr .ne. 0) then
     call CCTK_WARN(0, "Allocation problems with trivial_rp")
   end if
   
-!!$ 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)
+       STAT=ierr)
 
   if (ierr .ne. 0) then
     call CCTK_WARN(0, "Allocation problems with lbeta")
@@ -108,32 +103,62 @@ subroutine GRHydro_TVDReconstruct_drv(CCTK_ARGUMENTS)
   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
+!!$ Initialize variables that store reconstructed quantities:
+
+  !$OMP PARALLEL DO PRIVATE(i,j,k)
+  do k=1,cctk_lsh(3)
+     do j=1,cctk_lsh(2)
+        do i=1,cctk_lsh(1)
+           trivial_rp(i,j,k) = .false.
+           psi4(i,j,k) = 1.0d0
+           rhoplus(i,j,k) = 0.0d0
+           rhominus(i,j,k)= 0.0d0
+           epsplus(i,j,k) = 0.0d0
+           epsminus(i,j,k) = 0.0d0
+           velxplus(i,j,k) = 0.0d0
+           velxminus(i,j,k) = 0.0d0
+           velyplus(i,j,k) = 0.0d0
+           velyminus(i,j,k) = 0.0d0
+           velzplus(i,j,k) = 0.0d0
+           velzminus(i,j,k) = 0.0d0
+
+           if(evolve_mhd.ne.0) then
+             Bvecxplus(i,j,k) = 0.0d0
+             Bvecxminus(i,j,k) = 0.0d0
+             Bvecyplus(i,j,k) = 0.0d0
+             Bvecyminus(i,j,k) = 0.0d0
+             Bveczplus(i,j,k) = 0.0d0
+             Bveczminus(i,j,k) = 0.0d0
+             if(clean_divergence.ne.0) then
+                psidcplus(i,j,k) = 0.0d0
+                psidcminus(i,j,k) = 0.0d0
+             endif
+           endif
+
+           if (evolve_tracer .ne. 0) then
+              tracerplus(i,j,k,:) = 0.0d0
+              tracerminus(i,j,k,:) = 0.0d0
+           endif
+
+           if (evolve_Y_e .ne. 0) then
+              Y_e_plus(i,j,k) = 0.0d0
+              Y_e_minus(i,j,k) = 0.0d0
+           endif
+
+           if (shift_state .ne. 0) then
+             lbetax(i,j,k) = betax(i,j,k)
+             lbetay(i,j,k) = betay(i,j,k)
+             lbetaz(i,j,k) = betaz(i,j,k)
+           else
+             lbetax(i,j,k) = 0.d0
+             lbetay(i,j,k) = 0.d0
+             lbetaz(i,j,k) = 0.d0
+           end if
 
+        enddo
+     enddo
+  enddo
+  !$OMP END PARALLEL DO
 
 !!$ TVD starts:
      if (evolve_tracer .ne. 0) then
@@ -231,7 +256,6 @@ subroutine GRHydro_TVDReconstruct_drv(CCTK_ARGUMENTS)
 
   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 93c6b5d..a4e0873 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -52,6 +52,7 @@ SRCS = 	Utils.F90 \
 	GRHydro_EigenproblemM.F90 \
 	GRHydro_FluxM.F90 \
 	GRHydro_HLLEM.F90 \
+   GRHydro_PPMMReconstruct_drv.F90 \
    GRHydro_PPMM.F90 \
 	GRHydro_Prim2ConM.F90 \
 	GRHydro_RiemannSolveM.F90 \