Cleanup of PPM routine. Remove unused PPM parameter. Only apply stricter "beyond Colella" limiting scheme to epsilon (to avoid that correction in one step leads to negative values).

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

3 files changed, 138 insertions, 107 deletions

diff --git a/param.ccl b/param.ccl
index d3ecbdf..f5ad996 100644
--- a/param.ccl
+++ b/param.ccl
@@ -208,11 +208,6 @@ real enhanced_ppm_C3 "Parameter for enhancecd ppm limiter. Default from McCorquo
    0:* :: "must be greater than 0."
 } 0.1
 
-real GRHydro_ppm_atmo_tolerance "A relative value for rho below which we switch to standard ppm. This may be necessary for very sharp transitions in the density to atmosphere values." STEERABLE=RECOVER
-{
-   0:* :: "anything positive"
-} 1e4
-
 
 int eno_order "The order of accuracy of the ENO reconstruction"
 {
diff --git a/src/GRHydro_PPM.F90 b/src/GRHydro_PPM.F90
index 6e9c37d..b58f62a 100644
--- a/src/GRHydro_PPM.F90
+++ b/src/GRHydro_PPM.F90
@@ -57,7 +57,7 @@ subroutine SimplePPM_1d(handle,poly,&
      velxminus,velyminus,velzminus,epsminus,rhoplus,velxplus,velyplus,&
      velzplus,epsplus,trivial_rp, hydro_excision_mask,&
      gxx, gxy, gxz, gyy, gyz, gzz, beta, alp, w_lorentz, &
-     dir, ni, nj, nrx, nry, nrz, ev_l, ev_r, xw, rho_min)
+     dir, ni, nj, nrx, nry, nrz, ev_l, ev_r, xw)
 
   USE GRHydro_Scalars
   USE GRHydro_Eigenproblem
