Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 24759)
+++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 24760)
@@ -1423,4 +1423,5 @@
 	IssmDouble* Zcum=NULL;
 	IssmDouble* dzMin2=NULL;
+	IssmDouble* dzMax2=NULL;
 	IssmDouble zY2=zY;
 	bool lastCellFlag = false;
@@ -1689,4 +1690,5 @@
 		celldelete(&W,n,D,D_size);
 		celldelete(&d,n,D,D_size);
+		celldelete(&dz,n,D,D_size);
 		celldelete(&T,n,D,D_size);
 		celldelete(&a,n,D,D_size);
@@ -1710,87 +1712,60 @@
 	Zcum=xNew<IssmDouble>(n);
 	dzMin2=xNew<IssmDouble>(n);
-
+	dzMax2=xNew<IssmDouble>(n);
+
+	X=0;
 	Zcum[0]=dz[0]; // Compute a cumulative depth vector
-
-	for (int i=1;i<n;i++){
-		Zcum[i]=Zcum[i-1]+dz[i];
-	}
-
 	for (int i=0;i<n;i++){
-		if (Zcum[i]<=zTop+Dtol || i<1){
+		if (i==0){
 			dzMin2[i]=dzMin;
 		}
 		else{
-			dzMin2[i]=zY2*dzMin2[i-1];
-		}
-	}
-
-	// check if depth is too small
-	X = 0;
-	for(int i=n-1;i>=0;i--){
-		if(dz[i]<dzMin2[i]-Dtol){
-			X=i;
-			break;
-		}
-	}
-
-	//Last cell has to be treated separately if has to be merged (no underlying cell so need to merge with overlying cell)
-	if(X==n-1){
-		lastCellFlag = true;
-		X = X-1;
-	}
-	else{
-		lastCellFlag = false;
-	}
-
-	for (int i = 0; i<=X;i++){
-		if (dz[i] < dzMin2[i]-Dtol && d[i+1]<dIce-Dtol){  // merge top cells, don't merge with ice
-			// _printf_("dz > dzMin * 2')
+			Zcum[i]=Zcum[i-1]+dz[i];
+			if (Zcum[i]<=zTop+Dtol){
+				dzMin2[i]=dzMin;
+				X=i;
+			}
+			else{
+				dzMin2[i]=zY2*dzMin2[i-1];
+			}
+		}
+	}
+
+	// Check to see if any cells are too small and need to be merged
+	for (int i=0; i<n; i++){
+		if ( (i<=X && dz[i]<dzMin-Dtol) || (i>X && dz[i]<dzMin2[i]-Dtol) ) {
+
+			if (i==n-1){
+				X2=i;
+				//find closest cell to merge with
+				for(int j=n-2;j>=0;j--){
+					if(m[j]!=Delflag){
+						X1=j;
+						break;
+					}
+				}
+			}
+			else{
+				X1=i+1;
+				X2=i;
+			}
+
 			// adjust variables as a linearly weighted function of mass
-			IssmDouble m_new = m[i] + m[i+1];
-			T[i+1] = (T[i]*m[i] + T[i+1]*m[i+1]) / m_new;
-			a[i+1] = (a[i]*m[i] + a[i+1]*m[i+1]) / m_new;
-			re[i+1] = (re[i]*m[i] + re[i+1]*m[i+1]) / m_new;
-			gdn[i+1] = (gdn[i]*m[i] + gdn[i+1]*m[i+1]) / m_new;
-			gsp[i+1] = (gsp[i]*m[i] + gsp[i+1]*m[i+1]) / m_new;
+			IssmDouble m_new = m[X2] + m[X1];
+			T[X1] = (T[X2]*m[X2] + T[X1]*m[X1]) / m_new;
+			a[X1] = (a[X2]*m[X2] + a[X1]*m[X1]) / m_new;
+			re[X1] = (re[X2]*m[X2] + re[X1]*m[X1]) / m_new;
+			gdn[X1] = (gdn[X2]*m[X2] + gdn[X1]*m[X1]) / m_new;
+			gsp[X1] = (gsp[X2]*m[X2] + gsp[X1]*m[X1]) / m_new;
 
 			// merge with underlying grid cell and delete old cell
-			dz[i+1] = dz[i] + dz[i+1];                 // combine cell depths
-			d[i+1] = m_new / dz[i+1];                   // combine top densities
-			W[i+1] = W[i+1] + W[i];                     // combine liquid water
-			m[i+1] = m_new;                             // combine top masses
+			dz[X1] = dz[X2] + dz[X1];                 // combine cell depths
+			d[X1] = m_new / dz[X1];                   // combine top densities
+			W[X1] = W[X1] + W[X2];                     // combine liquid water
+			m[X1] = m_new;                             // combine top masses
 
 			// set cell to -99999 for deletion
-			m[i] = Delflag;
-		}
-	}
-
-	//If last cell has to be merged
-	if(lastCellFlag){
-      //find closest cell to merge with
-		for(int i=n-2;i>=0;i--){
-			if(m[i]!=Delflag){
-				X2=i;
-				X1=n-1;
-				break;
-			}
-		}
-
-		// adjust variables as a linearly weighted function of mass
-		IssmDouble m_new = m[X2] + m[X1];
-		T[X1] = (T[X2]*m[X2] + T[X1]*m[X1]) / m_new;
-		a[X1] = (a[X2]*m[X2] + a[X1]*m[X1]) / m_new;
-		re[X1] = (re[X2]*m[X2] + re[X1]*m[X1]) / m_new;
-		gdn[X1] = (gdn[X2]*m[X2] + gdn[X1]*m[X1]) / m_new;
-		gsp[X1] = (gsp[X2]*m[X2] + gsp[X1]*m[X1]) / m_new;
-
-		// merge with underlying grid cell and delete old cell
-  		dz[X1] = dz[X2] + dz[X1];                 // combine cell depths
-		d[X1] = m_new / dz[X1];                   // combine top densities
- 		W[X1] = W[X1] + W[X2];                     // combine liquid water
-		m[X1] = m_new;                             // combine top masses
-
-		// set cell to -99999 for deletion
-		m[X2] = Delflag;
+			m[X2] = Delflag;
+		}
 	}
 
