diff --git a/Exec/science/subch_planar/Problem_Derive.H b/Exec/science/subch_planar/Problem_Derive.H
index 2110a535a1..9173a634e9 100644
--- a/Exec/science/subch_planar/Problem_Derive.H
+++ b/Exec/science/subch_planar/Problem_Derive.H
@@ -8,3 +8,18 @@ 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*/);
+
+void ca_dergradpx
+  (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_dergradpy
+  (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_dergradrhooverrho
+  (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 6a67c0532f..2e97187e44 100644
--- a/Exec/science/subch_planar/Problem_Derive.cpp
+++ b/Exec/science/subch_planar/Problem_Derive.cpp
@@ -210,6 +210,7 @@ void ca_dergradpoverp(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*nco
 
 }
 
+
 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*/)
@@ -411,3 +412,306 @@ void ca_dergradpoverp1(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*nc
   });
 
 }
+
+
+void ca_dergradpx(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 / p . xhat
+
+    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
+
+    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)
+    {
+
+        // 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,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;
+
+
+        // we need to compute p in the full stencil
+
+        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
+
+            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
+
+            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);
+
+        der(i,j,k,0) = std::abs(dp_x) / p_zone;
+
+    });
+
+}
+
+
+void ca_dergradpy(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 / p . yhat
+
+    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
+
+    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)
+    {
+
+        // 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,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;
+
+
+
+        // we need to compute p in the full stencil
+
+        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
+
+            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_y = 0.5_rt * (p_jp1 - p_jm1);
+
+        der(i,j,k,0) = std::abs(dp_y) / p_zone;
+
+  });
+
+}
+
+
+void ca_dergradrhooverrho(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/,
+                          const FArrayBox& datfab, const Geometry& geomdata,
+                          Real /*time*/, const int* /*bcrec*/, int /*level*/)
+{
+
+    // compute grad rho / rho
+
+    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
+
+    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;
+
+        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 drho_x = 0.5_rt * (dat(i+1,j,k,URHO) - dat(i-1,j,k,URHO));
+        Real drho_y = 0.5_rt * (dat(i,j+1,k,URHO) - dat(i,j-1,k,URHO));
+
+        Real vel = std::sqrt(u0 * u0 + v0 * v0);
+
+        Real gradrhodx_over_rho{0.0_rt};
+        if (vel != 0.0) {
+            gradrhodx_over_rho = std::abs(drho_x * u0 + drho_y * v0) / vel;
+        }
+        gradrhodx_over_rho /= dat(i,j,k,URHO);
+
+        der(i,j,k,0) = gradrhodx_over_rho;
+
+    });
+
+}
diff --git a/Exec/science/subch_planar/Problem_Derives.H b/Exec/science/subch_planar/Problem_Derives.H
index ae147a39c1..748a2edafd 100644
--- a/Exec/science/subch_planar/Problem_Derives.H
+++ b/Exec/science/subch_planar/Problem_Derives.H
@@ -6,3 +6,12 @@
 
     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);
+
+    derive_lst.add("gradp_x", IndexType::TheCellType(), 1, ca_dergradpx, grow_box_by_one);
+    derive_lst.addComponent("gradp_x", desc_lst, State_Type, URHO, NUM_STATE);
+
+    derive_lst.add("gradp_y", IndexType::TheCellType(), 1, ca_dergradpy, grow_box_by_one);
+    derive_lst.addComponent("gradp_y", desc_lst, State_Type, URHO, NUM_STATE);
+
+    derive_lst.add("gradrho_over_rho", IndexType::TheCellType(), 1, ca_dergradrhooverrho, grow_box_by_one);
+    derive_lst.addComponent("gradrho_over_rho", desc_lst, State_Type, URHO, NUM_STATE);