Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 27240)
+++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 27241)
@@ -101,5 +101,5 @@
 	}
 	#endif
-	if(dzTop < 0.05){
+	if(dzTop < 0.05-Dtol){
 		_printf_("initial top grid cell length (dzTop) is < 0.05 m\n");
 	}
@@ -108,5 +108,5 @@
 	//initialize top grid depth vector
 	dzT = xNew<IssmDouble>(gpTop); 
-	for (i=0;i<gpTop;i++)dzT[i]=dzTop;
+	for (i=0;i<gpTop-Dtol;i++)dzT[i]=dzTop;
 
 	//bottom grid cell depth = x*zY^(cells from to structure)
@@ -115,5 +115,5 @@
 	z0 = zTop;
 	gpBottom = 0;
-	while (zMax > z0){
+	while (zMax > z0+Dtol){
 		gp0= gp0 * zY;
 		z0 = z0 + gp0;
@@ -131,5 +131,5 @@
 	//combine top and bottom dz vectors
 	dz = xNew<IssmDouble>(gpTop+gpBottom);
-	for(i=0;i<gpTop;i++){
+	for(i=0;i<gpTop-Dtol;i++){
 		dz[i]=dzT[i];
 	}
@@ -215,6 +215,6 @@
 
 	//effective solar zenith angle
-	IssmDouble x = min(pow(t/(3*cos(Pi*SZA/180.0)),0.5), 1.0);
-	IssmDouble u = 0.64*x + (1-x)*cos(Pi*SZA/180.0);
+	IssmDouble x = min(pow(t/(3.0*cos(Pi*SZA/180.0)),0.5), 1.0);
+	IssmDouble u = 0.64*x + (1.0-x)*cos(Pi*SZA/180.0);
 
 	// pure snow albedo
@@ -243,5 +243,5 @@
 
 		// determine the effective change due to finite depth and soot loading
-		IssmDouble A = min(1.0, (2.1 * pow(z1,1.35*(1-as) - 0.1*c1 - 0.13)));
+		IssmDouble A = min(1.0, (2.1 * pow(z1,1.35*(1.0-as) - 0.1*c1 - 0.13)));
 
 		dac =  (as2 + dac2 - as) + A*((as + dac) - (as2 + dac2));
@@ -249,8 +249,8 @@
 
 	// change in albedo due to solar zenith angle
-	IssmDouble dasz = 0.53*as*(1 - (as + dac))*pow(1 - u,1.2);
+	IssmDouble dasz = 0.53*as*(1.0 - (as + dac))*pow(1.0 - u,1.2);
 
 	// change in albedo due to cloud (apart from change in diffuse fraction)
-	IssmDouble dat = (0.1*t*pow(as + dac,1.3)) / (pow(1 + 1.5*t,as));
+	IssmDouble dat = (0.1*t*pow(as + dac,1.3)) / (pow(1.0 + 1.5*t,as));
 
 	// Broadband albedo
@@ -551,6 +551,6 @@
 
 	//some constants:
-	const IssmDouble dSnow = 300;   // density of fresh snow [kg m-3]       
-	const IssmDouble dPHC = 830;  //Pore closeoff density
+	const IssmDouble dSnow = 300.0;   // density of fresh snow [kg m-3]       
+	const IssmDouble dPHC = 830.0;  //Pore closeoff density
 	const IssmDouble ai_max = 0.58;  //maximum ice albedo, from Lefebre,2003
 	const IssmDouble ai_min = aIce;  //minimum ice albedo
@@ -565,5 +565,5 @@
 			// clabSnow, IssmDouble clabIce, IssmDouble SZA, IssmDouble COT, int m
 			a[0]=gardnerAlb(re, dz, d, clabSnow, clabIce, SZA, COT, dPHC, m);
-			adiff[0]=gardnerAlb(re, dz, d, clabSnow, clabIce, 50, COT, dPHC, m);
+			adiff[0]=gardnerAlb(re, dz, d, clabSnow, clabIce, 50.0, COT, dPHC, m);
 		}
 		else if(aIdx==2){
@@ -688,5 +688,5 @@
 					//  before run-off in kg m^-2 (melt per GEMB timestep, i.e. 3 hourly)
 					IssmDouble M = Msurf+W[0];
-					a[0]=max(ai_min + (ai_max - ai_min)*exp(-1*(M/200)), ai_min);
+					a[0]=max(ai_min + (ai_max - ai_min)*exp(-1.0*(M/200.0)), ai_min);
 
 				}
@@ -696,6 +696,6 @@
 
 	// Check for erroneous values
-	if (a[0] > 1) _printf_("albedo > 1.0\n");
-	else if (a[0] < 0) _printf_("albedo is negative\n");
+	if (a[0] > 1.0+Ttol) _printf_("albedo > 1.0\n");
+	else if (a[0] < 0.0-Dtol) _printf_("albedo is negative\n");
 	else if (xIsNan(a[0])) _error_("albedo == NAN\n");
 
@@ -959,5 +959,5 @@
 
 	/* CALCULATE ENERGY SOURCES AND DIFFUSION FOR EVERY TIME STEP [dt]*/
-	for (IssmDouble i=1;i<=dt0;i+=dt){
+	for (IssmDouble i=1.0;i<=dt0;i+=dt){
 
 		// PART OF ENERGY CONSERVATION CHECK
@@ -1068,6 +1068,6 @@
             }
             //zL = min(zL, 0.5); //Sjoblom, 2014
-            zLM=max(zL/Vz*z0,1.e-3);
-            zLT=max(zL/Tz*zT,1.e-3);
+            zLM=max(zL/Vz*z0,1e-3);
+            zLT=max(zL/Tz*zT,1e-3);
 
             //Ding et al. 2020, from Beljaars and Holtslag (1991)
@@ -1466,5 +1466,5 @@
 	// MAIN FUNCTION
 	// specify constants
-	IssmDouble dSnow = 150;    // density of snow [kg m-3]
+	IssmDouble dSnow = 150.0;    // density of snow [kg m-3]
 	IssmDouble reNew = 0.05;    // new snow grain size [mm]
 	IssmDouble gdnNew = 1.0;     // new snow dendricity 
@@ -1474,7 +1474,7 @@
 	IssmDouble* mInit=NULL;
 	bool        top=true;
-	IssmDouble  mass=0;
-	IssmDouble  massinit=0;
-	IssmDouble  mass_diff=0;
+	IssmDouble  mass=0.0;
+	IssmDouble  massinit=0.0;
+	IssmDouble  mass_diff=0.0;
 
 	/*output: */
@@ -1512,10 +1512,10 @@
 
 		case 1: // Density of Antarctica snow
-			dSnow = 350;
-			//dSnow = 360; //FirnMICE Lundin et al., 2017
+			dSnow = 350.0;
+			//dSnow = 360.0; //FirnMICE Lundin et al., 2017
 			break;
 
 		case 2: // Density of Greenland snow, Fausto et al., 2018
-			dSnow = 315;
+			dSnow = 315.0;
 			//From Vionnet et al., 2012 (Crocus)
 			gdnNew = min(max(1.29 - 0.17*V,0.20),1.0);
@@ -1527,5 +1527,5 @@
 			//dSnow = alpha1 + beta1*T + delta1*C + epsilon1*W
 			//     7.36x10-2  1.06x10-3  6.69x10-2  4.77x10-3 
-			dSnow=(7.36e-2 + 1.06e-3*min(Tmean,CtoK-Ttol) + 6.69e-2*C/1000 + 4.77e-3*Vmean)*1000;
+			dSnow=(7.36e-2 + 1.06e-3*min(Tmean,CtoK-Ttol) + 6.69e-2*C/1000. + 4.77e-3*Vmean)*1000.;
 			break;
 
@@ -1728,9 +1728,9 @@
 	// specify irreducible water content saturation [fraction]
 	const IssmDouble Swi = 0.07;                     // assumed constant after Colbeck, 1974
-	const IssmDouble dPHC = 830;                     //Pore closeoff density
+	const IssmDouble dPHC = 830.0;                     //Pore closeoff density
 
 	//// REFREEZE PORE WATER
 	// check if any pore water
-	if (cellsum(W,n) >0.0+Wtol){
+	if (cellsum(W,n) > 0.0+Wtol){
 		if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("      pore water refreeze\n");
 		// calculate maximum freeze amount, maxF [kg]
@@ -1896,5 +1896,5 @@
 					dMax = (dIce - d[i])*dz_0;          // maximum refreeze                                             
 					IssmDouble maxF2 = min(dMax, maxF[i]-F1);      // maximum refreeze
-					F2 = min(-1*dW[i], maxF2);            // pore water refreeze
+					F2 = min(-1.0*dW[i], maxF2);            // pore water refreeze
 					m[i] = m[i] + F2;                   // mass after refreeze
 					d[i] = m[i]/dz_0;
@@ -2156,14 +2156,14 @@
 	for (int i=0;i<n;i++){
 		if (i==0){
-			dzMax2[i]=dzMin*2;
+			dzMax2[i]=dzMin*2.0;
 		}
 		else{
 			Zcum[i]=Zcum[i-1]+dz[i];
 			if (Zcum[i]<=zTop+Dtol){
-				dzMax2[i]=dzMin*2;
+				dzMax2[i]=dzMin*2.0;
 				X=i;
 			}
 			else{
-				dzMax2[i]=max(zY2*dzMin2[i-1],dzMin*2);
+				dzMax2[i]=max(zY2*dzMin2[i-1],dzMin*2.0);
 			}
 		}
@@ -2416,5 +2416,5 @@
 				// ERA5 new aIdx=1, swIdx=0
 				if (aIdx==1 && swIdx==0){
-					if (fabs(adThresh - 820) < Dtol){
+					if (fabs(adThresh - 820.0) < Dtol){
                   // ERA5 v4
                   M0 = max(1.5131 - (0.1317 * log(C)),0.25);
@@ -2461,5 +2461,5 @@
 				// ERA5 new aIdx=1, swIdx=0
 				if (aIdx==1 && swIdx==0){
-					if (fabs(adThresh - 820) < Dtol){
+					if (fabs(adThresh - 820.0) < Dtol){
 						// ERA5 v4
 						M0 = max(1.3566 - (0.1350 * log(C)),0.25);