* Optimize: remove support for shift_state = 0 (except for shock tubes and

  Cowling calculations of spherically symmetric objects, there is no reason
  not to have storage for the shift).




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

14 files changed, 162 insertions, 302 deletions

diff --git a/src/GRHydro_FluxSplit.F90 b/src/GRHydro_FluxSplit.F90
index d54e446..1c701e8 100644
--- a/src/GRHydro_FluxSplit.F90
+++ b/src/GRHydro_FluxSplit.F90
@@ -73,11 +73,7 @@ subroutine GRHydro_FSAlpha(CCTK_ARGUMENTS)
                 gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),&
                 gyz(i,j,k),gzz(i,j,k))        
 
-          if (shift_state .eq. 0) then
-            beta = 0.d0
-          else
-            beta = betax(i,j,k)
-          end if
+           beta = betax(i,j,k)
           
           call eigenvalues(GRHydro_eos_handle, &
                rho      (i,j,k), &
@@ -119,11 +115,7 @@ subroutine GRHydro_FSAlpha(CCTK_ARGUMENTS)
                 gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),&
                 gyz(i,j,k),gzz(i,j,k))        
 
-          if (shift_state .eq. 0) then
-            beta = 0.d0
-          else
-            beta = betay(i,j,k)
-          end if
+           beta = betay(i,j,k)
           
           call eigenvalues(GRHydro_eos_handle, &
                rho      (i,j,k), &
@@ -165,11 +157,7 @@ subroutine GRHydro_FSAlpha(CCTK_ARGUMENTS)
                 gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),&
                 gyz(i,j,k),gzz(i,j,k))        
 
-          if (shift_state .eq. 0) then
-            beta = 0.d0
-          else
-            beta = betaz(i,j,k)
-          end if
+           beta = betaz(i,j,k)
           
           call eigenvalues(GRHydro_eos_handle, &
                rho      (i,j,k), &
@@ -310,11 +298,8 @@ subroutine GRHydro_SplitFlux(CCTK_ARGUMENTS)
     do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + 2
       do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + 2
 
-        if (shift_state .eq. 0) then
-          dummy = 0.d0
-        else
-          dummy = betax(:,j,k)
-        end if
+         dummy = betax(:,j,k)
+
         do i = 1, cctk_lsh(1)
            det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),\
                 gyy(i,j,k),gyz(i,j,k),gzz(i,j,k))
@@ -347,11 +332,8 @@ subroutine GRHydro_SplitFlux(CCTK_ARGUMENTS)
     do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + 2
       do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + 2
 
-        if (shift_state .eq. 0) then
-          dummy = 0.d0
-        else
-          dummy = betay(i,:,k)
-        end if
+         dummy = betay(i,:,k)
+
         do j = 1, cctk_lsh(2)
            det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),\
                 gyy(i,j,k),gyz(i,j,k),gzz(i,j,k))
@@ -384,12 +366,8 @@ subroutine GRHydro_SplitFlux(CCTK_ARGUMENTS)
     do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + 2
       do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + 2
 
-        if (shift_state .eq. 0) then
-          dummy = 0.d0
-        else
-          dummy = betaz(i,j,:)
-        end if
-
+         dummy = betaz(i,j,:)
+         
         do k = 1, cctk_lsh(3)
            det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),\
                 gyy(i,j,k),gyz(i,j,k),gzz(i,j,k))
diff --git a/src/GRHydro_HLLE.F90 b/src/GRHydro_HLLE.F90
index 467959a..2074b60 100644
--- a/src/GRHydro_HLLE.F90
+++ b/src/GRHydro_HLLE.F90
@@ -103,22 +103,18 @@ subroutine GRHydro_HLLE(CCTK_ARGUMENTS)
 !!$        Note also need the average of the lapse at the 
 !!$        left and right points.
 