@@ -101,23 +101,21 @@ subroutine SimplePPM_1d(handle,poly,&
   CCTK_REAL :: dupw, dloc, delta
   CCTK_REAL :: agxx, agxy, agxz, agyy, agyz, agzz
   CCTK_REAL, dimension(nx) :: xwind, l_ev_l, l_ev_r
-  !CCTK_INT, dimension(nx) :: use_std_ppm
 
-  CCTK_REAL :: D2a, D2aL, D2aR, D2aC, D2aLim, rhi, threshold, rho_min, daplus, daminus, D3a, D3aLL, D3aL, D3aR, D2aLL, D2aRR, D3aMin, D3aMax
+  CCTK_REAL :: D2a, D2aL, D2aR, D2aC, D2aLim, rhi, daplus, daminus, D3a, D3aLL, D3aL, D3aR, D2aLL, D2aRR, D3aMin, D3aMax
 
   logical :: cond
 
 
 !!$  Initially, all the Riemann problems will be trivial
-
 trivial_rp = .true.
-!use_std_ppm = 1
+
 
 #define STEEP(x,dx,dmx)                                              \
-	 if ( (x(i+1) - x(i)) * (x(i) - x(i-1)) > 0 ) then           &&\
+	 if ( (x(i+1) - x(i)) * (x(i) - x(i-1)) > 0.d0 ) then           &&\
 	    dmx(i) = sign(one, dx(i)) *                                   \
-	       min(abs(dx(i)), 2 * abs(x(i) - x(i-1)),                \
-				 2 * abs(x(i+1) - x(i)))              &&\
+	       min(abs(dx(i)), 2.d0 * abs(x(i) - x(i-1)),                \
+				 2.d0 * abs(x(i+1) - x(i)))              &&\
 	 else                                                           &&\
 	    dmx(i) = 0.d0                                                &&\
 	 end if
@@ -195,11 +193,6 @@ trivial_rp = .true.
 #define APPROX_AT_CELL_INTERFACE(a, ah)  \
       ah = 7.0d0/12.0d0*(a(i)+a(i+1)) - 1.0d0/12.0d0*(a(i-1)+a(i+2))
       
-      !! A threshold value of rho
-      threshold = rho_min * GRHydro_ppm_atmo_tolerance
-
-! .and. rho(i-2) .gt. threshold .and. rho(i+3) .gt. threshold
-
 #define LIMIT(a,ah,C, alim) \
       if ((min(a(i),a(i+1)) .le. ah) .and. (ah .le. max(a(i),a(i+1)))) then       &&\
 	 alim = ah   &&\
@@ -208,24 +201,34 @@ trivial_rp = .true.
 	 D2aL =         (a(i-1) - 2.0d0*a(i)   + a(i+1))                                     &&\
 	 D2aR =         (a(i)   - 2.0d0*a(i+1) + a(i+2))                                     &&\
 	 D2aLim = sign(one, D2a)*min(C*abs(D2aL), C*abs(D2aR), abs(D2a))/3.0d0             &&\
-	 if (D2a*D2aR .ge. 0 .and. D2a*D2aL .ge. 0 .and. abs(D2aLim) .lt. abs(0.5d0*(a(i)+a(i+1)))) then                     &&\
+	 if (D2a*D2aR .ge. 0 .and. D2a*D2aL .ge. 0) then                     &&\
 	    alim = 0.5d0*(a(i)+a(i+1)) - D2aLim &&\
 	 else                                                                             &&\
 	    alim = 0.5d0*(a(i)+a(i+1))                                                    &&\
 	 end if                                                                           &&\
       endif
 
-! .and. abs(D2aLim) .lt. abs(0.5d0*(a(i)+a(i+1)))
+#define LIMIT_EPS(a,ah,C, alim) \
+      if ((min(a(i),a(i+1)) .le. ah) .and. (ah .le. max(a(i),a(i+1)))) then       &&\
+	 alim = ah   &&\
+      else &&\
+	 D2a  = 3.0d0 * (a(i)   - 2.0d0*ah     + a(i+1))                                     &&\
+	 D2aL =         (a(i-1) - 2.0d0*a(i)   + a(i+1))                                     &&\
+	 D2aR =         (a(i)   - 2.0d0*a(i+1) + a(i+2))                                     &&\
+	 D2aLim = sign(one, D2a)*min(C*abs(D2aL), C*abs(D2aR), abs(D2a))/3.0d0             &&\
+	 if (D2a*D2aR .ge. 0 .and. D2a*D2aL .ge. 0 .and. abs(D2aLim) .lt. abs(0.5d0*(a(i)+a(i+1)))) then                     &&\
+	    alim = 0.5d0*(a(i)+a(i+1)) - D2aLim &&\
+	 else                                                                             &&\
+	    alim = 0.5d0*(a(i)+a(i+1))                                                    &&\
+	 end if                                                                           &&\
+      endif
 
-! .and. abs(D2aLim) .lt. abs(0.5d0*(a(i)+a(i+1)))
-!sign(1.0d0, D2a)*min(abs(D2aLim), abs(0.5*(a(i)+a(i+1))))
 
    if (PPM3) then
       
       !! We initialize "plus" \equiv a_j+1/2 with (16) via APPROX_AT_CELL_INTERFACE, 
       !! then checking for (13) of Colella & Sekora 2008 and applying
       !! (18) and (19) if (13) is not satisfied. This is done with LIMIT.
-      !! There, we also switch to lower order if we are near atmosphere values.
       do i = 2, nx-2
          APPROX_AT_CELL_INTERFACE(rho, rhoplus(i))
          LIMIT(rho, rhoplus(i), enhanced_ppm_C2, rhoplus(i))
@@ -246,7 +249,7 @@ trivial_rp = .true.
       if (poly .eq. 0) then
 	 do i = 2, nx-2
 	    APPROX_AT_CELL_INTERFACE(eps, epsplus(i))
-	    LIMIT(eps, epsplus(i), enhanced_ppm_C2, epsplus(i))
+	    LIMIT_EPS(eps, epsplus(i), enhanced_ppm_C2, epsplus(i))
 	    epsminus(i+1) = epsplus(i)
 	 end do
       end if
@@ -559,6 +562,60 @@ trivial_rp = .true.
 	    D3aMin = min(D3aLL, D3aL, D3a, D3aR)   &&\
 	    D3aMax = max(D3aLL, D3aL, D3a, D3aR)   &&\
 	    if (C3 * max(abs(D3aMin), abs(D3aMax)) .le. D3aMax-D3aMin) then &&\
+		  if (daplus*daminus .lt. 0) then        &&\
+		     aplus(i)  = a(i) + daplus * rhi     &&\
+		     aminus(i) = a(i) - daminus * rhi    &&\
+		  else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then   &&\
+		     aminus(i)  = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus   &&\
+		  else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\
+		     aplus(i)  = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus   &&\
+		  endif &&\
+	    endif &&\
+	 endif &&\
+	 trivial_rp(i-1) = .false.                                 &&\
+	 trivial_rp(i) = .false.                                   &&\
+      else                                                         &&\
+	    if (abs(daplus) .ge. 2.0d0*abs(daminus)) then     &&\
+	       aplus(i) = a(i) + 2.0d0*daminus               &&\
+	    else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\
+	       aminus(i) = a(i) - 2.0d0*daplus               &&\
+	    end if                                                    &&\
+	 trivial_rp(i-1) = .false.  &&\
+	 trivial_rp(i) = .false.    &&\
+      endif
+
+!! Monotonicity of PPM 2011 (McCorquodale & Colella 2011), Sec. 2.4.1, Eq. 23-34.
+!! This contains an additional limiter in case the correction becomes larger than
+!! the corrected value. This is to avoid negative values for epsilon.
+!! This requires 4 stencil points.
+#define MON_WITH_LOCAL_EXTREMUM_STENCIL4_EPS(aminus, a, aplus, C, C3) \
+      daplus = aplus(i)-a(i)   &&\
+      daminus = a(i)-aminus(i) &&\
+      if (daplus*daminus .le. 0 .or. (a(i-2)-a(i))*(a(i)-a(i+2)) .le. 0) then &&\
+	 D2a  = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i)))                                           && \
+	 D2aC = a(i-1) - 2.0d0*a(i)   + a(i+1)                                                     &&\
+	 D2aL = a(i-2) - 2.0d0*a(i-1) + a(i)                                                       &&\
+	 D2aLL = a(i-3) - 2.0d0*a(i-2) + a(i-1)                                                       &&\
+	 D2aR = a(i)   - 2.0d0*a(i+1) + a(i+2)                                                     &&\
+	 D2aRR = a(i+1)   - 2.0d0*a(i+2) + a(i+3)                                                     &&\
+	 D3a = D2aR - D2aC   &&\
+	 D3aL = D2aC - D2aL  &&\
+	 D3aR = D2aRR - D2aR &&\
+	 D3aLL = D2aL - D2aLL &&\
+	 if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\
+	    D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a))  &&\
+	 else                                                      &&\
+	    D2aLim = 0  &&\
+	 end if                                                    &&\
+	 if (abs(D2a) .le. 1.d-12*max(abs(a(i-2)), abs(a(i-1)), abs(a(i)), abs(a(i+1)), abs(a(i+2)))) then  &&\
+	    rhi = 0  &&\
+	 else        &&\
+	    rhi = D2aLim / D2a   &&\
+	 endif  &&\
+	 if (.not. (rhi .ge. 1.0d0 - 1.d-12)) then   &&\
+	    D3aMin = min(D3aLL, D3aL, D3a, D3aR)   &&\
+	    D3aMax = max(D3aLL, D3aL, D3a, D3aR)   &&\
+	    if (C3 * max(abs(D3aMin), abs(D3aMax)) .le. D3aMax-D3aMin) then &&\
 	       if (abs(daplus) .le. abs(a(i)) .and. abs(daminus) .le. abs(a(i))) then &&\
 		  if (daplus*daminus .lt. 0) then        &&\
 		     aplus(i)  = a(i) + daplus * rhi     &&\
