diff --git a/src/prob/incflo_prob_I.H b/src/prob/incflo_prob_I.H
index b4b65712..0cb105c9 100644
--- a/src/prob/incflo_prob_I.H
+++ b/src/prob/incflo_prob_I.H
@@ -165,8 +165,9 @@
                     amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& problo,
                     amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& probhi) const;
 
-    void column_collapse_granular (amrex::Box const& vbx,
+    void column_collapse_granular (amrex::Box const& vbx, amrex::Box const& nbx,
                                    amrex::Array4<amrex::Real> const& density,
+                                   amrex::Array4<amrex::Real> const& pressure,
                                    amrex::Box const& domain,
                                    amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& dx,
                                    amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& problo,
diff --git a/src/prob/prob_init_fluid.cpp b/src/prob/prob_init_fluid.cpp
index a22eb55b..41a36eba 100644
--- a/src/prob/prob_init_fluid.cpp
+++ b/src/prob/prob_init_fluid.cpp
@@ -26,6 +26,7 @@ void incflo::prob_init_fluid (int lev)
     {
         const Box& vbx = mfi.validbox();
         const Box& gbx = mfi.fabbox();
+        const Box& nbx = mfi.tilebox(IntVect::TheNodeVector());
         if (0 == m_probtype || 114 == m_probtype )
         { }
         else if (1 == m_probtype)
@@ -187,7 +188,9 @@ void incflo::prob_init_fluid (int lev)
         }
         else if (531 == m_probtype)
         {
-            column_collapse_granular(vbx, ld.density.array(mfi),
+            column_collapse_granular(vbx, nbx,
+                                     ld.density.array(mfi),
+                                     ld.p_nd.array(mfi),
                                      domain, dx, problo, probhi);
         }
         else
@@ -1126,8 +1129,9 @@ void incflo::init_burggraf (Box const& vbx, Box const& /*gbx*/,
     });
 }
 
-void incflo::column_collapse_granular (Box const& vbx,
+void incflo::column_collapse_granular (Box const& vbx, Box const& nbx,
                                  Array4<Real> const& density,
+                                 Array4<Real> const& pressure,
                                  Box const& /*domain*/,
                                  GpuArray<Real, AMREX_SPACEDIM> const& dx,
                                  GpuArray<Real, AMREX_SPACEDIM> const& problo,
@@ -1135,14 +1139,6 @@ void incflo::column_collapse_granular (Box const& vbx,
 {
     // Ensure it is set to two_fluid
     if (!m_two_fluid) amrex::Abort("probtype 531 involves two fluids");
-    // Use m_gp0 through m_delp instead of gravity
-    // This will ensure p_nd has both static and dynamic contributions
-    if (m_gravity[0]*m_gravity[0] + m_gravity[1]*m_gravity[1]
-        + m_gravity[2]*m_gravity[2] > Real(0.0)) 
-        amrex::Abort("Use delp instead of gravity for probtype 531");
-
-    if (m_gp0[0]*m_gp0[0] + m_gp0[1]*m_gp0[1] + m_gp0[2]*m_gp0[2] 
-        == Real(0.0)) amrex::Abort("Use delp for probtype 531");
 
     Vector<Real> granlen_vec{AMREX_D_DECL(Real(0.),Real(0.),Real(0.))};
     ParmParse pp;
@@ -1151,16 +1147,17 @@ void incflo::column_collapse_granular (Box const& vbx,
         granlen_vec[0],granlen_vec[1],granlen_vec[2])};
     granLen[0] += problo[0];
     granLen[1] += problo[1];
-    granLen[2] += problo[2];
-    if (granLen[0] > probhi[0]) 
+    if (granLen[0] > probhi[0])
         amrex::Abort("Granular length along x is larger than setup");
-    if (granLen[1] > probhi[1]) 
+    if (granLen[1] > probhi[1])
         amrex::Abort("Granular length along y is larger than setup");
-    if (granLen[2] > probhi[2]) 
+#if (AMREX_SPACEDIM==3)
+    granLen[2] += problo[2];
+    if (granLen[2] > probhi[2])
         amrex::Abort("Granular length along z is larger than setup");
-    if (m_initial_iterations < 10) 
-        amrex::Abort("Set large value for m_initial_iterations as pressure is NOT initialized");
+#endif
     Real rho_1 = m_ro_0; Real rho_2 = m_ro_0_second;
+    // Density
     amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
     {
         Real x = problo[0] + (Real(i)+Real(0.5))*dx[0];
@@ -1168,14 +1165,53 @@ void incflo::column_collapse_granular (Box const& vbx,
 #if (AMREX_SPACEDIM == 3)
         Real z = problo[2] + (Real(k)+Real(0.5))*dx[2];
 #endif
-        if ( x <= granLen[0] and y <= granLen[1]
+        if ( (x <= granLen[0]) and (y <= granLen[1])
 #if (AMREX_SPACEDIM == 3)
-            and z <= granLen[3]
+            and (z <= granLen[2])
 #endif
            ) {
             density(i,j,k) = rho_2;
         } else {
             density(i,j,k) = rho_1;
         }
-    });  
+    });
+    // Pressure
+    GpuArray<Real,AMREX_SPACEDIM> grav{0.,0.,0.};
+    if (m_gravity[0]*m_gravity[0] + m_gravity[1]*m_gravity[1]+
+        m_gravity[2]*m_gravity[2] > Real(0.)) {
+        grav[0] = std::abs(m_gravity[0]);
+        grav[1] = std::abs(m_gravity[1]);
+        grav[2] = std::abs(m_gravity[2]);
+    } else if (m_gp0[0]*m_gp0[0] + m_gp0[1]*m_gp0[1]+
+        m_gp0[2]*m_gp0[2] > Real(0.)) {
+        grav[0] = std::abs(m_gp0[0]);
+        grav[1] = std::abs(m_gp0[1]);
+        grav[2] = std::abs(m_gp0[2]);
+    }
+    else {
+        amrex::Abort("Body force required in probtype 531");
+    }
+
+    amrex::ParallelFor(nbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+    {
+        // Taking pressure based on the last dimension
+        Real x = problo[0] + Real(i)*dx[0];
+        Real y = problo[1] + Real(j)*dx[1];
+        Real p_h = y;
+#if (AMREX_SPACEDIM == 3)
+        Real z = problo[2] + Real(k)*dx[2];
+        p_h = z;
+#endif
+        int p_dim = AMREX_SPACEDIM - 1;
+        if ( (x <= granLen[0]) and (y <= granLen[1])
+#if (AMREX_SPACEDIM == 3)
+            and (z <= granLen[2])
+#endif
+           ) {
+            pressure(i,j,k) = rho_1*grav[p_dim]*(probhi[p_dim]-granLen[p_dim])
+                              + rho_2*grav[p_dim]*(granLen[p_dim]-p_h);
+        } else {
+            pressure(i,j,k) = rho_1*grav[p_dim]*(probhi[p_dim]-p_h);
+        }
+    });
 }