-        if (shift_state .ne. 0) then
-          if (flux_direction == 1) then
-            avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
-                 betax(i,j,k))
-          else if (flux_direction == 2) then
-            avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
-                 betay(i,j,k))
-          else if (flux_direction == 3) then
-            avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
-                 betaz(i,j,k))
-          else
-            call CCTK_WARN(0, "Flux direction not x,y,z")
-          end if
+        if (flux_direction == 1) then
+           avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
+                betax(i,j,k))
+        else if (flux_direction == 2) then
+           avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
+                betay(i,j,k))
+        else if (flux_direction == 3) then
+           avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
+                betaz(i,j,k))
         else
-          avg_beta = 0.d0
-        end if
+            call CCTK_WARN(0, "Flux direction not x,y,z")
+         end if
 
         avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset))
         
@@ -470,21 +466,17 @@ subroutine GRHydro_HLLE_Tracer(CCTK_ARGUMENTS)
 !!$        Note also need the average of the lapse at the 
 !!$        left and right points.
 
-        if (shift_state .ne. 0) then
-          if (flux_direction == 1) then
-            avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
-                 betax(i,j,k))
-          else if (flux_direction == 2) then
-            avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
-                 betay(i,j,k))
-          else if (flux_direction == 3) then
-            avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
-                 betaz(i,j,k))
-          else
-            call CCTK_WARN(0, "Flux direction not x,y,z")
-          end if
+        if (flux_direction == 1) then
+           avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
+                betax(i,j,k))
+        else if (flux_direction == 2) then
+           avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
+                betay(i,j,k))
+        else if (flux_direction == 3) then
+           avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
+                betaz(i,j,k))
         else
-          avg_beta = 0.d0
+           call CCTK_WARN(0, "Flux direction not x,y,z")
         end if
 
         avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset))
diff --git a/src/GRHydro_HLLEM.F90 b/src/GRHydro_HLLEM.F90
index b85efb6..200575f 100644
--- a/src/GRHydro_HLLEM.F90
+++ b/src/GRHydro_HLLEM.F90
@@ -144,27 +144,20 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS)
 !!$        
 !!$        In MHD, we need all three shift components regardless of the flux direction
 
-        if (shift_state .ne. 0) then
-           avg_betax = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
-                betax(i,j,k))
-           avg_betay = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
-                betay(i,j,k))
-           avg_betaz = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
-                betaz(i,j,k))
-           if (flux_direction == 1) then
-              avg_beta=avg_betax
-           else if (flux_direction == 2) then
-              avg_beta=avg_betay
-           else if (flux_direction == 3) then
-              avg_beta=avg_betaz
-           else
-              call CCTK_WARN(0, "Flux direction not x,y,z")
-           end if
+        avg_betax = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
+             betax(i,j,k))
+        avg_betay = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
+             betay(i,j,k))
+        avg_betaz = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
+             betaz(i,j,k))
+        if (flux_direction == 1) then
+           avg_beta=avg_betax
+        else if (flux_direction == 2) then
+           avg_beta=avg_betay
+        else if (flux_direction == 3) then
+           avg_beta=avg_betaz
         else
-           avg_beta = 0.d0
-           avg_betax = 0.d0
-           avg_betay = 0.d0
-           avg_betaz = 0.d0
+           call CCTK_WARN(0, "Flux direction not x,y,z")
         end if
 
         avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset))
@@ -580,21 +573,17 @@ subroutine GRHydro_HLLE_TracerM(CCTK_ARGUMENTS)
 !!$        Note also need the average of the lapse at the 
 !!$        left and right points.
 
-        if (shift_state .ne. 0) then
-          if (flux_direction == 1) then
-            avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
-                 betax(i,j,k))
-          else if (flux_direction == 2) then
-            avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
-                 betay(i,j,k))
-          else if (flux_direction == 3) then
-            avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
-                 betaz(i,j,k))
-          else
-            call CCTK_WARN(0, "Flux direction not x,y,z")
-          end if
+        if (flux_direction == 1) then
+           avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
+                betax(i,j,k))
+        else if (flux_direction == 2) then
+           avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
+                betay(i,j,k))
+        else if (flux_direction == 3) then
+           avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
+                betaz(i,j,k))
         else