@@ -599,35 +656,8 @@ trivial_rp = .true.
 	 trivial_rp(i-1) = .false.  &&\
 	 trivial_rp(i) = .false.    &&\
       endif
-            
 
-! We use (20) of Colella & Sekora 2008 to check wether we are at a local maximum.
-! If so, we make use of (22) and (23) to preserve accuracy there. Otherwise we use
-! enhanced monotocinity preservation (26)
-! #define MON_WITH_LOCAL_EXTREMUM(aminus, a, aplus, C) \
-! 	 if ((aplus(i)-a(i))*(a(i)-aminus(i)) .le. 0 .or. (a(i-1)-a(i))*(a(i)-a(i+1)) .le. 0) then &&\
-! 	    D2a  = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i)))                                           && \
-! 	    D2aC = a(i-1) - 2.0d0*a(i)   + a(i+1)                                                     &&\
-! 	    D2aL = a(i-2) - 2.0d0*a(i-1) + a(i)                                                       &&\
-! 	    D2aR = a(i)   - 2.0d0*a(i+1) + a(i+2)                                                     &&\
-! 	    if (D2a .ne. 0 .and. sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\
-! 	       aplus(i)  = a(i) + (aplus(i) -a(i)) * sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) / D2a                   &&\
-! 	       aminus(i) = a(i) + (aminus(i)-a(i)) * sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) / D2a                   &&\
-! 	    else                                                      &&\
-! 	       aplus(i)  = a(i)                                       &&\
-! 	       aminus(i) = a(i)                                       &&\
-! 	    end if                                                    &&\
-! 	    trivial_rp(i-1) = .false.                                 &&\
-! 	    trivial_rp(i) = .false.                                   &&\
-! 	 else                                                         &&\
-! 	    if (abs(aplus(i)-a(i)) .ge. 2.0d0*abs(aminus(i)-a(i))) then     &&\
-! 	       aplus(i) = a(i) - 2.0d0*(aminus(i) - a(i))               &&\
-! 	    else if (abs(aminus(i)-a(i)) .ge. 2.0d0*abs(aplus(i)-a(i))) then &&\
-! 	       aminus(i) = a(i) - 2.0d0*(aplus(i) - a(i))               &&\
-! 	    end if                                                    &&\
-! 	    trivial_rp(i-1) = .false.  &&\
-! 	    trivial_rp(i) = .false.    &&\
-! 	 endif
+
 
 !! Monotonicity of PPM 2011 (McCorquodale & Colella 2011), Sec. 2.4.1, Eq. 23-34.
 !! This does not use the check for deviations from a cubic. Thus, it gets away with only 3 stencil points!
