diff --git a/Source/Particles/ERFPCEvolve.cpp b/Source/Particles/ERFPCEvolve.cpp
index 0aed01002..47a0efd27 100644
--- a/Source/Particles/ERFPCEvolve.cpp
+++ b/Source/Particles/ERFPCEvolve.cpp
@@ -231,10 +231,16 @@ void ERFPC::AdvectWithGravity (  int                                 a_lev,
             iv[1] += domain.smallEnd()[1];
             ParticleReal zlo, zhi;
             if (use_terrain) {
-                zlo = Real(0.25) * (zheight(iv[0],iv[1]  ,iv[2]  ) + zheight(iv[0]+1,iv[1]  ,iv[2]  ) +
-                                    zheight(iv[0],iv[1]+1,iv[2]  ) + zheight(iv[0]+1,iv[1]+1,iv[2]  ));
-                zhi = Real(0.25) * (zheight(iv[0],iv[1]  ,iv[2]+1) + zheight(iv[0]+1,iv[1]  ,iv[2]+1) +
-                                    zheight(iv[0],iv[1]+1,iv[2]+1) + zheight(iv[0]+1,iv[1]+1,iv[2]+1));
+                Real lx = (p.pos(0)-plo[0])*dxi[0] - static_cast<Real>(iv[0]-domain.smallEnd()[0]);
+                Real ly = (p.pos(1)-plo[1])*dxi[1] - static_cast<Real>(iv[1]-domain.smallEnd()[1]);
+                zlo =  zheight(iv[0]  ,iv[1]  ,iv[2]  ) * (1.0-lx) * (1.0-ly) +
+                       zheight(iv[0]+1,iv[1]  ,iv[2]  ) *      lx  * (1.0-ly) +
+                       zheight(iv[0]  ,iv[1]+1,iv[2]  ) * (1.0-lx) * ly +
+                       zheight(iv[0]+1,iv[1]+1,iv[2]  ) *      lx  * ly;
+                zhi =  zheight(iv[0]  ,iv[1]  ,iv[2]+1) * (1.0-lx) * (1.0-ly) +
+                       zheight(iv[0]+1,iv[1]  ,iv[2]+1) *      lx  * (1.0-ly) +
+                       zheight(iv[0]  ,iv[1]+1,iv[2]+1) * (1.0-lx) * ly +
+                       zheight(iv[0]+1,iv[1]+1,iv[2]+1) *      lx  * ly;
             } else {
                 zlo =  iv[2]    * dx[2];
                 zhi = (iv[2]+1) * dx[2];