-          avg_beta = 0.d0
+           call CCTK_WARN(0, "Flux direction not x,y,z")
         end if
 
         avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset))
diff --git a/src/GRHydro_Marquina.F90 b/src/GRHydro_Marquina.F90
index 58dd742..ac54b0a 100644
--- a/src/GRHydro_Marquina.F90
+++ b/src/GRHydro_Marquina.F90
@@ -111,19 +111,15 @@ subroutine GRHydro_Marquina(CCTK_ARGUMENTS)
         
 !!$        Set metric terms at interface
           
-          if (shift_state .ne. 0) then
-            if (flux_direction == 1) then
-              avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
-                   betax(i,j,k))
-            else if (flux_direction == 2) then
-              avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
-                   betay(i,j,k))
-            else
-              avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
-                   betaz(i,j,k))
-            end if
+          if (flux_direction == 1) then
+             avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
+                  betax(i,j,k))
+          else if (flux_direction == 2) then
+             avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
+                  betay(i,j,k))
           else
-            avg_beta = 0.d0
+             avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
+                  betaz(i,j,k))
           end if
 
           avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset))
diff --git a/src/GRHydro_ParamCheck.F90 b/src/GRHydro_ParamCheck.F90
index 2c29d13..ba6e641 100644
--- a/src/GRHydro_ParamCheck.F90
+++ b/src/GRHydro_ParamCheck.F90
@@ -92,6 +92,10 @@ subroutine GRHydro_ParamCheck(CCTK_ARGUMENTS)
      evolve_MHD = 0
   endif
 
+  if(shift_state.eq.0) then
+    call CCTK_PARAMWARN("shift_state = 0 (no shift storage) no longer supported!d");
+  endif
+
   if (CCTK_EQUALS(Y_e_evolution_method,"GRHydro")) then
      evolve_Y_e = 1
   else
diff --git a/src/GRHydro_Reconstruct.F90 b/src/GRHydro_Reconstruct.F90
index 09cbf16..cbcfd48 100644
--- a/src/GRHydro_Reconstruct.F90
+++ b/src/GRHydro_Reconstruct.F90
@@ -41,11 +41,6 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
   CCTK_INT  :: i,j,k
   CCTK_REAL :: local_min_tracer
 
-  ! set things up
-  if (shift_state .eq. 0) then
-     call CCTK_WARN(0,"This code no longer supports shift_state = 0");
-  endif
-
   if (CCTK_EQUALS(recon_method,"tvd")) then
     ! this handles MHD and non-MHD
     call GRHydro_TVDReconstruct_drv(CCTK_PASS_FTOF)
diff --git a/src/GRHydro_ReconstructM.F90 b/src/GRHydro_ReconstructM.F90
index a48d67f..9eb44a1 100644
--- a/src/GRHydro_ReconstructM.F90
+++ b/src/GRHydro_ReconstructM.F90
@@ -107,15 +107,9 @@ subroutine ReconstructionM(CCTK_ARGUMENTS)
 
   psi4 = 1.d0
 
-  if (shift_state .ne. 0) then
-    lbetax = betax
-    lbetay = betay
-    lbetaz = betaz
-  else
-    lbetax = 0.d0
-    lbetay = 0.d0
-    lbetaz = 0.d0
-  end if
+  lbetax = betax
+  lbetay = betay
+  lbetaz = betaz
 
 !!$ Initialize variables that store reconstructed quantities
 
diff --git a/src/GRHydro_ReconstructPoly.F90 b/src/GRHydro_ReconstructPoly.F90
index 9d1f860..c4ed3af 100644
--- a/src/GRHydro_ReconstructPoly.F90
+++ b/src/GRHydro_ReconstructPoly.F90
@@ -118,19 +118,12 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS)
   !$OMP WORKSHARE
   psi4 = 1.0d0
   !$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 WORKSHARE