@@ -650,6 +680,51 @@ trivial_rp = .true.
 	    rhi = D2aLim / D2a   &&\
 	 endif  &&\
 	 if (.not. (rhi .ge. 1.0d0 - 1.d-12)) then   &&\
+	       if (daplus*daminus .lt. 0) then        &&\
+		  aplus(i)  = a(i) + daplus * rhi     &&\
+		  aminus(i) = a(i) - daminus * rhi    &&\
+	       else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then   &&\
+		  aminus(i)  = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus   &&\
+	       else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\
+		  aplus(i)  = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus   &&\
+	       endif &&\
+	 endif &&\
+	 trivial_rp(i-1) = .false.                                 &&\
+	 trivial_rp(i) = .false.                                   &&\
+      else                                                         &&\
+	    if (abs(daplus) .ge. 2.0d0*abs(daminus)) then     &&\
+	       aplus(i) = a(i) + 2.0d0*daminus               &&\
+	    else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\
+	       aminus(i) = a(i) - 2.0d0*daplus               &&\
+	    end if                                                    &&\
+	 trivial_rp(i-1) = .false.  &&\
+	 trivial_rp(i) = .false.    &&\
+      endif   
+      
+
+!! Monotonicity of PPM 2011 (McCorquodale & Colella 2011), Sec. 2.4.1, Eq. 23-34.
+!! This does not use the check for deviations from a cubic. Thus, it gets away with only 3 stencil points!
+!! This contains an additional limiter in case the correction becomes larger than
+!! the corrected value. This is to avoid negative values for epsilon.
+#define MON_WITH_LOCAL_EXTREMUM_EPS(aminus, a, aplus, C) \
+      daplus = aplus(i)-a(i)   &&\
+      daminus = a(i)-aminus(i) &&\
+      if (daplus*daminus .le. 0 .or. (a(i-2)-a(i))*(a(i)-a(i+2)) .le. 0) then &&\
+	 D2a  = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i)))                                           && \
+	 D2aC = a(i-1) - 2.0d0*a(i)   + a(i+1)                                                     &&\
+	 D2aL = a(i-2) - 2.0d0*a(i-1) + a(i)                                                       &&\
+	 D2aR = a(i)   - 2.0d0*a(i+1) + a(i+2)                                                     &&\
+	 if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\
+	    D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a))  &&\
+	 else                                                      &&\
+	    D2aLim = 0  &&\
+	 end if                                                    &&\
+	 if (abs(D2a) .le. 1.d-12*max(abs(a(i-2)), abs(a(i-1)), abs(a(i)), abs(a(i+1)), abs(a(i+2)))) then  &&\
+	    rhi = 0  &&\
+	 else        &&\
+	    rhi = D2aLim / D2a   &&\
+	 endif  &&\
+	 if (.not. (rhi .ge. 1.0d0 - 1.d-12)) then   &&\
 	    if (abs(daplus) .le. abs(a(i)) .and. abs(daminus) .le. abs(a(i))) then &&\
 	       if (daplus*daminus .lt. 0) then        &&\
 		  aplus(i)  = a(i) + daplus * rhi     &&\