@@ -1814,17 +1789,25 @@
 	xDelete<int>(D);
 
-	// check if any of the top 10 cell depths are too large
+	// check if any of the cell depths are too large
 	X=0;
-	for(int i=min(9,n-1);i>=0;i--){
-		if(dz[i]> 2.0*dzMin+Dtol){
-			X=i;
-			break;
-		}
-	}
-
-	int j=0;
-	while(j<=X){
-		if (dz[j] > dzMin*2.0+Dtol){
-
+	Zcum[0]=dz[0]; // Compute a cumulative depth vector
+	for (int i=0;i<n;i++){
+		if (i==0){
+			dzMax2[i]=dzMin*2;
+		}
+		else{
+			Zcum[i]=Zcum[i-1]+dz[i];
+			if (Zcum[i]<=zTop+Dtol){
+				dzMax2[i]=dzMin*2;
+				X=i;
+			}
+			else{
+				dzMax2[i]=max(zY2*dzMin2[i-1],dzMin*2);
+			}
+		}
+	}
+
+	for (int j=n-1;j>=0;j--){
+		if ((j<X && dz[j] > dzMax2[j]+Dtol) || (dz[j] > dzMax2[j]*zY2+Dtol)){
 			// _printf_("dz > dzMin * 2");
 			// split in two
@@ -1841,7 +1824,5 @@
 			cellsplit(&gsp, n, j,1.0);
 			n++;
-			X=X+1;
-		}
-		else j++;
+		}
 	}
 
Index: /issm/trunk-jpl/src/m/classes/SMBgemb.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/SMBgemb.m	(revision 24759)
+++ /issm/trunk-jpl/src/m/classes/SMBgemb.m	(revision 24760)
@@ -344,5 +344,5 @@
 			fielddisplay(self,'Aini','Initial albedo when restart [-]');
 			fielddisplay(self,'Tini','Initial snow temperature when restart [K]');
-			fielddisplay(self,'Sizeini','Initial number of layers when restart [K]');
+			fielddisplay(self,'Sizeini','Initial number of layers when restart [-]');
 
 
Index: /issm/trunk-jpl/src/m/classes/SMBgemb.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/SMBgemb.py	(revision 24759)
+++ /issm/trunk-jpl/src/m/classes/SMBgemb.py	(revision 24760)
@@ -158,5 +158,5 @@
         string = "%s\n%s" % (string, fielddisplay(self, 'Aini', 'Initial albedo when restart [-]'))
         string = "%s\n%s" % (string, fielddisplay(self, 'Tini', 'Initial snow temperature when restart [K]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'Sizeini', 'Initial number of layers when restart [K]'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'Sizeini', 'Initial number of layers when restart [-]'))
 
         #additional albedo parameters:
Index: /issm/trunk-jpl/test/NightlyRun/test243.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test243.m	(revision 24759)
+++ /issm/trunk-jpl/test/NightlyRun/test243.m	(revision 24760)
@@ -9,5 +9,5 @@
 % Use of Gemb method for SMB computation
 md.smb = SMBgemb(md.mesh,md.geometry);
-md.smb.dsnowIdx = 0;
+md.smb.dsnowIdx = 1;
 
 %load hourly surface forcing date from 1979 to 2009:
@@ -46,9 +46,12 @@
 md=solve(md,'Transient');
 
-nlayers=size(md.results.TransientSolution(end).SmbT,2);
+nlayers=size(md.results.TransientSolution(1).SmbT,2);
+for i=2:length(md.results.TransientSolution)
+   nlayers=min(size(md.results.TransientSolution(i).SmbT,2), nlayers);
+end
 
 %Fields and tolerances to track changes
 field_names      ={'Layers','SmbDz','SmbT' ,'SmbD' ,'SmbRe','SmbGdn','SmbGsp','SmbA' ,'SmbEC','SmbMassBalance','SmbMAdd','SmbDzAdd','SmbFAC','SmbMeanSHF','SmbMeanLHF','SmbMeanULW','SmbNetLW','SmbNetSW'};
-field_tolerances ={1e-12,1e-12,1e-12,1e-11,1e-11,1e-11,1e-11,1e-12,1e-11,1e-12,1e-12,1e-12,1e-11,2e-11,2e-11,1e-11,9e-10,2e-11};
+field_tolerances ={1e-12,2e-12,1e-12,1e-11,1e-11,2e-11,1e-11,1e-12,1e-11,1e-12,1e-12,1e-12,1e-11,2e-11,2e-11,1e-11,9e-10,2e-11};
 
 field_values={...
Index: /issm/trunk-jpl/test/NightlyRun/test243.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test243.py	(revision 24759)
+++ /issm/trunk-jpl/test/NightlyRun/test243.py	(revision 24760)
@@ -21,5 +21,5 @@
 md.smb = SMBgemb()
 md.smb.setdefaultparameters(md.mesh, md.geometry)
-md.smb.dsnowIdx = 0
+md.smb.dsnowIdx = 1
 
 #load hourly surface forcing date from 1979 to 2009:
@@ -60,9 +60,11 @@
 md = solve(md, 'Transient')
 
-nlayers=np.size(md.results.TransientSolution[-1].SmbT,1)
+nlayers=np.size(md.results.TransientSolution[0].SmbT,1)
+for i in range(1, np.size(md.results.TransientSolution,0)):
+    nlayers=np.minimum(np.size(md.results.TransientSolution[i].SmbT,1),nlayers)
 
 #Fields and tolerances to track changes
 field_names = ['Layers','SmbDz', 'SmbT', 'SmbD', 'SmbRe', 'SmbGdn', 'SmbGsp', 'SmbA', 'SmbEC', 'SmbMassBalance', 'SmbMAdd', 'SmbDzAdd', 'SmbFAC', 'SmbMeanSHF', 'SmbMeanLHF', 'SmbMeanULW', 'SmbNetLW', 'SmbNetSW']
-field_tolerances = [1e-12,1e-12,1e-12,1e-11,1e-11,1e-11,1e-11,1e-12,1e-11,1e-12,1e-12,1e-12,1e-11,2e-11,2e-11,1e-11,9e-10,2e-11] 
+field_tolerances = [1e-12,2e-12,1e-12,1e-11,1e-11,2e-11,1e-11,1e-12,1e-11,1e-12,1e-12,1e-12,1e-11,2e-11,2e-11,1e-11,9e-10,2e-11] 
 #shape is different in python solution (fixed using reshape) which can cause test failure:
 field_values = [nlayers,md.results.TransientSolution[-1].SmbDz[0, 0:nlayers].reshape(1, -1),
Index: /issm/trunk-jpl/test/NightlyRun/test252.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test252.m	(revision 24759)
+++ /issm/trunk-jpl/test/NightlyRun/test252.m	(revision 24760)
@@ -9,5 +9,5 @@
 % Use of Gemb method for SMB computation
 md.smb = SMBgemb(md.mesh,md.geometry);
-md.smb.dsnowIdx = 0;
+md.smb.dsnowIdx = 4;
 
 %load hourly surface forcing date from 1979 to 2009:
@@ -54,5 +54,8 @@
 md=solve(md,'Transient');
 
-nlayers=size(md.results.TransientSolution(end).SmbT,2);
+nlayers=size(md.results.TransientSolution(1).SmbT,2);
+for i=2:length(md.results.TransientSolution)
+   nlayers=min(size(md.results.TransientSolution(i).SmbT,2), nlayers);
+end
 
 %Fields and tolerances to track changes
@@ -62,6 +65,6 @@
 	'SmbDz4','SmbT4' ,'SmbD4' ,'SmbRe4','SmbGdn4','SmbGsp4','SmbA4' ,'SmbEC4','SmbMassBalance4','SmbMAdd4','SmbDzAdd4','SmbFAC4'};
 field_tolerances ={1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,...
-                   1e-12,1e-12,1e-11,1e-10,1e-11,1e-11,1e-12,4e-12,1e-12,1e-12,1e-12,1e-11,...
-                   1e-12,1e-12,2e-12,2e-11,1e-11,1e-11,1e-12,1e-11,1e-11,1e-12,1e-12,1e-11,...
+                   1e-12,1e-12,1e-11,1e-10,2e-11,1e-11,1e-12,4e-12,1e-12,1e-12,1e-12,1e-11,...
+                   1e-12,1e-12,2e-12,2e-11,4e-11,1e-11,1e-12,1e-11,1e-11,1e-12,1e-12,1e-11,...
                    1e-11,1e-11,4e-11,4e-11,1e-12,3e-11,1e-12,1e-12,1e-10,1e-12,1e-12,2e-11};
 
Index: /issm/trunk-jpl/test/NightlyRun/test252.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test252.py	(revision 24759)
+++ /issm/trunk-jpl/test/NightlyRun/test252.py	(revision 24760)
@@ -21,5 +21,5 @@
 md.smb = SMBgemb()
 md.smb.setdefaultparameters(md.mesh, md.geometry)
-md.smb.dsnowIdx = 0
+md.smb.dsnowIdx = 4
 
 #load hourly surface forcing date from 1979 to 2009:
@@ -69,5 +69,7 @@
 md = solve(md, 'Transient')
 
-nlayers=np.size(md.results.TransientSolution[-1].SmbT,1)
+nlayers=np.size(md.results.TransientSolution[0].SmbT,1)
+for i in range(1, np.size(md.results.TransientSolution,0)):
+    nlayers=np.minimum(np.size(md.results.TransientSolution[i].SmbT,1),nlayers)
 
 #Fields and tolerances to track changes
@@ -78,6 +80,6 @@
                'SmbDz4', 'SmbT4', 'SmbD4', 'SmbRe4', 'SmbGdn4', 'SmbGsp4', 'SmbA4', 'SmbEC4', 'SmbMassBalance4', 'SmbMAdd4', 'SmbDzAdd4', 'SmbFAC4']
 field_tolerances = [1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,
-                   1e-12,1e-12,1e-11,1e-10,1e-11,1e-11,1e-12,4e-12,1e-12,1e-12,1e-12,1e-11,
-                   1e-12,1e-12,2e-12,2e-11,1e-11,1e-11,1e-12,1e-11,1e-11,1e-12,1e-12,1e-11,
+                   1e-12,1e-12,1e-11,1e-10,2e-11,1e-11,1e-12,4e-12,1e-12,1e-12,1e-12,1e-11,
+                   1e-12,1e-12,2e-12,2e-11,4e-11,1e-11,1e-12,1e-11,1e-11,1e-12,1e-12,1e-11,
                    1e-11,1e-11,4e-11,4e-11,1e-12,3e-11,1e-12,1e-12,1e-10,1e-12,1e-12,2e-11]
 
Index: /issm/trunk-jpl/test/NightlyRun/test253.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test253.m	(revision 24759)
+++ /issm/trunk-jpl/test/NightlyRun/test253.m	(revision 24760)
@@ -58,4 +58,7 @@
 
 nlayers=size(md.results.TransientSolution(1).SmbT,2);
+for i=2:length(md.results.TransientSolution)
+	nlayers=min(size(md.results.TransientSolution(i).SmbT,2), nlayers);
+end
 
 %Fields and tolerances to track changes
Index: /issm/trunk-jpl/test/NightlyRun/test253.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test253.py	(revision 24759)
+++ /issm/trunk-jpl/test/NightlyRun/test253.py	(revision 24760)
@@ -73,4 +73,6 @@
 
 nlayers=np.size(md.results.TransientSolution[0].SmbT,1)
+for i in range(1, np.size(md.results.TransientSolution,0)):
+    nlayers=np.minimum(np.size(md.results.TransientSolution[i].SmbT,1),nlayers)
 
 #Fields and tolerances to track changes