+  lbetax = betax
+  lbetay = betay
+  lbetaz = betaz
+  !$OMP END WORKSHARE NOWAIT
 
 !!$ Initialize variables that store reconstructed quantities
 
diff --git a/src/GRHydro_ReconstructPolyM.F90 b/src/GRHydro_ReconstructPolyM.F90
index 9b4bcd6..6ab9bfd 100644
--- a/src/GRHydro_ReconstructPolyM.F90
+++ b/src/GRHydro_ReconstructPolyM.F90
@@ -107,15 +107,9 @@ subroutine ReconstructionPolytypeM(CCTK_ARGUMENTS)
 
   psi4 = 1.d0
 
-  if (shift_state .ne. 0) then
-    lbetax = betax
-    lbetay = betay
-    lbetaz = betaz
-  else
-    lbetax = 0.d0
-    lbetay = 0.d0
-    lbetaz = 0.d0
-  end if
+  lbetax = betax
+  lbetay = betay
+  lbetaz = betaz
 
 !!$ Initialize variables that store reconstructed quantities
 
diff --git a/src/GRHydro_RoeSolver.F90 b/src/GRHydro_RoeSolver.F90
index 26f9999..18d2f4d 100644
--- a/src/GRHydro_RoeSolver.F90
+++ b/src/GRHydro_RoeSolver.F90
@@ -107,21 +107,17 @@ subroutine GRHydro_RoeSolve(CCTK_ARGUMENTS)
 
 !!$        Set metric terms at interface
         
-        if (shift_state .ne. 0) then
-          if (flux_direction == 1) then
-            avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
-                 betax(i,j,k))
-          else if (flux_direction == 2) then
-            avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
+        if (flux_direction == 1) then
+           avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
+                betax(i,j,k))
+        else if (flux_direction == 2) then
+           avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
                  betay(i,j,k))
-          else if (flux_direction == 3) then
-            avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
-                 betaz(i,j,k))
-          else
-            call CCTK_WARN(0, "Flux direction not x,y,z")
-          end if
+        else if (flux_direction == 3) then
+           avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
+                betaz(i,j,k))
         else
-          avg_beta = 0.d0
+           call CCTK_WARN(0, "Flux direction not x,y,z")
         end if
         
         avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset))
diff --git a/src/GRHydro_Source.F90 b/src/GRHydro_Source.F90
index 1a981db..f8fd9eb 100644
--- a/src/GRHydro_Source.F90
+++ b/src/GRHydro_Source.F90
@@ -175,62 +175,40 @@ subroutine SourceTerms(CCTK_ARGUMENTS)
         rhoenthalpyW2 = (rho(i,j,k)*(one + eps(i,j,k)) + press(i,j,k))*&
              w_lorentz(i,j,k)**2
         
-        if (shift_state .ne. 0) then 
+        shiftx = betax(i,j,k)
+        shifty = betay(i,j,k)
+        shiftz = betaz(i,j,k)
 