@@ -689,8 +764,8 @@ trivial_rp = .true.
 	 trivial_rp(i-1) = .false.  &&\
 	 trivial_rp(i) = .false.    &&\
       endif   
-      
-   
+
+
 
    if (use_enhanced_ppm .eq. 1) then
       !! Constrain parabolic profiles, PPM 2011/2008
@@ -703,7 +778,7 @@ trivial_rp = .true.
 	    end do
 	    if (poly .eq. 0) then
 		  do i = 3, nx - 2
-		     MON_WITH_LOCAL_EXTREMUM(epsminus,eps,epsplus, enhanced_ppm_C2)
+		     MON_WITH_LOCAL_EXTREMUM_EPS(epsminus,eps,epsplus, enhanced_ppm_C2)
 		  end do
 	    end if
 	    
@@ -716,7 +791,7 @@ trivial_rp = .true.
 	    end do
 	    if (poly .eq. 0) then
 		  do i = 4, nx - 3
-		     MON_WITH_LOCAL_EXTREMUM_STENCIL4(epsminus,eps,epsplus, enhanced_ppm_C2, enhanced_ppm_C3)
+		     MON_WITH_LOCAL_EXTREMUM_STENCIL4_EPS(epsminus,eps,epsplus, enhanced_ppm_C2, enhanced_ppm_C3)
 		  end do
 	    endif
       end if 
@@ -744,7 +819,7 @@ trivial_rp = .true.
 	 end do
       else  !!$ Really implement C&W, page 197; which requires stencil 4.
 	 do i = 4, nx - 3
-	    s=sign(1.d0, -dpress(i))
+	    s=sign(one, -dpress(i))
 	    flatten = max(tilde_flatten(i), tilde_flatten(i+s))  
 	    if (abs(1.d0 - flatten) > 0.d0) then
 	       trivial_rp(i-1) = .false.
