diff --git a/src/simulation/Simulation.cpp b/src/simulation/Simulation.cpp
index 2588ad7e7..87603f2a9 100644
--- a/src/simulation/Simulation.cpp
+++ b/src/simulation/Simulation.cpp
@@ -88,6 +88,7 @@ int Simulation::Load(int fullX, int fullY, GameSave * save)
 		}
 	}
 	parts_lastActiveIndex = NPART-1;
+	force_stacking_check = 1;
 	for(int i = 0; i < save->signs.size() && signs.size() < MAXSIGNS; i++)
 	{
 		sign tempSign = save->signs[i];
@@ -2914,6 +2915,7 @@ void Simulation::update_particles_i(int start, int inc)
 	int lighting_ok=1;
 	unsigned int elem_properties;
 	float pGravX, pGravY, pGravD;
+	int excessive_stacking_found = 0;
 
 	if (lighting_recreate>0)
     {
@@ -2939,6 +2941,66 @@ void Simulation::update_particles_i(int start, int inc)
 	//if (sys_pause&&!framerender)//do nothing if paused
 	//	return;
 
+    if (force_stacking_check || (rand()%10)==0)
+    {
+    	force_stacking_check = 0;
+    	excessive_stacking_found = 0;
+    	for (y=0; y<YRES; y++)
+    	{
+    		for (x=0; x<XRES; x++)
+    		{
+    			// Use a threshold, since some particle stacking can be normal (e.g. BIZR + FILT)
+    			// Setting pmap_count[y][x] > NPART means BHOL will form in that spot
+    			if (pmap_count[y][x]>5)
+    			{
+    				if (bmap[y/CELL][x/CELL]==WL_EHOLE)
+    				{
+    					// Allow more stacking in E-hole
+    					if (pmap_count[y][x]>1500)
+    					{
+    						pmap_count[y][x] = pmap_count[y][x] + NPART;
+    						excessive_stacking_found = 1;
+    					}
+    				}
+    				else if (pmap_count[y][x]>1500 || (rand()%1600)<=(pmap_count[y][x]+100))
+    				{
+    					pmap_count[y][x] = pmap_count[y][x] + NPART;
+    					excessive_stacking_found = 1;
+    				}
+    			}
+    		}
+    	}
+    	if (excessive_stacking_found)
+    	{
+    		for (i=0; i<=parts_lastActiveIndex; i++)
+    		{
+    			if (parts[i].type)
+    			{
+    				t = parts[i].type;
+    				x = (int)(parts[i].x+0.5f);
+    				y = (int)(parts[i].y+0.5f);
+    				if (x>=0 && y>=0 && x<XRES && y<YRES && !(elements[t].Properties&TYPE_ENERGY))
+    				{
+    					if (pmap_count[y][x]>=NPART)
+    					{
+    						if (pmap_count[y][x]>NPART)
+    						{
+    							create_part(i, x, y, PT_NBHL);
+    							parts[i].temp = MAX_TEMP;
+    							parts[i].tmp = pmap_count[y][x]-NPART;//strength of grav field
+    							pmap_count[y][x] = NPART;
+    						}
+    						else
+    						{
+    							kill_part(i);
+    						}
+    					}
+    				}
+    			}
+    		}
+    	}
+    }
+
 	//wire!
 	if(elementCount[PT_WIRE] > 0)
 	{
@@ -4073,6 +4135,7 @@ void Simulation::update_particles()//doesn't update the particles themselves, bu
 	}
 
 	memset(pmap, 0, sizeof(pmap));
+	memset(pmap_count, 0, sizeof(pmap_count));
 	memset(photons, 0, sizeof(photons));
 	NUM_PARTS = 0;
 	for (i=0; i<=parts_lastActiveIndex; i++)//the particle loop that resets the pmap/photon maps every frame, to update them.
@@ -4087,7 +4150,10 @@ void Simulation::update_particles()//doesn't update the particles themselves, bu
 				if (elements[t].Properties & TYPE_ENERGY)
 					photons[y][x] = t|(i<<8);
 				else
+				{
 					pmap[y][x] = t|(i<<8);
+					pmap_count[y][x]++;
+				}
 			}
 			lastPartUsed = i;
 			NUM_PARTS ++;
diff --git a/src/simulation/Simulation.h b/src/simulation/Simulation.h
index b4b2c9346..8e25daad0 100644
--- a/src/simulation/Simulation.h
+++ b/src/simulation/Simulation.h
@@ -71,6 +71,7 @@ public:
 	int NUM_PARTS;
 	int elementCount[PT_NUM];
 	int ISWIRE;
+	int force_stacking_check;
 	//Gol sim
 	int CGOL;
 	int ISGOL;
@@ -96,6 +97,7 @@ public:
 	Particle parts[NPART];
 	int pmap[YRES][XRES];
 	int photons[YRES][XRES];
+	int pmap_count[YRES][XRES];
 	//
 	int gravityMode;
 	//int airMode;
diff --git a/src/simulation/elements/NBHL.cpp b/src/simulation/elements/NBHL.cpp
index 88d6ec22c..feb42ef4a 100644
--- a/src/simulation/elements/NBHL.cpp
+++ b/src/simulation/elements/NBHL.cpp
@@ -49,9 +49,12 @@ Element_NBHL::Element_NBHL()
 //#TPT-Directive ElementHeader Element_NBHL static int update(UPDATE_FUNC_ARGS)
 int Element_NBHL::update(UPDATE_FUNC_ARGS)
  {
-	sim->gravmap[(y/CELL)*(XRES/CELL)+(x/CELL)] += 0.1f;
+	if (parts[i].tmp)
+		sim->gravmap[(y/CELL)*(XRES/CELL)+(x/CELL)] += restrict_flt(0.001f*parts[i].tmp, 0.1f, 51.2f);
+	else
+		sim->gravmap[(y/CELL)*(XRES/CELL)+(x/CELL)] += 0.1f;
 	return 0;
 }
 
 
-Element_NBHL::~Element_NBHL() {}
\ No newline at end of file
+Element_NBHL::~Element_NBHL() {}