diff --git a/Exec/science/subch_planar/Problem_Derive.H b/Exec/science/subch_planar/Problem_Derive.H
index 075c1828db..2110a535a1 100644
--- a/Exec/science/subch_planar/Problem_Derive.H
+++ b/Exec/science/subch_planar/Problem_Derive.H
@@ -3,3 +3,8 @@ void ca_dergradpoverp
   (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
    const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
    amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
+
+void ca_dergradpoverp1
+  (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
+   const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+   amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
diff --git a/Exec/science/subch_planar/Problem_Derive.cpp b/Exec/science/subch_planar/Problem_Derive.cpp
index 3cb25d2cb4..6a67c0532f 100644
--- a/Exec/science/subch_planar/Problem_Derive.cpp
+++ b/Exec/science/subch_planar/Problem_Derive.cpp
@@ -12,207 +12,401 @@ void ca_dergradpoverp(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*nco
                       Real /*time*/, const int* /*bcrec*/, int /*level*/)
 {
 
-  // derive the dynamic pressure
+    // Compute grad p . U / (p |U|) -- this is what we do in the shock
+    // detection algorithm
 
-  const auto dx = geomdata.CellSizeArray();
-  const int coord_type = geomdata.Coord();
+    const auto dx = geomdata.CellSizeArray();
+    const int coord_type = geomdata.Coord();
 
-  auto const dat = datfab.array();
-  auto const der = derfab.array();
+    auto const dat = datfab.array();
+    auto const der = derfab.array();
 
 #if AMREX_SPACEDIM == 3
-  amrex::Error("3D not supported");
+    amrex::Error("3D not supported");
 #endif
 #if AMREX_SPACEDIM == 1
-  amrex::Error("1D not supported");
+    amrex::Error("1D not supported");
 #endif
 
-  Real dxinv = 1.0_rt / dx[0];
-  Real dyinv = 1.0_rt / dx[1];
+    Real dxinv = 1.0_rt / dx[0];
+    Real dyinv = 1.0_rt / dx[1];
 
-  // Compute grad p . U / (p |U|) -- this is what we do in the shock
-  // detection algorithm
+    amrex::ParallelFor(bx,
+    [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k)
+    {
 
-  amrex::ParallelFor(bx,
-  [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k)
-  {
+        Real div_u = 0.0_rt;
 
-      Real div_u = 0.0_rt;
+        // get the pressure
+        eos_t eos_state;
 
-      // get the pressure
-      eos_t eos_state;
+        Real rhoInv = 1.0 / dat(i,j,k,URHO);
 
-      Real rhoInv = 1.0 / dat(i,j,k,URHO);
-
-      eos_state.rho = dat(i,j,k,URHO);
-      eos_state.T = dat(i,j,k,UTEMP);
-      eos_state.e = dat(i,j,k,UEINT) * rhoInv;
-      for (int n = 0; n < NumSpec; n++) {
-          eos_state.xn[n] = dat(i,j,k,UFS+n) * rhoInv;
-      }
+        eos_state.rho = dat(i,j,k,URHO);
+        eos_state.T = dat(i,j,k,UTEMP);
+        eos_state.e = dat(i,j,k,UEINT) * rhoInv;
+        for (int n = 0; n < NumSpec; n++) {
+            eos_state.xn[n] = dat(i,j,k,UFS+n) * rhoInv;
+        }
 #if NAUX_NET > 0
-      for (int n = 0; n < NumAux; n++) {
-          eos_state.aux[n] = dat(i,j,k,UFX+n) * rhoInv;
-      }
+        for (int n = 0; n < NumAux; n++) {
+            eos_state.aux[n] = dat(i,j,k,UFX+n) * rhoInv;
+        }
 #endif
 
-      if (eos_state.e <= 0.0_rt) {
-          eos(eos_input_rt, eos_state);
-      } else {
-          eos(eos_input_re, eos_state);
-      }
-      Real p_zone = eos_state.p;
+        if (eos_state.e <= 0.0_rt) {
+            eos(eos_input_rt, eos_state);
+        } else {
+            eos(eos_input_re, eos_state);
+        }
+        Real p_zone = eos_state.p;
 
-      // construct div{U}
-      if (coord_type == 0) {
+        Real up = dat(i+1,j,k,UMX) / dat(i+1,j,k,URHO);
+        Real um = dat(i-1,j,k,UMX) / dat(i-1,j,k,URHO);
+        Real u0 = dat(i,j,k,UMX) / dat(i,j,k,URHO);
 
-          // Cartesian
-          div_u += 0.5_rt * (dat(i+1,j,k,UMX) / dat(i+1,j,k,URHO) -
-                             dat(i-1,j,k,UMX) / dat(i-1,j,k,URHO)) * dxinv;
-#if (AMREX_SPACEDIM >= 2)
-          div_u += 0.5_rt * (dat(i,j+1,k,UMY) / dat(i,j+1,k,URHO) -
-                             dat(i,j-1,k,UMY) / dat(i,j-1,k,URHO)) * dyinv;
-#endif
+        Real vp = dat(i,j+1,k,UMY) / dat(i,j+1,k,URHO);
+        Real vm = dat(i,j-1,k,UMY) / dat(i,j-1,k,URHO);
+        Real v0 = dat(i,j,k,UMY) / dat(i,j,k,URHO);
 
-      } else if (coord_type == 1) {
+        // construct div{U}
+        if (coord_type == 0) {
 
-          // r-z
-          Real rc = (i + 0.5_rt) * dx[0];
-          Real rm = (i - 1 + 0.5_rt) * dx[0];
-          Real rp = (i + 1 + 0.5_rt) * dx[0];
+            // Cartesian
+            div_u += 0.5_rt * (up - um) * dxinv;
+            div_u += 0.5_rt * (vp - vm) * dyinv;
 
-          div_u += 0.5_rt * (rp * dat(i+1,j,k,UMX) / dat(i+1,j,k,URHO) -
-                             rm * dat(i-1,j,k,UMX) / dat(i-1,j,k,URHO)) / (rc * dx[0]);
-          div_u += 0.5_rt * (rp * dat(i+1,j,k,UMX) / dat(i+1,j,k,URHO) -
-                             rm * dat(i-1,j,k,UMX) / dat(i-1,j,k,URHO)) / (rc * dx[0]) +
-              0.5_rt * (dat(i,j+1,k,UMY) / dat(i,j+1,k,URHO) -
-                        dat(i,j-1,k,UMY) / dat(i,j-1,k,URHO)) * dyinv;
+        } else if (coord_type == 1) {
 
-#ifndef AMREX_USE_GPU
+            // r-z
+            Real rc = (i + 0.5_rt) * dx[0];
+            Real rm = (i - 1 + 0.5_rt) * dx[0];
+            Real rp = (i + 1 + 0.5_rt) * dx[0];
 
-      } else {
-          amrex::Error("ERROR: invalid coord_type in shock");
+            div_u += 0.5_rt * (rp * up - rm * um) / (rc * dx[0]) +
+                0.5_rt * (vp - vm) * dyinv;
+
+#ifndef AMREX_USE_GPU
+        } else {
+            amrex::Error("ERROR: invalid coord_type in shock");
 #endif
-      }
+        }
 
 
-      // now compute (grad P - rho g) . dx
+        // we need to compute p in the full stencil
 
-      // we need to compute p in the full stencil
+        Real p_ip1{};
+        {
+            Real lrhoInv = 1.0 / dat(i+1,j,k,URHO);
 
-      Real p_ip1{};
-      {
-          Real lrhoInv = 1.0 / dat(i+1,j,k,URHO);
+            eos_state.rho = dat(i+1,j,k,URHO);
+            eos_state.T = dat(i+1,j,k,UTEMP);
+            eos_state.e = dat(i+1,j,k,UEINT) * lrhoInv;
+            for (int n = 0; n < NumSpec; n++) {
+                eos_state.xn[n] = dat(i+1,j,k,UFS+n) * lrhoInv;
+            }
+#if NAUX_NET > 0
+            for (int n = 0; n < NumAux; n++) {
+                eos_state.aux[n] = dat(i+1,j,k,UFX+n) * lrhoInv;
+            }
+#endif
 
-          eos_state.rho = dat(i+1,j,k,URHO);
-          eos_state.T = dat(i+1,j,k,UTEMP);
-          eos_state.e = dat(i+1,j,k,UEINT) * lrhoInv;
-          for (int n = 0; n < NumSpec; n++) {
-              eos_state.xn[n] = dat(i+1,j,k,UFS+n) * lrhoInv;
-          }
+            if (eos_state.e <= 0.0_rt) {
+                eos(eos_input_rt, eos_state);
+            } else {
+                eos(eos_input_re, eos_state);
+            }
+            p_ip1 = eos_state.p;
+        }
+
+        Real p_im1{};
+        {
+            Real lrhoInv = 1.0 / dat(i-1,j,k,URHO);
+
+            eos_state.rho = dat(i-1,j,k,URHO);
+            eos_state.T = dat(i-1,j,k,UTEMP);
+            eos_state.e = dat(i-1,j,k,UEINT) * lrhoInv;
+            for (int n = 0; n < NumSpec; n++) {
+                eos_state.xn[n] = dat(i-1,j,k,UFS+n) * lrhoInv;
+            }
 #if NAUX_NET > 0
-          for (int n = 0; n < NumAux; n++) {
-              eos_state.aux[n] = dat(i+1,j,k,UFX+n) * lrhoInv
-          }
+            for (int n = 0; n < NumAux; n++) {
+                eos_state.aux[n] = dat(i-1,j,k,UFX+n) * lrhoInv;
+            }
 #endif
 
-          if (eos_state.e <= 0.0_rt) {
-              eos(eos_input_rt, eos_state);
-          } else {
-              eos(eos_input_re, eos_state);
-          }
-          p_ip1 = eos_state.p;
+            if (eos_state.e <= 0.0_rt) {
+                eos(eos_input_rt, eos_state);
+            } else {
+                eos(eos_input_re, eos_state);
+            }
+            p_im1 = eos_state.p;
       }
 
-      Real p_im1{};
-      {
-          Real lrhoInv = 1.0 / dat(i-1,j,k,URHO);
+        Real dP_x = 0.5_rt * (p_ip1 - p_im1);
+
+        Real p_jp1{};
+        {
+            Real lrhoInv = 1.0 / dat(i,j+1,k,URHO);
 
-          eos_state.rho = dat(i-1,j,k,URHO);
-          eos_state.T = dat(i-1,j,k,UTEMP);
-          eos_state.e = dat(i-1,j,k,UEINT) * lrhoInv;
-          for (int n = 0; n < NumSpec; n++) {
-              eos_state.xn[n] = dat(i-1,j,k,UFS+n) * lrhoInv;
-          }
+            eos_state.rho = dat(i,j+1,k,URHO);
+            eos_state.T = dat(i,j+1,k,UTEMP);
+            eos_state.e = dat(i,j+1,k,UEINT) * lrhoInv;
+            for (int n = 0; n < NumSpec; n++) {
+                eos_state.xn[n] = dat(i,j+1,k,UFS+n) * lrhoInv;
+            }
 #if NAUX_NET > 0
-          for (int n = 0; n < NumAux; n++) {
-              eos_state.aux[n] = dat(i-1,j,k,UFX+n) * lrhoInv
-          }
+            for (int n = 0; n < NumAux; n++) {
+                eos_state.aux[n] = dat(i,j+1,k,UFX+n) * lrhoInv;
+            }
 #endif
 
-          if (eos_state.e <= 0.0_rt) {
-              eos(eos_input_rt, eos_state);
-          } else {
-              eos(eos_input_re, eos_state);
-          }
-          p_im1 = eos_state.p;
-      }
+            if (eos_state.e <= 0.0_rt) {
+                eos(eos_input_rt, eos_state);
+            } else {
+                eos(eos_input_re, eos_state);
+            }
+            p_jp1 = eos_state.p;
+        }
+
+        Real p_jm1{};
+        {
+            Real lrhoInv = 1.0 / dat(i,j-1,k,URHO);
+
+            eos_state.rho = dat(i,j-1,k,URHO);
+            eos_state.T = dat(i,j-1,k,UTEMP);
+            eos_state.e = dat(i,j-1,k,UEINT) * lrhoInv;
+            for (int n = 0; n < NumSpec; n++) {
+                eos_state.xn[n] = dat(i,j-1,k,UFS+n) * lrhoInv;
+            }
+#if NAUX_NET > 0
+            for (int n = 0; n < NumAux; n++) {
+                eos_state.aux[n] = dat(i,j-1,k,UFX+n) * lrhoInv;
+            }
+#endif
+
+            if (eos_state.e <= 0.0_rt) {
+                eos(eos_input_rt, eos_state);
+            } else {
+                eos(eos_input_re, eos_state);
+            }
+            p_jm1 = eos_state.p;
+        }
 
-      Real dP_x = 0.5_rt * (p_ip1 - p_im1);
+        Real dP_y = 0.5_rt * (p_jp1 - p_jm1);
 
-      Real p_jp1{};
-      {
-          Real lrhoInv = 1.0 / dat(i,j+1,k,URHO);
+        //Real gradPdx_over_P = std::sqrt(dP_x * dP_x + dP_y * dP_y + dP_z * dP_z) / dat(i,j,k,QPRES);
+        Real vel = std::sqrt(u0 * u0 + v0 * v0);
+
+        Real gradPdx_over_P{0.0_rt};
+        if (vel != 0.0) {
+            gradPdx_over_P = std::abs(dP_x * u0 + dP_y * v0) / vel;
+        }
+        gradPdx_over_P /= p_zone;
+
+        der(i,j,k,0) = gradPdx_over_P;
+
+    });
+
+}
+
+void ca_dergradpoverp1(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
+                       const FArrayBox& datfab, const Geometry& geomdata,
+                       Real /*time*/, const int* /*bcrec*/, int /*level*/)
+{
+
+    // compute grad p projected onto div{U} as an alternative to the shock detection
+
+    const auto dx = geomdata.CellSizeArray();
+    const int coord_type = geomdata.Coord();
+
+    auto const dat = datfab.array();
+    auto const der = derfab.array();
+
+#if AMREX_SPACEDIM == 3
+    amrex::Error("3D not supported");
+#endif
+#if AMREX_SPACEDIM == 1
+    amrex::Error("1D not supported");
+#endif
 
-          eos_state.rho = dat(i,j+1,k,URHO);
-          eos_state.T = dat(i,j+1,k,UTEMP);
-          eos_state.e = dat(i,j+1,k,UEINT) * lrhoInv;
-          for (int n = 0; n < NumSpec; n++) {
-              eos_state.xn[n] = dat(i,j+1,k,UFS+n) * lrhoInv;
-          }
+    Real dxinv = 1.0_rt / dx[0];
+    Real dyinv = 1.0_rt / dx[1];
+
+
+    amrex::ParallelFor(bx,
+    [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k)
+    {
+
+        Real div_u = 0.0_rt;
+
+        // get the pressure
+        eos_t eos_state;
+
+        Real rhoInv = 1.0 / dat(i,j,k,URHO);
+
+        eos_state.rho = dat(i,j,k,URHO);
+        eos_state.T = dat(i,j,k,UTEMP);
+        eos_state.e = dat(i,j,k,UEINT) * rhoInv;
+        for (int n = 0; n < NumSpec; n++) {
+            eos_state.xn[n] = dat(i,j,k,UFS+n) * rhoInv;
+        }
 #if NAUX_NET > 0
-          for (int n = 0; n < NumAux; n++) {
-              eos_state.aux[n] = dat(i,j+1,k,UFX+n) * lrhoInv
-          }
+        for (int n = 0; n < NumAux; n++) {
+            eos_state.aux[n] = dat(i,j,k,UFX+n) * rhoInv;
+        }
 #endif
 
-          if (eos_state.e <= 0.0_rt) {
-              eos(eos_input_rt, eos_state);
-          } else {
-              eos(eos_input_re, eos_state);
-          }
-          p_jp1 = eos_state.p;
-      }
+        if (eos_state.e <= 0.0_rt) {
+            eos(eos_input_rt, eos_state);
+        } else {
+            eos(eos_input_re, eos_state);
+        }
+        Real p_zone = eos_state.p;
+
+        Real up = dat(i+1,j,k,UMX) / dat(i+1,j,k,URHO);
+        Real um = dat(i-1,j,k,UMX) / dat(i-1,j,k,URHO);
+        Real u0 = dat(i,j,k,UMX) / dat(i,j,k,URHO);
+
+        Real vp = dat(i,j+1,k,UMY) / dat(i,j+1,k,URHO);
+        Real vm = dat(i,j-1,k,UMY) / dat(i,j-1,k,URHO);
+        Real v0 = dat(i,j,k,UMY) / dat(i,j,k,URHO);
+
+        // construct div{U}
+        if (coord_type == 0) {
+
+            // Cartesian
+            div_u += 0.5_rt * (up - um) * dxinv;
+            div_u += 0.5_rt * (vp - vm) * dyinv;
+
+        } else if (coord_type == 1) {
+
+            // r-z
+            Real rc = (i + 0.5_rt) * dx[0];
+            Real rm = (i - 1 + 0.5_rt) * dx[0];
+            Real rp = (i + 1 + 0.5_rt) * dx[0];
+
+            div_u += 0.5_rt * (rp * up - rm * um) / (rc * dx[0]) +
+                0.5_rt * (vp - vm) * dyinv;
+
+#ifndef AMREX_USE_GPU
+        } else {
+            amrex::Error("ERROR: invalid coord_type in shock");
+#endif
+        }
 
-      Real p_jm1{};
-      {
-          Real lrhoInv = 1.0 / dat(i,j-1,k,URHO);
+        // we need to compute p in the full stencil
 
-          eos_state.rho = dat(i,j-1,k,URHO);
-          eos_state.T = dat(i,j-1,k,UTEMP);
-          eos_state.e = dat(i,j-1,k,UEINT) * lrhoInv;
-          for (int n = 0; n < NumSpec; n++) {
-              eos_state.xn[n] = dat(i,j-1,k,UFS+n) * lrhoInv;
-          }
+        Real p_ip1{};
+        {
+            Real lrhoInv = 1.0 / dat(i+1,j,k,URHO);
+
+            eos_state.rho = dat(i+1,j,k,URHO);
+            eos_state.T = dat(i+1,j,k,UTEMP);
+            eos_state.e = dat(i+1,j,k,UEINT) * lrhoInv;
+            for (int n = 0; n < NumSpec; n++) {
+                eos_state.xn[n] = dat(i+1,j,k,UFS+n) * lrhoInv;
+            }
 #if NAUX_NET > 0
-          for (int n = 0; n < NumAux; n++) {
-              eos_state.aux[n] = dat(i,j-1,k,UFX+n) * lrhoInv
-          }
+            for (int n = 0; n < NumAux; n++) {
+                eos_state.aux[n] = dat(i+1,j,k,UFX+n) * lrhoInv;
+            }
 #endif
 
-          if (eos_state.e <= 0.0_rt) {
-              eos(eos_input_rt, eos_state);
-          } else {
-              eos(eos_input_re, eos_state);
-          }
-          p_jm1 = eos_state.p;
-      }
+            if (eos_state.e <= 0.0_rt) {
+                eos(eos_input_rt, eos_state);
+            } else {
+                eos(eos_input_re, eos_state);
+            }
+            p_ip1 = eos_state.p;
+        }
+
+        Real p_im1{};
+        {
+            Real lrhoInv = 1.0 / dat(i-1,j,k,URHO);
+
+            eos_state.rho = dat(i-1,j,k,URHO);
+            eos_state.T = dat(i-1,j,k,UTEMP);
+            eos_state.e = dat(i-1,j,k,UEINT) * lrhoInv;
+            for (int n = 0; n < NumSpec; n++) {
+                eos_state.xn[n] = dat(i-1,j,k,UFS+n) * lrhoInv;
+            }
+#if NAUX_NET > 0
+            for (int n = 0; n < NumAux; n++) {
+                eos_state.aux[n] = dat(i-1,j,k,UFX+n) * lrhoInv;
+            }
+#endif
 
-      Real dP_y = 0.5_rt * (p_jp1 - p_jm1);
+            if (eos_state.e <= 0.0_rt) {
+                eos(eos_input_rt, eos_state);
+            } else {
+                eos(eos_input_re, eos_state);
+            }
+            p_im1 = eos_state.p;
+        }
+
+        Real dP_x = 0.5_rt * (p_ip1 - p_im1);
+
+        Real p_jp1{};
+        {
+            Real lrhoInv = 1.0 / dat(i,j+1,k,URHO);
+
+            eos_state.rho = dat(i,j+1,k,URHO);
+            eos_state.T = dat(i,j+1,k,UTEMP);
+            eos_state.e = dat(i,j+1,k,UEINT) * lrhoInv;
+            for (int n = 0; n < NumSpec; n++) {
+                eos_state.xn[n] = dat(i,j+1,k,UFS+n) * lrhoInv;
+            }
+#if NAUX_NET > 0
+            for (int n = 0; n < NumAux; n++) {
+                eos_state.aux[n] = dat(i,j+1,k,UFX+n) * lrhoInv;
+            }
+#endif
 
-      //Real gradPdx_over_P = std::sqrt(dP_x * dP_x + dP_y * dP_y + dP_z * dP_z) / dat(i,j,k,QPRES);
-      Real vel = std::sqrt(dat(i,j,k,UMX) * dat(i,j,k,UMX) +
-                           dat(i,j,k,UMY) * dat(i,j,k,UMY) +
-                           dat(i,j,k,UMZ) * dat(i,j,k,UMZ)) / dat(i,j,k,URHO);
+            if (eos_state.e <= 0.0_rt) {
+                eos(eos_input_rt, eos_state);
+            } else {
+                eos(eos_input_re, eos_state);
+            }
+            p_jp1 = eos_state.p;
+        }
+
+        Real p_jm1{};
+        {
+            Real lrhoInv = 1.0 / dat(i,j-1,k,URHO);
+
+            eos_state.rho = dat(i,j-1,k,URHO);
+            eos_state.T = dat(i,j-1,k,UTEMP);
+            eos_state.e = dat(i,j-1,k,UEINT) * lrhoInv;
+            for (int n = 0; n < NumSpec; n++) {
+                eos_state.xn[n] = dat(i,j-1,k,UFS+n) * lrhoInv;
+            }
+#if NAUX_NET > 0
+            for (int n = 0; n < NumAux; n++) {
+                eos_state.aux[n] = dat(i,j-1,k,UFX+n) * lrhoInv;
+            }
+#endif
 
-      Real gradPdx_over_P{0.0_rt};
-      if (vel != 0.0) {
-          gradPdx_over_P = std::abs(dP_x * dat(i,j,k,UMX) / dat(i,j,k,URHO) +
-                                    dP_y * dat(i,j,k,UMY) / dat(i,j,k,URHO)) / vel;
-      }
-      gradPdx_over_P /= p_zone;
+            if (eos_state.e <= 0.0_rt) {
+                eos(eos_input_rt, eos_state);
+            } else {
+                eos(eos_input_re, eos_state);
+            }
+            p_jm1 = eos_state.p;
+        }
+
+        Real dP_y = 0.5_rt * (p_jp1 - p_jm1);
+
+        //Real gradPdx_over_P = std::sqrt(dP_x * dP_x + dP_y * dP_y + dP_z * dP_z) / dat(i,j,k,QPRES);
+        Real du_x = std::min(up - um, 0.0);
+        Real dv_y = std::min(vp - vm, 0.0);
+
+        Real divu_mag = std::sqrt(du_x * du_x + dv_y * dv_y + 1.e-30);
+
+        Real gradPdx_over_P = std::abs(dP_x * du_x + dP_y * dv_y) / divu_mag;
+        gradPdx_over_P /= p_zone;
 
-      der(i,j,k,0) = gradPdx_over_P;
+        der(i,j,k,0) = gradPdx_over_P;
 
   });
 
diff --git a/Exec/science/subch_planar/Problem_Derives.H b/Exec/science/subch_planar/Problem_Derives.H
index 300db54eda..ae147a39c1 100644
--- a/Exec/science/subch_planar/Problem_Derives.H
+++ b/Exec/science/subch_planar/Problem_Derives.H
@@ -3,3 +3,6 @@
     //
     derive_lst.add("gradp_over_p", IndexType::TheCellType(), 1, ca_dergradpoverp, grow_box_by_one);
     derive_lst.addComponent("gradp_over_p", desc_lst, State_Type, URHO, NUM_STATE);
+
+    derive_lst.add("gradp_over_p1", IndexType::TheCellType(), 1, ca_dergradpoverp1, grow_box_by_one);
+    derive_lst.addComponent("gradp_over_p1", desc_lst, State_Type, URHO, NUM_STATE);