@@ -1037,7 +1112,7 @@ subroutine SimplePPM_tracer_1d(nx,dx,rho,velx,vely,velz, &
   else  !!$ Really implement C&W, page 197; which requires stencil 4.
      do itracer=1,number_of_tracers
         do i = 4, nx - 3
-           s=sign(1.d0, -dpress(i))
+           s=sign(one, -dpress(i))
            flatten = max(tilde_flatten(i), tilde_flatten(i+s))  
            tracerplus(i,itracer) = flatten * tracerplus(i,itracer) + &
                 (1.d0 - flatten) * tracer(i,itracer)
@@ -1130,7 +1205,7 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, &
   CCTK_REAL :: dpress2,w,flatten,dvel
   CCTK_REAL :: eta, etatilde
 
-  CCTK_REAL :: D2a, D2aL, D2aR, D2aC, D2aLim, rhi, threshold, rho_min, daplus, daminus, D3a, D3aLL, D3aL, D3aR, D2aLL, D2aRR, D3aMin, D3aMax
+  CCTK_REAL :: D2a, D2aL, D2aR, D2aC, D2aLim, rhi, daplus, daminus, D3a, D3aLL, D3aL, D3aR, D2aLL, D2aRR, D3aMin, D3aMax
 
   
   do i = 2, nx - 1
@@ -1199,7 +1274,7 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, &
 	 D2a  = 3.0d0 * (a(i)   - 2.0d0*ah     + a(i+1))                                     &&\
 	 D2aL =         (a(i-1) - 2.0d0*a(i)   + a(i+1))                                     &&\
 	 D2aR =         (a(i)   - 2.0d0*a(i+1) + a(i+2))                                     &&\
-	 if (D2a*D2aR .ge. 0 .and. D2a*D2aL .ge. 0 .and. abs(D2aLim) .lt. abs(0.5d0*(a(i)+a(i+1)))) then                     &&\
+	 if (D2a*D2aR .ge. 0 .and. D2a*D2aL .ge. 0) then                     &&\
 	    alim = 0.5d0*(a(i)+a(i+1)) - sign(one, D2a)*min(C*abs(D2aL), C*abs(D2aR), abs(D2a))/3.0d0 &&\
 	 else                                                                             &&\
 	    alim = 0.5d0*(a(i)+a(i+1))                                                    &&\
@@ -1270,7 +1345,7 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, &
 	 end do
       else  !!$ Really implement C&W, page 197; which requires stencil 4.
 	 do i = 4, nx - 3
-	    s=sign(1.d0, -dpress(i))
+	    s=sign(one, -dpress(i))
 	    flatten = max(tilde_flatten(i), tilde_flatten(i+s))  
 	    Y_e_plus(i) = flatten * Y_e_plus(i) + &
 		  (1.d0 - flatten) * Y_e(i)
@@ -1312,7 +1387,6 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, &
 	    D3aMin = min(D3aLL, D3aL, D3a, D3aR)   &&\
 	    D3aMax = max(D3aLL, D3aL, D3a, D3aR)   &&\
 	    if (C3 * max(abs(D3aMin), abs(D3aMax)) .le. D3aMax-D3aMin) then &&\
-	       if (abs(daplus) .le. abs(a(i)) .and. abs(daminus) .le. abs(a(i))) then &&\
 		  if (daplus*daminus .lt. 0) then        &&\
 		     aplus(i)  = a(i) + daplus * rhi     &&\
 		     aminus(i) = a(i) - daminus * rhi    &&\
@@ -1321,32 +1395,14 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, &
 		  else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\
 		     aplus(i)  = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus   &&\
 		  endif &&\
-	       else &&\
-		  if (daplus*daminus .lt. 0) then        &&\
-		     aplus(i)  = a(i) &&\
-		     aminus(i) = a(i) &&\
-		  else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then   &&\
-		     aminus(i)  = a(i) &&\
-		  else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\
-		     aplus(i)  = a(i) &&\
-		  endif &&\
-	       endif &&\
 	    endif &&\
 	 endif &&\
       else                                                         &&\
-	 if (abs(daplus) .le. abs(a(i)) .and. abs(daminus) .le. abs(a(i))) then &&\
 	    if (abs(daplus) .ge. 2.0d0*abs(daminus)) then     &&\
 	       aplus(i) = a(i) + 2.0d0*daminus               &&\
 	    else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\
 	       aminus(i) = a(i) - 2.0d0*daplus               &&\
 	    end if                                                    &&\
-	 else &&\
-	    if (abs(daplus) .ge. 2.0d0*abs(daminus)) then     &&\
-	       aplus(i) = a(i)               &&\
-	    else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\
-	       aminus(i) = a(i)               &&\
-	    end if                                                    &&\
-	 endif &&\
       endif
             
 
@@ -1375,7 +1431,6 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, &
 	    rhi = D2aLim / D2a   &&\
 	 endif  &&\
 	 if (.not. (rhi .ge. 1.0d0 - 1.d-12)) then   &&\
-	    if (abs(daplus) .le. abs(a(i)) .and. abs(daminus) .le. abs(a(i))) then &&\
 	       if (daplus*daminus .lt. 0) then        &&\
 		  aplus(i)  = a(i) + daplus * rhi     &&\
 		  aminus(i) = a(i) - daminus * rhi    &&\
@@ -1384,31 +1439,13 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, &
 	       else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\
 		  aplus(i)  = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus   &&\
 	       endif &&\
-	    else &&\
-	       if (daplus*daminus .lt. 0) then        &&\
-		  aplus(i)  = a(i) &&\
-		  aminus(i) = a(i) &&\
-	       else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then   &&\
-		  aminus(i)  = a(i) &&\
-	       else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\
-		  aplus(i)  = a(i) &&\
-	       endif &&\
-	    endif &&\
 	 endif      &&\
       else                                                         &&\
-	 if (abs(daplus) .le. abs(a(i)) .and. abs(daminus) .le. abs(a(i))) then &&\
 	    if (abs(daplus) .ge. 2.0d0*abs(daminus)) then     &&\
 	       aplus(i) = a(i) - 2.0d0*(aminus(i) - a(i))               &&\
 	    else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\
 	       aminus(i) = a(i) - 2.0d0*(aplus(i) - a(i))               &&\
 	    end if  &&\
-	 else &&\
-	    if (abs(daplus) .ge. 2.0d0*abs(daminus)) then     &&\
-	       aplus(i) = a(i)               &&\
-	    else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\
-	       aminus(i) = a(i)               &&\
-	    end if                                                    &&\
-	 endif &&\
       endif
 
 
@@ -1438,7 +1475,7 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, &
 	 end do
       else  !!$ Really implement C&W, page 197; which requires stencil 4.
 	 do i = 4, nx - 3
-	    s=sign(1.d0, -dpress(i))
+	    s=sign(one, -dpress(i))
 	    flatten = max(tilde_flatten(i), tilde_flatten(i+s))  
 	    Y_e_plus(i) = flatten * Y_e_plus(i) + &
 		  (1.d0 - flatten) * Y_e(i)
diff --git a/src/GRHydro_PPMReconstruct_drv.F90 b/src/GRHydro_PPMReconstruct_drv.F90
index 3b51e2d..02486d0 100644
--- a/src/GRHydro_PPMReconstruct_drv.F90
+++ b/src/GRHydro_PPMReconstruct_drv.F90
@@ -151,7 +151,6 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS)
   enddo
   !$OMP END PARALLEL DO
 
