Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 22613)
+++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 22614)
@@ -1247,5 +1247,4 @@
 	IssmDouble* surpT=NULL;
 	IssmDouble* surpE=NULL;
-	IssmDouble* F=NULL;
 	IssmDouble* flxDn=NULL;
 	IssmDouble  ER=0.0;
@@ -1311,6 +1310,6 @@
 	/*Allocations: */
 	M=xNewZeroInit<IssmDouble>(n); 
-	maxF=xNew<IssmDouble>(n); 
-	dW=xNew<IssmDouble>(n); 
+	maxF=xNewZeroInit<IssmDouble>(n); 
+	dW=xNewZeroInit<IssmDouble>(n); 
 
 	// store initial mass [kg] and energy [J]
@@ -1355,5 +1354,5 @@
 
 		// if pore water froze in ice then adjust d and dz thickness
-		for(int i=0;i<n;i++)if(d[i]>dIce)d[i]=dIce;
+		for(int i=0;i<n;i++)if(d[i]>dIce-Dtol)d[i]=dIce;
 		for(int i=0;i<n;i++) dz[i]= m[i]/d[i];
 	}
@@ -1375,7 +1374,7 @@
 		// (maximum T of snow before entire grid cell melts is a constant
 		// LF/CI = 159.1342)
-		surpT=xNew<IssmDouble>(n); for(int i=0;i<n;i++)surpT[i] = fmax(0.0, exsT[i]- 159.1342);
+		surpT=xNew<IssmDouble>(n); for(int i=0;i<n;i++)surpT[i] = fmax(0.0, exsT[i]- LF/CI);
         
-		if (cellsum(surpT,n) > 0.0 + Ttol ){
+		if (cellsum(surpT,n) > 0.0 + Ttol){
 			// _printf_("T Surplus");
 			// calculate surplus energy
@@ -1390,9 +1389,9 @@
 				T[i+1] = fmin(CtoK, T[i+1]);
                 
-				surpT[i+1] = fmax(0.0, exsT[i+1] - 159.1342);
+				surpT[i+1] = fmax(0.0, exsT[i+1] - LF/CI);
 				surpE[i+1] = surpT[i+1] * CI * m[i+1];
                 
 				// adjust current cell properties (again 159.1342 is the max T)
-				exsT[i] = 159.1342;
+				exsT[i] = LF/CI;
 				surpE[i] = 0.0;
 				i = i + 1;
@@ -1401,6 +1400,6 @@
 
 		// convert temperature excess to melt [kg]
-		for(int i=0;i<n;i++) M[i] = exsT[i] * d[i] * dz[i] * CI / LF;      // melt
-		sumM = cellsum(M,n);                                               // total melt [kg]
+		for(int i=0;i<n;i++) M[i] = fmin(exsT[i] * d[i] * dz[i] * CI / LF, m[i]);  // melt
+		sumM = cellsum(M,n);                                                       // total melt [kg]
 
 		// calculate maximum refreeze amount, maxF [kg]
@@ -1409,8 +1408,8 @@
 		// initialize refreeze, runoff, flxDn and dW vectors [kg]
  		IssmDouble* F = xNewZeroInit<IssmDouble>(n);
-		IssmDouble* R=xNewZeroInit<IssmDouble>(n);
+		IssmDouble* R = xNewZeroInit<IssmDouble>(n);
 
 		for(int i=0;i<n;i++)dW[i] = 0.0;
-		flxDn=xNewZeroInit<IssmDouble>(n+1); for(int i=0;i<n;i++)flxDn[i+1]=F[i];
+		flxDn=xNewZeroInit<IssmDouble>(n+1);
 
 		// determine the deepest grid cell where melt/pore water is generated
@@ -1442,6 +1441,7 @@
 				m[i] = m[i] - M[i];                     // mass after melt
 				Wi = (dIce-d[i]) * Swi * (m[i]/d[i]);    // irreducible water 
-				dW[i] = fmin(inM, Wi - W[i]);            // change in pore water
+				dW[i] = fmax(fmin(inM, Wi - W[i]),-1*W[i]);            // change in pore water
 				R[i] = fmax(0.0, inM - dW[i]);             // runoff
+				F[i] = 0.0;
 			}
 			// check if no energy to refreeze meltwater
@@ -1453,6 +1453,7 @@
 				m[i] = m[i] - M[i];                     // mass after melt
 				Wi = (dIce-d[i]) * Swi * (m[i]/d[i]);    // irreducible water 
-				dW[i] = fmin(inM, Wi-W[i]);              // change in pore water
-				flxDn[i+1] = fmax(0.0, inM-dW[i]);         // meltwater out
+				dW[i] = fmax(fmin(inM, Wi - W[i]),-1*W[i]);              // change in pore water
+				flxDn[i+1] = fmax(0.0, inM - dW[i]);         // meltwater out
+				R[i] = 0.0;
 				F[i] = 0.0;                               // no freeze 
 			}