-          shiftx = betax(i,j,k)
-          shifty = betay(i,j,k)
-          shiftz = betaz(i,j,k)
-
-          if (local_spatial_order .eq. 2) then
-
-            dx_betax = DIFF_X_2(betax)
-            dx_betay = DIFF_X_2(betay)
-            dx_betaz = DIFF_X_2(betaz)
-            
-            dy_betax = DIFF_Y_2(betax)
-            dy_betay = DIFF_Y_2(betay)
-            dy_betaz = DIFF_Y_2(betaz)
-            
-            dz_betax = DIFF_Z_2(betax)
-            dz_betay = DIFF_Z_2(betay)
-            dz_betaz = DIFF_Z_2(betaz)
-
-          else
-
-            dx_betax = DIFF_X_4(betax)
-            dx_betay = DIFF_X_4(betay)
-            dx_betaz = DIFF_X_4(betaz)
-            
-            dy_betax = DIFF_Y_4(betax)
-            dy_betay = DIFF_Y_4(betay)
-            dy_betaz = DIFF_Y_4(betaz)
-            
-            dz_betax = DIFF_Z_4(betax)
-            dz_betay = DIFF_Z_4(betay)
-            dz_betaz = DIFF_Z_4(betaz)
-
-          end if
-          
+        if (local_spatial_order .eq. 2) then
+           
+           dx_betax = DIFF_X_2(betax)
+           dx_betay = DIFF_X_2(betay)
+           dx_betaz = DIFF_X_2(betaz)
+           
+           dy_betax = DIFF_Y_2(betax)
+           dy_betay = DIFF_Y_2(betay)
+           dy_betaz = DIFF_Y_2(betaz)
+           
+           dz_betax = DIFF_Z_2(betax)
+           dz_betay = DIFF_Z_2(betay)
+           dz_betaz = DIFF_Z_2(betaz)
+           
         else
 
-          shiftx = 0.0d0
-          shifty = 0.0d0
-          shiftz = 0.0d0
-
-          dx_betax = 0.0d0
-          dx_betay = 0.0d0
-          dx_betaz = 0.0d0
-          
-          dy_betax = 0.0d0
-          dy_betay = 0.0d0
-          dy_betaz = 0.0d0
-          
-          dz_betax = 0.0d0
-          dz_betay = 0.0d0
-          dz_betaz = 0.0d0
+           dx_betax = DIFF_X_4(betax)
+           dx_betay = DIFF_X_4(betay)
+           dx_betaz = DIFF_X_4(betaz)
+           
+           dy_betax = DIFF_Y_4(betax)
+           dy_betay = DIFF_Y_4(betay)
+           dy_betaz = DIFF_Y_4(betaz)
+            
+           dz_betax = DIFF_Z_4(betax)
+           dz_betay = DIFF_Z_4(betay)
+           dz_betaz = DIFF_Z_4(betaz)
+           
+        end if
           
-        endif
-
         invalp = 1.0d0 / alp(i,j,k)
         invalp2 = invalp**2
         velxshift = velx(i,j,k) - shiftx*invalp
diff --git a/src/GRHydro_SourceM.F90 b/src/GRHydro_SourceM.F90
index ba2f383..5245df9 100644
--- a/src/GRHydro_SourceM.F90
+++ b/src/GRHydro_SourceM.F90
@@ -199,62 +199,40 @@ subroutine SourceTermsM(CCTK_ARGUMENTS)
              localgxy, localgxz, localgyy, localgyz, localgzz)
         
 
-        if (shift_state .ne. 0) then 
-
-          shiftx = betax(i,j,k)
-          shifty = betay(i,j,k)
-          shiftz = betaz(i,j,k)
-
-          if (local_spatial_order .eq. 2) then
-
-            dx_betax = DIFF_X_2(betax)
-            dx_betay = DIFF_X_2(betay)
-            dx_betaz = DIFF_X_2(betaz)
-            
-            dy_betax = DIFF_Y_2(betax)
-            dy_betay = DIFF_Y_2(betay)
-            dy_betaz = DIFF_Y_2(betaz)
-            
-            dz_betax = DIFF_Z_2(betax)
-            dz_betay = DIFF_Z_2(betay)
-            dz_betaz = DIFF_Z_2(betaz)
-
-          else
-
-            dx_betax = DIFF_X_4(betax)
-            dx_betay = DIFF_X_4(betay)
-            dx_betaz = DIFF_X_4(betaz)
-            
-            dy_betax = DIFF_Y_4(betax)
-            dy_betay = DIFF_Y_4(betay)
-            dy_betaz = DIFF_Y_4(betaz)
+        shiftx = betax(i,j,k)
+        shifty = betay(i,j,k)
+        shiftz = betaz(i,j,k)
+        
+        if (local_spatial_order .eq. 2) then
+           
+           dx_betax = DIFF_X_2(betax)
+           dx_betay = DIFF_X_2(betay)
+           dx_betaz = DIFF_X_2(betaz)
+           
+           dy_betax = DIFF_Y_2(betax)
+           dy_betay = DIFF_Y_2(betay)
+           dy_betaz = DIFF_Y_2(betaz)
             