-
 !!$ PPM starts:
     if (flux_direction == 1) then
          !$OMP PARALLEL DO PRIVATE(i, j, k)
@@ -168,7 +167,7 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS)
                 w_lorentz(:,j,k), &
                 1, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_x_left, &
                 GRHydro_mppm_eigenvalue_x_right, &
-                GRHydro_mppm_xwind, GRHydro_rho_min)
+                GRHydro_mppm_xwind)
              do i = 1, nx
                 if (trivial_rp(i,j,k)) then
                   SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, trivialx)
@@ -198,7 +197,7 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS)
                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), GRHydro_rho_min)
+                    press(:,j,k))
             end do
          end do
          !$OMP END PARALLEL DO
@@ -219,7 +218,7 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS)
                  w_lorentz(j,:,k), &
                  2, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_y_left, &
                  GRHydro_mppm_eigenvalue_y_right, &
-                 GRHydro_mppm_xwind, GRHydro_rho_min)
+                 GRHydro_mppm_xwind)
             do i = 1, ny
               if (trivial_rp(j,i,k)) then
                  SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, trivialy)
@@ -249,7 +248,7 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS)
                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), GRHydro_rho_min)
+                    press(j,:,k))
             end do
          end do
          !$OMP END PARALLEL DO
@@ -270,7 +269,7 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS)
                  w_lorentz(j,k,:), &
                  3, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_z_left, &
                  GRHydro_mppm_eigenvalue_z_right, &
-                 GRHydro_mppm_xwind, GRHydro_rho_min)
+                 GRHydro_mppm_xwind)
             do i = 1, nz
               if (trivial_rp(j,k,i)) then
                 SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, trivialz)
@@ -301,7 +300,7 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS)
                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,:), GRHydro_rho_min)
+                    press(j,k,:))
             end do
          end do
          !$OMP END PARALLEL DO