@@ -1462,4 +1463,5 @@
 				// _printf_("MELT REFREEZE");
 				//-----------------------melt water-----------------------------
+				m[i] = m[i] - M[i];
 				IssmDouble dz_0 = m[i]/d[i];          
 				IssmDouble dMax = (dIce - d[i])*dz_0;              // d max = dIce
@@ -1471,5 +1473,5 @@
 				Wi = (dIce-d[i])* Swi * dz_0;            // irreducible water 
 				dW[i] = fmin(inM - F1, Wi-W[i]);         // change in pore water
-				if (-1*dW[i]>W[i]-Wtol ){
+				if (dW[i] < 0.0-Wtol && -1*dW[i]>W[i]-Wtol ){
 					dW[i]= -1*W[i];
 				}
@@ -1482,7 +1484,10 @@
 					m[i] = m[i] + F2;                   // mass after refreeze
 					d[i] = m[i]/dz_0;
+					dW[i] = dW[i] - F2;
 				}
 
-				flxDn[i+1] = inM - F1 - dW[i] - F2;     // meltwater out        
+				F[i] = F1 + F2;
+
+				flxDn[i+1] = inM - F1 - dW[i];     // meltwater out        
 				T[i] = T[i] + ((F1+F2)*(LF+(CtoK - T[i])*CI)/(m[i]*CI));// change in temperature
 
@@ -1507,9 +1512,12 @@
 		for(int i=0;i<n;i++)W[i] += dW[i];
 
+		//calculate Rsum:
+		Rsum=cellsum(R,n);
+
 		// delete all cells with zero mass
-		D_size=0; for(int i=0;i<n;i++)if(m[i]!=0)D_size++; 
+		D_size=0; for(int i=0;i<n;i++)if(m[i]>0.0)D_size++; 
 		D=xNew<int>(D_size);
-		D_size=0; for(int i=0;i<n;i++)if(m[i]!=0){ D[D_size] = i; D_size++;}
-		
+		D_size=0; for(int i=0;i<n;i++)if(m[i]>0.0){ D[D_size] = i; D_size++;}
+
 		celldelete(&m,n,D,D_size);
 		celldelete(&W,n,D,D_size);
@@ -1522,5 +1530,4 @@
 		celldelete(&EI,n,D,D_size);
 		celldelete(&EW,n,D,D_size);
-		celldelete(&R,n,D,D_size);
 		n=D_size;
 		xDelete<int>(D);
@@ -1528,7 +1535,4 @@
 		// calculate new grid lengths
 		for(int i=0;i<n;i++)dz[i] = m[i] / d[i];
-
-		//calculate Rsum:
-		Rsum=cellsum(R,n);
 
 		/*Free ressources:*/
@@ -1592,5 +1596,5 @@
             
 			// set cell to 99999 for deletion
-			m[i] = 99999;
+			m[i] = -99999;
 		}
 	}
@@ -1600,5 +1604,5 @@
          //find closest cell to merge with
 		for(int i=n-2;i>=0;i--){
-			if(m[i]!=99999){
+			if(m[i]!=-99999){
 				X2=i;
 				X1=n-1;
@@ -1622,11 +1626,11 @@
         
 		// set cell to 99999 for deletion
-		m[X2] = 99999;
+		m[X2] = -99999;
 	}
 
 	// delete combined cells
-	D_size=0; for(int i=0;i<n;i++)if(m[i]!=99999)D_size++; 
+	D_size=0; for(int i=0;i<n;i++)if(m[i]!=-99999)D_size++; 
 	D=xNew<int>(D_size); 
-	D_size=0; for(int i=0;i<n;i++)if(m[i]!=99999){ D[D_size] = i; D_size++;}
+	D_size=0; for(int i=0;i<n;i++)if(m[i]!=-99999){ D[D_size] = i; D_size++;}
 
 	celldelete(&m,n,D,D_size);
@@ -1723,5 +1727,5 @@
 		dz_add=-(dz[n-1]);
         
-		// add a grid cell of the same size and temperature to the bottom
+		// remove a grid cell from the bottom
 		D_size=n-1;
 		D=xNew<int>(D_size);
@@ -1756,5 +1760,5 @@
 	/*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_ADOLC_