-            dz_betax = DIFF_Z_4(betax)
-            dz_betay = DIFF_Z_4(betay)
-            dz_betaz = DIFF_Z_4(betaz)
+           dz_betax = DIFF_Z_2(betax)
+           dz_betay = DIFF_Z_2(betay)
+           dz_betaz = DIFF_Z_2(betaz)
 
-          end if
-          
         else
 
-          shiftx = 0.0d0
-          shifty = 0.0d0
-          shiftz = 0.0d0
-
-          dx_betax = 0.0d0
-          dx_betay = 0.0d0
-          dx_betaz = 0.0d0
-          
-          dy_betax = 0.0d0
-          dy_betay = 0.0d0
-          dy_betaz = 0.0d0
-          
-          dz_betax = 0.0d0
-          dz_betay = 0.0d0
-          dz_betaz = 0.0d0
+           dx_betax = DIFF_X_4(betax)
+           dx_betay = DIFF_X_4(betay)
+           dx_betaz = DIFF_X_4(betaz)
+           
+           dy_betax = DIFF_Y_4(betax)
+           dy_betay = DIFF_Y_4(betay)
+           dy_betaz = DIFF_Y_4(betaz)
+           
+           dz_betax = DIFF_Z_4(betax)
+           dz_betay = DIFF_Z_4(betay)
+           dz_betaz = DIFF_Z_4(betaz)
+           
+        end if
           
-        endif
-
         invalp = 1.0d0 / alp(i,j,k)
         invalp2 = invalp**2
         velxshift = velx(i,j,k) - shiftx*invalp
diff --git a/src/GRHydro_Tmunu.F90 b/src/GRHydro_Tmunu.F90
index 6109fb7..1892f53 100644
--- a/src/GRHydro_Tmunu.F90
+++ b/src/GRHydro_Tmunu.F90
@@ -92,24 +92,12 @@
 
 !!$       Calculate lower components and square of shift vector.
 
-         if (shift_state .ne. 0) then
-
-            betaxlow = gxx(i,j,k)*betax(i,j,k) + gxy(i,j,k)*betay(i,j,k) + gxz(i,j,k)*betaz(i,j,k)
-            betaylow = gxy(i,j,k)*betax(i,j,k) + gyy(i,j,k)*betay(i,j,k) + gyz(i,j,k)*betaz(i,j,k)
-            betazlow = gxz(i,j,k)*betax(i,j,k) + gyz(i,j,k)*betay(i,j,k) + gzz(i,j,k)*betaz(i,j,k)
-
-            beta2 = betax(i,j,k)*betaxlow + betay(i,j,k)*betaylow + betaz(i,j,k)*betazlow 
+          betaxlow = gxx(i,j,k)*betax(i,j,k) + gxy(i,j,k)*betay(i,j,k) + gxz(i,j,k)*betaz(i,j,k)
+          betaylow = gxy(i,j,k)*betax(i,j,k) + gyy(i,j,k)*betay(i,j,k) + gyz(i,j,k)*betaz(i,j,k)
+          betazlow = gxz(i,j,k)*betax(i,j,k) + gyz(i,j,k)*betay(i,j,k) + gzz(i,j,k)*betaz(i,j,k)
+          
+          beta2 = betax(i,j,k)*betaxlow + betay(i,j,k)*betaylow + betaz(i,j,k)*betazlow 
            
