Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 27226)
+++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 27227)
@@ -1828,4 +1828,110 @@
 	xDelete<IssmDouble>(F);
 
+	//Manage the layering to match the user defined requirements
+	managelayers(&mAdd, &dz_add, &addE, &m, &EI, &EW, &T, &d, &dz, &W, &a, &adiff, &re, &gdn, &gsp, &n, dzMin, zMax, zMin, zTop, zY);
+
+	//// CHECK FOR MASS AND ENERGY CONSERVATION
+
+	// calculate final mass [kg] and energy [J]
+	sumER = Rsum * (LF + CtoK * CI);
+	for(int i=0;i<n;i++)EI[i] = m[i] * T[i] * CI;
+	for(int i=0;i<n;i++)EW[i] = W[i] * (LF + CtoK * CI);
+
+	mSum1 = cellsum(W,n) + cellsum(m,n) + Rsum;
+	sumE1 = cellsum(EI,n) + cellsum(EW,n);
+
+	/*checks: */
+	for(int i=0;i<n;i++) if (W[i]<0.0-Wtol) _error_("negative pore water generated in melt equations\n");
+
+	/*only in forward mode! avoid round in AD mode as it is not differentiable: */
+	#ifndef _HAVE_AD_
+	dm = round((mSum0 - mSum1 + mAdd)*100.0)/100.0;
+	dE = round(sumE0 - sumE1 - sumER +  addE - surplusE);
+	if (dm !=0  || dE !=0) _error_("mass or energy are not conserved in melt equations\n"
+			<< "dm: " << dm << " dE: " << dE << "\n");
+	#endif
+
+	
+	/*Free resources:*/
+	if(m)xDelete<IssmDouble>(m);
+	if(EI)xDelete<IssmDouble>(EI);
+	if(EW)xDelete<IssmDouble>(EW);
+	if(maxF)xDelete<IssmDouble>(maxF);
+	if(dW)xDelete<IssmDouble>(dW);
+	if(exsW)xDelete<IssmDouble>(exsW);
+	if(exsT)xDelete<IssmDouble>(exsT);
+	if(surpT)xDelete<IssmDouble>(surpT);
+	if(surpE)xDelete<IssmDouble>(surpE);
+	if(flxDn)xDelete<IssmDouble>(flxDn);
+	if(D)xDelete<int>(D);
+	if(M)xDelete<IssmDouble>(M);
+
+	/*Assign output pointers:*/
+	*pMs=Msurf;
+	*pM=sumM;
+	*pR=Rsum;
+	*pF=Fsum;
+	*pmAdd=mAdd;
+	*pdz_add=dz_add;
+
+	*pT=T;
+	*pd=d;
+	*pdz=dz;
+	*pW=W;
+	*pa=a;
+	*padiff=adiff;
+	*pre=re;
+	*pgdn=gdn;
+	*pgsp=gsp;
+	*pn=n;
+
+} /*}}}*/ 
+void managelayers(IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble* paddE, IssmDouble** pm, IssmDouble** pEI, IssmDouble** pEW, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** padiff, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, IssmDouble zY){ /*{{{*/
+
+	/*intermediary:*/
+	IssmDouble* Zcum=NULL;
+	IssmDouble* dzMin2=NULL;
+	IssmDouble* dzMax2=NULL;
+	int*        D=NULL;
+
+	IssmDouble zY2=zY;
+	IssmDouble X=0.0;
+	int X1=0;
+	int X2=0;
+	int D_size = 0;
+
+	IssmDouble Ztot=0.0;
+	IssmDouble T_bot=0.0;
+	IssmDouble m_bot=0.0;
+	IssmDouble dz_bot=0.0;
+	IssmDouble d_bot=0.0;
+	IssmDouble W_bot=0.0;
+	IssmDouble a_bot=0.0;
+	IssmDouble adiff_bot=0.0;
+	IssmDouble re_bot=0.0;
+	IssmDouble gdn_bot=0.0;
+	IssmDouble gsp_bot=0.0;
+	IssmDouble EI_bot=0.0;
+	IssmDouble EW_bot=0.0;
+	bool       top=false;
+
+	/*outputs:*/
+	IssmDouble  mAdd = 0.0;
+	IssmDouble  addE = 0.0;
+	IssmDouble  dz_add = 0.0;
+	IssmDouble* T=*pT;
+	IssmDouble* d=*pd;
+	IssmDouble* dz=*pdz;
+	IssmDouble* W=*pW;
+	IssmDouble* a=*pa;
+	IssmDouble* adiff=*padiff;
+	IssmDouble* re=*pre;
+	IssmDouble* gdn=*pgdn;
+	IssmDouble* gsp=*pgsp;
+	IssmDouble* m=*pm;
+	IssmDouble* EI=*pEI;
+	IssmDouble* EW=*pEW;
+	int         n=*pn;
+
 	//Merging of cells as they are burried under snow.
 	Zcum=xNew<IssmDouble>(n);
@@ -2019,38 +2125,5 @@
 	}
 