-         else
-
-           betaxlow = 0.0D0
-           betaylow = 0.0D0
-           betazlow = 0.0D0
-
-           beta2 = 0.0D0
-
-         end if
-
 !!$       Calculate the specific relativistic enthalpy times rho times the
 !!$       square of the lorentz factor.
 
diff --git a/src/GRHydro_TmunuM.F90 b/src/GRHydro_TmunuM.F90
index dfb7ab5..e4df1e0 100644
--- a/src/GRHydro_TmunuM.F90
+++ b/src/GRHydro_TmunuM.F90
@@ -100,35 +100,20 @@
           
 !!$       Calculate lower components and square of shift vector.
           
-          if (shift_state .ne. 0) then
-             
-             betaxlow = gxx(i,j,k)*betax(i,j,k) + gxy(i,j,k)*betay(i,j,k) + gxz(i,j,k)*betaz(i,j,k)
-             betaylow = gxy(i,j,k)*betax(i,j,k) + gyy(i,j,k)*betay(i,j,k) + gyz(i,j,k)*betaz(i,j,k)
-             betazlow = gxz(i,j,k)*betax(i,j,k) + gyz(i,j,k)*betay(i,j,k) + gzz(i,j,k)*betaz(i,j,k)
-             beta2 = betax(i,j,k)*betaxlow + betay(i,j,k)*betaylow + betaz(i,j,k)*betazlow 
-            
-             bdotbeta = betaxlow*Bvecx(i,j,k)+betaylow*Bvecy(i,j,k)+betazlow*Bvecz(i,j,k)
-             vdotbeta = betaxlow*velx(i,j,k)+betaylow*vely(i,j,k)+betazlow*velz(i,j,k)
-
-!!$   u0 low is missing the w_lorentz factor (see below)!!
-             utlow = -1.d0*alp(i,j,k) + vdotbeta
-
-             btlow = -1.0d0*w_lorentz(i,j,k)*alp(i,j,k)*bdotv + &
-                  bdotbeta/w_lorentz(i,j,k) + w_lorentz(i,j,k)*bdotv*vdotbeta
-
-            
-         else
-            
-            betaxlow = 0.0D0
-            betaylow = 0.0D0
-            betazlow = 0.0D0            
-            beta2 = 0.0D0
+          
+          betaxlow = gxx(i,j,k)*betax(i,j,k) + gxy(i,j,k)*betay(i,j,k) + gxz(i,j,k)*betaz(i,j,k)
+          betaylow = gxy(i,j,k)*betax(i,j,k) + gyy(i,j,k)*betay(i,j,k) + gyz(i,j,k)*betaz(i,j,k)
+          betazlow = gxz(i,j,k)*betax(i,j,k) + gyz(i,j,k)*betay(i,j,k) + gzz(i,j,k)*betaz(i,j,k)
+          beta2 = betax(i,j,k)*betaxlow + betay(i,j,k)*betaylow + betaz(i,j,k)*betazlow 
+          
+          bdotbeta = betaxlow*Bvecx(i,j,k)+betaylow*Bvecy(i,j,k)+betazlow*Bvecz(i,j,k)
+          vdotbeta = betaxlow*velx(i,j,k)+betaylow*vely(i,j,k)+betazlow*velz(i,j,k)
 
 !!$   u0 low is missing the w_lorentz factor (see below)!!
-            utlow = -1.0*alp(i,j,k)
-            btlow = utlow*w_lorentz(i,j,k)*bdotv
-            
-         end if
+          utlow = -1.d0*alp(i,j,k) + vdotbeta
+          
+          btlow = -1.0d0*w_lorentz(i,j,k)*alp(i,j,k)*bdotv + &
+               bdotbeta/w_lorentz(i,j,k) + w_lorentz(i,j,k)*bdotv*vdotbeta
 
 !!$       Calculate the specific relativistic enthalpy times rho + the mag. field contribution times the
 !!$       square of the lorentz factor.