-	//// CHECK FOR MASS AND ENERGY CONSERVATION
-
-	// calculate final mass [kg] and energy [J]
-	sumER = Rsum * (LF + CtoK * CI);
-	for(int i=0;i<n;i++)EI[i] = m[i] * T[i] * CI;
-	for(int i=0;i<n;i++)EW[i] = W[i] * (LF + CtoK * CI);
-
-	mSum1 = cellsum(W,n) + cellsum(m,n) + Rsum;
-	sumE1 = cellsum(EI,n) + cellsum(EW,n);
-
-	/*checks: */
-	for(int i=0;i<n;i++) if (W[i]<0.0-Wtol) _error_("negative pore water generated in melt equations\n");
-
-	/*only in forward mode! avoid round in AD mode as it is not differentiable: */
-	#ifndef _HAVE_AD_
-	dm = round((mSum0 - mSum1 + mAdd)*100.0)/100.0;
-	dE = round(sumE0 - sumE1 - sumER +  addE - surplusE);
-	if (dm !=0  || dE !=0) _error_("mass or energy are not conserved in melt equations\n"
-			<< "dm: " << dm << " dE: " << dE << "\n");
-	#endif
-
 	/*Free resources:*/
-	if(m)xDelete<IssmDouble>(m);
-	if(EI)xDelete<IssmDouble>(EI);
-	if(EW)xDelete<IssmDouble>(EW);
-	if(maxF)xDelete<IssmDouble>(maxF);
-	if(dW)xDelete<IssmDouble>(dW);
-	if(exsW)xDelete<IssmDouble>(exsW);
-	if(exsT)xDelete<IssmDouble>(exsT);
-	if(surpT)xDelete<IssmDouble>(surpT);
-	if(surpE)xDelete<IssmDouble>(surpE);
-	if(flxDn)xDelete<IssmDouble>(flxDn);
-	if(D)xDelete<int>(D);
-	if(M)xDelete<IssmDouble>(M);
  	xDelete<IssmDouble>(Zcum);
 	xDelete<IssmDouble>(dzMin2);
@@ -2058,11 +2131,4 @@
 
 	/*Assign output pointers:*/
-	*pMs=Msurf;
-	*pM=sumM;
-	*pR=Rsum;
-	*pF=Fsum;
-	*pmAdd=mAdd;
-	*pdz_add=dz_add;
-
 	*pT=T;
 	*pd=d;
@@ -2075,4 +2141,11 @@
 	*pgsp=gsp;
 	*pn=n;
+	*pm=m;
+	*pEI=EI;
+	*pEW=EW;
+
+	*pmAdd=mAdd;
+	*paddE=addE;
+	*pdz_add=dz_add;
 
 } /*}}}*/ 
Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h	(revision 27226)
+++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h	(revision 27227)
@@ -37,4 +37,5 @@
 void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** padiff, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, IssmDouble* pRa, int* pm, int aIdx, int dsnowIdx, IssmDouble Tmean, IssmDouble Ta, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, IssmDouble C, IssmDouble V, IssmDouble Vmean, IssmDouble dIce, int sid);
 void melt(IssmDouble* pM, IssmDouble* pMs, IssmDouble* pR, IssmDouble* pF, IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** padiff, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble Ra, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, IssmDouble zY, IssmDouble dIce, int sid);
+void managelayers(IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble* paddE, IssmDouble** pm, IssmDouble** pEI, IssmDouble** pEW, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** padiff, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, IssmDouble zY);
 void densification(IssmDouble** pd,IssmDouble** pdz, IssmDouble* T, IssmDouble* re, int denIdx, int aIdx, int swIdx, IssmDouble adThresh, IssmDouble C, IssmDouble dt, IssmDouble Tmean, IssmDouble dIce, int m, int sid);
 void turbulentFlux(IssmDouble* pshf, IssmDouble* plhf, IssmDouble* pEC, IssmDouble Ta, IssmDouble Ts, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble ds, IssmDouble Ws, IssmDouble Vz, IssmDouble Tz, IssmDouble dIce, int sid);