Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 22248)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 22249)
@@ -2521,28 +2521,47 @@
 
 	/*Intermediary variables: {{{*/
-	IssmDouble isinitialized;
-	IssmDouble zTop,dzTop,zMax,zMin,zY,dzMin;
-	IssmDouble Tmean; 
-	IssmDouble C; 
-	IssmDouble Tz,Vz; 
+	IssmDouble isinitialized=0.0;
+	IssmDouble zTop=0.0;
+	IssmDouble dzTop=0.0;
+	IssmDouble zMax=0.0;
+	IssmDouble zMin=0.0;
+	IssmDouble zY=0.0;
+	IssmDouble dzMin=0.0;
+	IssmDouble Tmean=0.0;
+	IssmDouble C=0.0;
+	IssmDouble Tz,Vz=0.0;
 	IssmDouble rho_ice, rho_water,aSnow,aIce;
 	IssmDouble time,dt;
 	IssmDouble t,smb_dt;
 	IssmDouble yts;
-	IssmDouble Ta,V,dlw,dsw,P,eAir,pAir;
+	IssmDouble Ta=0.0;
+	IssmDouble V=0.0;
+	IssmDouble dlw=0.0;
+	IssmDouble dsw=0.0;
+	IssmDouble P=0.0;
+	IssmDouble eAir=0.0;
+	IssmDouble pAir=0.0;
 	int        aIdx=0;
 	int        denIdx=0;
 	int        swIdx=0;
 	IssmDouble cldFrac,t0wet, t0dry, K;
-	IssmDouble ulw;
-	IssmDouble netSW;
-	IssmDouble netLW;
-	IssmDouble lhf, shf, dayEC;
-	IssmDouble initMass;
-	IssmDouble sumR, sumM, sumEC, sumP, sumW,sumMassAdd;
-	IssmDouble sumdz_add;
-	IssmDouble sumMass,dMass;
+	IssmDouble ulw=0.0;
+	IssmDouble netSW=0.0;
+	IssmDouble netLW=0.0;
+	IssmDouble lhf=0.0;
+	IssmDouble shf=0.0;
+	IssmDouble dayEC=0.0;
+	IssmDouble initMass=0.0;
+	IssmDouble sumR=0.0;
+	IssmDouble sumM=0.0;
+	IssmDouble sumEC=0.0;
+	IssmDouble sumP=0.0;
+	IssmDouble sumW=0.0;
+	IssmDouble sumMassAdd=0.0;
+	IssmDouble sumdz_add=0.0;
+	IssmDouble sumMass=0.0;
+	IssmDouble dMass=0.0;
 	bool isgraingrowth,isalbedo,isshortwave,isthermal,isaccumulation,ismelt,isdensification,isturbulentflux;
-	IssmDouble init_scaling;
+	IssmDouble init_scaling=0.0;
 
 	/*}}}*/
@@ -2553,15 +2572,15 @@
 	IssmDouble* gdn = NULL;
 	IssmDouble* gsp = NULL;
-	IssmDouble  EC = 0;
+	IssmDouble  EC = 0.0;
 	IssmDouble* W = NULL;
 	IssmDouble* a = NULL;
 	IssmDouble* swf=NULL;
 	IssmDouble* T = NULL;
-	IssmDouble  T_bottom;
-	IssmDouble  M;
-	IssmDouble  R; 
-	IssmDouble  mAdd;
-	IssmDouble  dz_add;
-    
+	IssmDouble  T_bottom = 0.0;
+	IssmDouble  M = 0.0;
+	IssmDouble  R = 0.0;
+	IssmDouble  mAdd = 0.0;
+	IssmDouble  dz_add = 0.0;
+
 	IssmDouble* dzini=NULL;
 	IssmDouble* dini = NULL;
@@ -2572,6 +2591,6 @@
 	IssmDouble* aini = NULL;
 	IssmDouble* Tini = NULL;
-    
-	int         m;
+
+	int         m=0;
 	int         count=0;
 	/*}}}*/
@@ -2606,5 +2625,5 @@
 	parameters->FindParam(&isturbulentflux,SmbIsturbulentfluxEnum);
 	parameters->FindParam(&init_scaling,SmbInitDensityScalingEnum);
-    
+
 	/*}}}*/
 	/*Retrieve inputs: {{{*/
@@ -2645,102 +2664,102 @@
 	/*First, check that the initial structures have been setup in GEMB. If not, initialize profile variables: layer thickness dz, * density d, temperature T, etc. {{{*/
 	if(isinitialized==0.0){
-        if(VerboseSmb() && this->Sid()==0)_printf0_("smb core: Initializing grid\n");
-        //if(this->Sid()==1) for(int i=0;i<m;i++)_printf_("z[" << i << "]=" <<
-        //dz[i] << "\n");
-        
-        DoubleArrayInput* dz_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDziniEnum)); _assert_(dz_input);
-        DoubleArrayInput* d_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDiniEnum));_assert_(d_input);
-        DoubleArrayInput* re_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbReiniEnum));_assert_(re_input);
-        DoubleArrayInput* gdn_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGdniniEnum));_assert_(gdn_input);
-        DoubleArrayInput* gsp_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGspiniEnum));_assert_(gsp_input);
-        DoubleInput* EC_input= dynamic_cast<DoubleInput*>(this->GetInput(SmbECiniEnum));_assert_(EC_input);
-        DoubleArrayInput* W_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbWiniEnum));_assert_(W_input);
-        DoubleArrayInput* a_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbAiniEnum));_assert_(a_input);
-        DoubleArrayInput* T_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbTiniEnum));_assert_(T_input);
-
-        dz_input->GetValues(&dzini,&m);
-        d_input->GetValues(&dini,&m);
-        re_input->GetValues(&reini,&m);
-        gdn_input->GetValues(&gdnini,&m);
-        gsp_input->GetValues(&gspini,&m);
-        EC_input->GetInputValue(&EC);
-        W_input->GetValues(&Wini,&m);
-        a_input->GetValues(&aini,&m);
-        T_input->GetValues(&Tini,&m);
-        
-        /*Retrive the correct value of m (without the zeroes at the end)*/
-        Input* Size_input=this->GetInput(SmbSizeiniEnum); _assert_(Size_input);
-        Size_input->GetInputValue(&m);
-        
-        if(m==2){ //Snow properties are initialized with default values. Vertical grid has to be initialized too
-//            if(VerboseSmb() && this->Sid()==0)_printf0_("Snow properties initialized w DEFAULT values\n");
-            
-            /*initialize profile variables:*/
-            GembgridInitialize(&dz, &m, zTop, dzTop, zMax, zY);
-            
-            d = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)d[i]=dini[0]; //ice density [kg m-3]
-            re = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)re[i]=reini[0];         //set grain size to old snow [mm]
-            gdn = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)gdn[i]=gdnini[0];         //set grain dentricity to old snow
-            gsp = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)gsp[i]=gspini[0];         //set grain sphericity to old snow
-            W = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)W[i]=Wini[0];             //set water content to zero [kg m-2]
-            a = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)a[i]=aini[0];         //set albedo equal to fresh snow [fraction]
-            T = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)T[i]=Tmean;         //set initial grid cell temperature to the annual mean temperature [K]
-/*            /!\ Default value of T can not be retrived from SMBgemb.m (like other snow properties) because don't know Tmean yet when set default values.
-            Default value of 0C given in SMBgemb.m is overwritten here with value of Tmean*/
-            
-            //fixed lower temperature bounday condition - T is fixed
-            T_bottom=T[m-1];
-        }
-        else{ //Retrieve snow properties from previous run. Need to provide values for all layers
-//            if(VerboseSmb() && this->Sid()==0)_printf0_("Snow properties initialized w RESTART values\n");
-            
-            dz = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)dz[i]=dzini[i];
-            d = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)d[i]=dini[i];
-            re = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)re[i]=reini[i];
-            gdn = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)gdn[i]=gdnini[i];
-            gsp = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)gsp[i]=gspini[i];
-            W = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)W[i]=Wini[i];
-            a = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)a[i]=aini[i];
-            T = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)T[i]=Tini[i];
-
-            //fixed lower temperature bounday condition - T is fixed
-            T_bottom=T[m-1];
-        }
-        
-        /*Flag the initialization:*/
-        this->AddInput(new DoubleInput(SmbIsInitializedEnum,1.0));
-    }
-    else{
-        /*Recover inputs: */
-        DoubleArrayInput* dz_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDzEnum)); _assert_(dz_input);
-        DoubleArrayInput* d_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDEnum));_assert_(d_input);
-        DoubleArrayInput* re_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbReEnum));_assert_(re_input);
-        DoubleArrayInput* gdn_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGdnEnum));_assert_(gdn_input);
-        DoubleArrayInput* gsp_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGspEnum));_assert_(gsp_input);
-        DoubleInput* EC_input= dynamic_cast<DoubleInput*>(this->GetInput(SmbECEnum));_assert_(EC_input);
-        DoubleArrayInput* W_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbWEnum));_assert_(W_input);
-        DoubleArrayInput* a_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbAEnum));_assert_(a_input);
-        DoubleArrayInput* T_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbTEnum));_assert_(T_input);
-
-        /*Recover arrays: */
-        dz_input->GetValues(&dz,&m);
-        d_input->GetValues(&d,&m);
-        re_input->GetValues(&re,&m);
-        gdn_input->GetValues(&gdn,&m);
-        gsp_input->GetValues(&gsp,&m);
-        EC_input->GetInputValue(&EC);
-        W_input->GetValues(&W,&m);
-        a_input->GetValues(&a,&m);
-        T_input->GetValues(&T,&m);
-
-        //fixed lower temperature bounday condition - T is fixed
-        T_bottom=T[m-1];
-
-    } /*}}}*/
+		if(VerboseSmb() && this->Sid()==0)_printf0_("smb core: Initializing grid\n");
+		//if(this->Sid()==1) for(int i=0;i<m;i++)_printf_("z[" << i << "]=" <<
+		//dz[i] << "\n");
+
+		DoubleArrayInput* dz_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDziniEnum)); _assert_(dz_input);
+		DoubleArrayInput* d_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDiniEnum));_assert_(d_input);
+		DoubleArrayInput* re_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbReiniEnum));_assert_(re_input);
+		DoubleArrayInput* gdn_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGdniniEnum));_assert_(gdn_input);
+		DoubleArrayInput* gsp_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGspiniEnum));_assert_(gsp_input);
+		DoubleInput* EC_input= dynamic_cast<DoubleInput*>(this->GetInput(SmbECiniEnum));_assert_(EC_input);
+		DoubleArrayInput* W_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbWiniEnum));_assert_(W_input);
+		DoubleArrayInput* a_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbAiniEnum));_assert_(a_input);
+		DoubleArrayInput* T_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbTiniEnum));_assert_(T_input);
+
+		dz_input->GetValues(&dzini,&m);
+		d_input->GetValues(&dini,&m);
+		re_input->GetValues(&reini,&m);
+		gdn_input->GetValues(&gdnini,&m);
+		gsp_input->GetValues(&gspini,&m);
+		EC_input->GetInputValue(&EC);
+		W_input->GetValues(&Wini,&m);
+		a_input->GetValues(&aini,&m);
+		T_input->GetValues(&Tini,&m);
+
+		/*Retrive the correct value of m (without the zeroes at the end)*/
+		Input* Size_input=this->GetInput(SmbSizeiniEnum); _assert_(Size_input);
+		Size_input->GetInputValue(&m);
+
+		if(m==2){ //Snow properties are initialized with default values. Vertical grid has to be initialized too
+			//            if(VerboseSmb() && this->Sid()==0)_printf0_("Snow properties initialized w DEFAULT values\n");
+
+			/*initialize profile variables:*/
+			GembgridInitialize(&dz, &m, zTop, dzTop, zMax, zY);
+
+			d = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)d[i]=dini[0]; //ice density [kg m-3]
+			re = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)re[i]=reini[0];         //set grain size to old snow [mm]
+			gdn = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)gdn[i]=gdnini[0];         //set grain dentricity to old snow
+			gsp = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)gsp[i]=gspini[0];         //set grain sphericity to old snow
+			W = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)W[i]=Wini[0];             //set water content to zero [kg m-2]
+			a = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)a[i]=aini[0];         //set albedo equal to fresh snow [fraction]
+			T = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)T[i]=Tmean;         //set initial grid cell temperature to the annual mean temperature [K]
+			/*            /!\ Default value of T can not be retrived from SMBgemb.m (like other snow properties) because don't know Tmean yet when set default values.
+							  Default value of 0C given in SMBgemb.m is overwritten here with value of Tmean*/
+
+			//fixed lower temperature bounday condition - T is fixed
+			T_bottom=T[m-1];
+		}
+		else{ //Retrieve snow properties from previous run. Need to provide values for all layers
+			//            if(VerboseSmb() && this->Sid()==0)_printf0_("Snow properties initialized w RESTART values\n");
+
+			dz = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)dz[i]=dzini[i];
+			d = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)d[i]=dini[i];
+			re = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)re[i]=reini[i];
+			gdn = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)gdn[i]=gdnini[i];
+			gsp = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)gsp[i]=gspini[i];
+			W = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)W[i]=Wini[i];
+			a = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)a[i]=aini[i];
+			T = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)T[i]=Tini[i];
+
+			//fixed lower temperature bounday condition - T is fixed
+			T_bottom=T[m-1];
+		}
+
+		/*Flag the initialization:*/
+		this->AddInput(new DoubleInput(SmbIsInitializedEnum,1.0));
+	}
+	else{
+		/*Recover inputs: */
+		DoubleArrayInput* dz_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDzEnum)); _assert_(dz_input);
+		DoubleArrayInput* d_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDEnum));_assert_(d_input);
+		DoubleArrayInput* re_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbReEnum));_assert_(re_input);
+		DoubleArrayInput* gdn_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGdnEnum));_assert_(gdn_input);
+		DoubleArrayInput* gsp_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGspEnum));_assert_(gsp_input);
+		DoubleInput* EC_input= dynamic_cast<DoubleInput*>(this->GetInput(SmbECEnum));_assert_(EC_input);
+		DoubleArrayInput* W_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbWEnum));_assert_(W_input);
+		DoubleArrayInput* a_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbAEnum));_assert_(a_input);
+		DoubleArrayInput* T_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbTEnum));_assert_(T_input);
+
+		/*Recover arrays: */
+		dz_input->GetValues(&dz,&m);
+		d_input->GetValues(&d,&m);
+		re_input->GetValues(&re,&m);
+		gdn_input->GetValues(&gdn,&m);
+		gsp_input->GetValues(&gsp,&m);
+		EC_input->GetInputValue(&EC);
+		W_input->GetValues(&W,&m);
+		a_input->GetValues(&a,&m);
+		T_input->GetValues(&T,&m);
+
+		//fixed lower temperature bounday condition - T is fixed
+		T_bottom=T[m-1];
+
+	} /*}}}*/
 
 	// determine initial mass [kg]
 	initMass=0; for(int i=0;i<m;i++) initMass += dz[i]*d[i] + W[i];
-    
-        // initialize cumulative variables
+
+	// initialize cumulative variables
 	sumR = 0; sumM = 0; sumEC = 0; sumP = 0; sumMassAdd = 0;
 	sumdz_add=0;
@@ -2771,30 +2790,30 @@
 
 		/*Snow, firn and ice albedo:*/
-		if(isalbedo)albedo(a,aIdx,re,d,cldFrac,aIce, aSnow,T,W,P,EC,t0wet,t0dry,K,smb_dt,m,this->Sid());
-		
-					
+		if(isalbedo)albedo(a,aIdx,re,d,cldFrac,aIce, aSnow,T,W,P,EC,t0wet,t0dry,K,smb_dt,rho_ice,m,this->Sid());
+
+
 		/*Distribution of absorbed short wave radation with depth:*/
-		if(isshortwave)shortwave(&swf, swIdx, aIdx, dsw, a[0], d, dz, re,m,this->Sid());
-		
+		if(isshortwave)shortwave(&swf, swIdx, aIdx, dsw, a[0], d, dz, re,rho_ice,m,this->Sid());
+
 		/*Calculate net shortwave [W m-2]*/
 		netSW = cellsum(swf,m);
 
 		/*Thermal profile computation:*/
-		if(isthermal)thermo(&EC, T, dz, d, swf, dlw, Ta, V, eAir, pAir, W[0], smb_dt, m, Vz, Tz,this->Sid());
+		if(isthermal)thermo(&EC, T, dz, d, swf, dlw, Ta, V, eAir, pAir, W[0], smb_dt, m, Vz, Tz,rho_ice,this->Sid());
 
 		/*Change in thickness of top cell due to evaporation/condensation  assuming same density as top cell. 
 		 * need to fix this in case all or more of cell evaporates */
 		dz[0] = dz[0] + EC / d[0];
-		
+
 		/*Add snow/rain to top grid cell adjusting cell depth, temperature and density*/
-		if(isaccumulation)accumulation(&T, &dz, &d, &W, &a, &re, &gdn, &gsp, &m, Ta, P, dzMin, aSnow,this->Sid());
+		if(isaccumulation)accumulation(&T, &dz, &d, &W, &a, &re, &gdn, &gsp, &m, Ta, P, dzMin, aSnow,rho_ice,this->Sid());
 
 		/*Calculate water production, M [kg m-2] resulting from snow/ice temperature exceeding 273.15 deg K
 		 * (> 0 deg C), runoff R [kg m-2] and resulting changes in density and determine wet compaction [m]*/
-		if(ismelt)melt(&M, &R, &mAdd, &dz_add, &T, &d, &dz, &W, &a, &re, &gdn, &gsp, &m, dzMin, zMax, zMin, zTop,this->Sid());
+		if(ismelt)melt(&M, &R, &mAdd, &dz_add, &T, &d, &dz, &W, &a, &re, &gdn, &gsp, &m, dzMin, zMax, zMin, zTop,rho_ice,this->Sid());
 
 		/*Allow non-melt densification and determine compaction [m]*/
 		if(isdensification)densification(d,dz, T, re, denIdx, C, smb_dt, Tmean,rho_ice,m,this->Sid());
-		
+
 		/*Calculate upward longwave radiation flux [W m-2] not used in energy balance. Calculated for every 
 		 * sub-time step in thermo equations*/
@@ -2803,23 +2822,23 @@
 		/*Calculate net longwave [W m-2]*/
 		netLW = dlw - ulw;
-		
+
 		/*Calculate turbulent heat fluxes [W m-2]*/
-		if(isturbulentflux)turbulentFlux(&shf, &lhf, &dayEC, Ta, T[0], V, eAir, pAir, d[0], W[0], Vz, Tz,this->Sid());
-		
+		if(isturbulentflux)turbulentFlux(&shf, &lhf, &dayEC, Ta, T[0], V, eAir, pAir, d[0], W[0], Vz, Tz,rho_ice,this->Sid());
+
 		/*Verbose some resuls in debug mode: {{{*/
 		if(VerboseSmb() && 0){ 
 			_printf_("smb log: count[" << count << "] m[" << m << "] " 
-				<< setprecision(16)   << "T[" << cellsum(T,m)  << "] " 
-					                  << "d[" << cellsum(d,m)  << "] "
-					                  << "dz[" << cellsum(dz,m)  << "] "
-					                  << "a[" << cellsum(a,m)  << "] "
-					                  << "W[" << cellsum(W,m)  << "] "
-					                  << "re[" << cellsum(re,m)  << "] "
-					                  << "gdn[" << cellsum(gdn,m)  << "] "
-					                  << "gsp[" << cellsum(gsp,m)  << "] "
-					                  << "swf[" << netSW << "] "
-									  << "\n");
+						<< setprecision(16)   << "T[" << cellsum(T,m)  << "] " 
+						<< "d[" << cellsum(d,m)  << "] "
+						<< "dz[" << cellsum(dz,m)  << "] "
+						<< "a[" << cellsum(a,m)  << "] "
+						<< "W[" << cellsum(W,m)  << "] "
+						<< "re[" << cellsum(re,m)  << "] "
+						<< "gdn[" << cellsum(gdn,m)  << "] "
+						<< "gsp[" << cellsum(gsp,m)  << "] "
+						<< "swf[" << netSW << "] "
+						<< "\n");
 		} /*}}}*/
-		
+
 		/*Sum component mass changes [kg m-2]*/
 		sumMassAdd = mAdd + sumMassAdd;
@@ -2841,8 +2860,8 @@
 		if (dMass != 0.0) _printf_("total system mass not conserved in MB function \n");
 		#endif
-		
+
 		/*Check bottom grid cell T is unchanged:*/
 		if (T[m-1]!=T_bottom) _printf_("T(end)~=T_bottom" << "\n");
-		
+
 		/*Free ressources: */
 		xDelete<IssmDouble>(swf);
@@ -2859,12 +2878,13 @@
 	this->AddInput(new DoubleArrayInput(SmbGspEnum,gsp,m));
 	this->AddInput(new DoubleArrayInput(SmbTEnum,T,m));
-	this->AddInput(new DoubleInput(SmbECEnum,sumEC/yts));
+	this->AddInput(new DoubleInput(SmbECEnum,sumEC/dt/rho_ice));
 	this->AddInput(new DoubleArrayInput(SmbWEnum,W,m));
 	this->AddInput(new DoubleArrayInput(SmbAEnum,a,m));
-	this->AddInput(new DoubleInput(SmbMassBalanceEnum,(sumP + sumEC -sumR)/yts));
-	this->AddInput(new DoubleInput(SmbRunoffEnum,sumR/yts));
-	this->AddInput(new DoubleInput(SmbPrecipitationEnum,sumP/yts));
+	this->AddInput(new DoubleInput(SmbMassBalanceEnum,(sumP + sumEC -sumR)/dt/rho_ice));
+	this->AddInput(new DoubleInput(SmbMeltEnum,sumM/dt/rho_ice));
+	this->AddInput(new DoubleInput(SmbRunoffEnum,sumR/dt/rho_ice));
+	this->AddInput(new DoubleInput(SmbPrecipitationEnum,sumP/dt/rho_ice));
 	this->AddInput(new DoubleInput(SmbDz_addEnum,sumdz_add/yts));
-	this->AddInput(new DoubleInput(SmbM_addEnum,sumMassAdd/yts));
+	this->AddInput(new DoubleInput(SmbM_addEnum,sumMassAdd/dt));
 
 	/*Free allocations:{{{*/
Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 22248)
+++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 22249)
@@ -8,4 +8,22 @@
 
 const double Pi = 3.141592653589793;
+const double CtoK = 273.15;             // Kelvin to Celcius conversion/ice melt. point T in K
+const double dts = 86400.0;              // Number of seconds in a day
+
+/* Tolerances have to be defined for if loops involving some specific values
+	like densitiy of ice or melting point temp because values are "rounded"
+	(e.g rho ice = 909.99... instead of 910.0) and we don't always go into the right loop */
+const double Ttol = 1e-10;
+const double Dtol = 1e-11;
+const double Gdntol = 1e-10;
+const double Wtol = 1e-13;
+
+const double CI = 2102.0;                       // heat capacity of snow/ice (J kg-1 k-1)
+const double LF = 0.3345E6;             // latent heat of fusion (J kg-1)
+const double LV = 2.495E6;               // latent heat of vaporization (J kg-1)
+const double LS = 2.8295E6;             // latent heat of sublimation (J kg-1)
+const double SB = 5.67E-8;                // Stefan-Boltzmann constant (W m-2 K-4)
+const double CA = 1005.0;                    // heat capacity of air (J kg-1 k-1)
+const double R = 8.314;                      // gas constant (mol-1 K-1)
 
 void Gembx(FemModel* femmodel){  /*{{{*/
@@ -25,8 +43,10 @@
 
 	/*intermediary:*/
-	IssmDouble dgpTop;
-	int gpTop, gpBottom;
-	int i;
-	IssmDouble gp0,z0;
+	IssmDouble dgpTop=0.0;
+	int gpTop=0;
+	int gpBottom=0;
+	int i=0;
+	IssmDouble gp0=0.0;
+	IssmDouble z0=0.0;
 	IssmDouble* dzT=NULL;
 	IssmDouble* dzB=NULL;
@@ -101,12 +121,12 @@
 
 	// initialize
-	IssmDouble F = 0, H=0, G=0;
+	IssmDouble F = 0.0, H=0.0, G=0.0;
 	const IssmDouble E = 0.09;        //[mm d-1] model time growth constant E
 	// convert T from K to ºC
-	T = T - 273.15;
+	T = T - CtoK;
 
 	// temperature coefficient F
 	if(T>-6.0) F =  0.7 + ((T/-6.0) * 0.3);
-	if(T<=-6.0 && T>-22.0) F =  1 - ((T+6.0)/-16.0 * 0.8);
+	if(T<=-6.0 && T>-22.0) F =  1.0 - ((T+6.0)/-16.0 * 0.8);
 	if(T<=-22.0 && T>-40.0) F =  0.2 - ((T+22.0)/-18.0 * 0.2);
 
@@ -114,5 +134,5 @@
 	if(d<150.0) H=1.0;
 
-	if(d>=150.0 & d<400.0) H = 1 - ((d-150.0)/250.0);
+	if(d>=150.0 & d<400.0) H = 1.0 - ((d-150.0)/250.0);
 
 	// temperature gradient coefficient G
@@ -121,5 +141,5 @@
 	if(dT >= 0.4  && dT < 0.5)  G = 0.67 + (((dT - 0.4)/0.1) * 0.23);
 	if(dT >= 0.5  && dT < 0.7)  G = 0.9 + (((dT - 0.5)/0.2) * 0.1);
-	if(dT >= .7              )  G = 1.0;
+	if(dT >= 0.7              )  G = 1.0;
 
 	// grouped coefficient Q
@@ -164,10 +184,10 @@
 
 	/*intermediary: */
-	IssmDouble  dt;
+	IssmDouble  dt=0.0;
 	IssmDouble* gsz=NULL;
 	IssmDouble* dT=NULL;
 	IssmDouble* zGPC=NULL;
 	IssmDouble* lwc=NULL;
-	IssmDouble  Q=0;
+	IssmDouble  Q=0.0;
 
 	if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   grain growth module\n");
@@ -183,5 +203,5 @@
 
 	/*Convert dt from seconds to day: */
-	dt=smb_dt/86400.0;
+	dt=smb_dt/dts;
 
 	/*Determine liquid-water content in percentage: */
@@ -189,5 +209,5 @@
 
 	//set maximum water content by mass to 9 percent (Brun, 1980)
-	for(int i=0;i<m;i++)if(lwc[i]>9.0) lwc[i]=9.0;
+	for(int i=0;i<m;i++)if(lwc[i]>9.0-Wtol) lwc[i]=9.0;
 
 
@@ -202,5 +222,5 @@
 	for(int i=0;i<m;i++){
 		for (int j=0;j<=i;j++) zGPC[i]+=dz[j];
-		zGPC[i]-=dz[i]/2;
+		zGPC[i]-=(dz[i]/2.0);
 	}
 
@@ -217,10 +237,10 @@
 	/*Snow metamorphism. Depends on value of snow dendricity and wetness of the snowpack: */
 	for(int i=0;i<m;i++){
-		if (gdn[i]>0){
+		if (gdn[i]>0.0+Gdntol){
 
 			//_printf_("Dendritic dry snow metamorphism\n");
 
 			//index for dentricity > 0 and T gradients < 5 degC m-1 and >= 5 degC m-1
-			if(dT[i]<5.0){
+			if(dT[i]<5.0-Ttol){
 				//determine coefficients
 				IssmDouble A = - 2e8 * exp(-6e3 / T[i]) * dt;
@@ -240,5 +260,5 @@
 
 			// wet snow metamorphism
-			if(W[i]>0.0){
+			if(W[i]>0.0+Wtol){
 
 				//_printf_("D}ritic wet snow metamorphism\n");
@@ -253,11 +273,11 @@
 
 			// dendricity and sphericity can not be > 1 or < 0
-			if (gdn[i]<0.0)gdn[i]=0.0;
-			if (gsp[i]<0.0)gsp[i]=0.0;
-			if (gdn[i]>1.0)gdn[i]=1.0;
-			if (gsp[i]>1.0)gsp[i]=1.0;
+			if (gdn[i]<0.0+Wtol)gdn[i]=0.0;
+			if (gsp[i]<0.0+Wtol)gsp[i]=0.0;
+			if (gdn[i]>1.0-Wtol)gdn[i]=1.0;
+			if (gsp[i]>1.0-Wtol)gsp[i]=1.0;
 
 			// determine new grain size (mm)
-			gsz[i] = 0.1 + (1-gdn[i])*0.25 + (0.5-gsp[i])*0.1;
+			gsz[i] = 0.1 + (1.0-gdn[i])*0.25 + (0.5-gsp[i])*0.1;
 
 		}
@@ -271,11 +291,11 @@
 
 			// calculate grain growth
-			gsz[i] += Q* dt;
+			gsz[i] += (Q*dt);
 
 			//Wet snow metamorphism (Brun, 1989)
-			if(W[i]>0.0){
+			if(W[i]>0.0+Wtol){
 				//_printf_("Nond}ritic wet snow metamorphism\n");
 				//wet rate of change coefficient
-				IssmDouble E = 1.28e-8 + (4.22e-10 * pow(lwc[i],3.0))* (dt *86400.0);   // [mm^3 s^-1]
+				IssmDouble E = 1.28e-8 + (4.22e-10 * pow(lwc[i],3.0))* (dt *dts);   // [mm^3 s^-1]
 
 				// calculate change in grain volume and convert to grain size
@@ -285,8 +305,8 @@
 
 			// grains with sphericity == 1 can not have grain sizes > 2 mm (Brun, 1992)
-			if (gsp[i]==1.0 && gsz[i]>2.0) gsz[i]=2.0;
+			if (fabs(gsp[i]-1.0)<Wtol && gsz[i]>2.0-Wtol) gsz[i]=2.0;
 
 			// grains with sphericity == 0 can not have grain sizes > 5 mm (Brun, 1992)
-			if (gsp[i]!=1.0 && gsz[i]>5.0) gsz[i]=5.0;
+			if (fabs(gsp[i]-1.0)>Wtol && gsz[i]>5.0-Wtol) gsz[i]=5.0;
 		}
 
@@ -302,5 +322,5 @@
 
 }  /*}}}*/
-void albedo(IssmDouble* a, int aIdx, IssmDouble* re, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble* TK, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, int m,int sid) { /*{{{*/
+void albedo(IssmDouble* a, int aIdx, IssmDouble* re, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble* TK, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m,int sid) { /*{{{*/
 
 	//// Calculates Snow, firn and ice albedo as a function of:
@@ -350,5 +370,4 @@
 	//some constants:
 	const IssmDouble dSnow = 300;   // density of fresh snow [kg m-3]       
-	const IssmDouble dIce = 910;    // density of ice [kg m-3]
 
 	if(aIdx==1){ 
@@ -356,8 +375,8 @@
         
         //convert effective radius to specific surface area [cm2 g-1]
-        IssmDouble S = 3.0 / (.091 * re[0]);
+        IssmDouble S = 3.0 / (0.091 * re[0]);
         
         //determine broadband albedo
-        a[0]= 1.48 - pow(S,-.07);
+        a[0]= 1.48 - pow(S,-0.07);
 	}
 	else if(aIdx==2){
@@ -369,5 +388,5 @@
         
         // convert effective radius to grain size in meters
-        IssmDouble gsz = (re[0] * 2) / 1000.0;
+        IssmDouble gsz = (re[0] * 2.0) / 1000.0;
         
         // spectral range:
@@ -398,5 +417,5 @@
         // where: t0 = timescale for albedo decay
         
-        dt = dt / 86400;    // convert from [s] to [d]
+        dt = dt / dts;    // convert from [s] to [d]
         
         // initialize variables
@@ -414,12 +433,12 @@
         // K = 7                // time scale temperature coef. (7) [d]
         // W0 = 300;            // 200 - 600 [mm]
-        const IssmDouble z_snow = 15;            // 16 - 32 [mm]
+        const IssmDouble z_snow = 15.0;            // 16 - 32 [mm]
         
         // determine timescale for albedo decay
-        for(int i=0;i<m;i++)if(W[i]>0)t0[i]=t0wet; // wet snow timescale
-        for(int i=0;i<m;i++)T[i]=TK[i] - 273.15; // change T from K to degC
+        for(int i=0;i<m;i++)if(W[i]>0.0+Wtol)t0[i]=t0wet; // wet snow timescale
+        for(int i=0;i<m;i++)T[i]=TK[i] - CtoK; // change T from K to degC
         for(int i=0;i<m;i++) t0warm[i]= fabs(T[i]) * K + t0dry; //// 'warm' snow timescale
-        for(int i=0;i<m;i++)if(W[i]==0.0 && T[i]>=-10)t0[i]= t0warm[i];
-        for(int i=0;i<m;i++)if(T[i]<-10) t0[i] =  10 * K + t0dry; // 'cold' snow timescale
+        for(int i=0;i<m;i++)if(fabs(W[i])<Wtol && T[i]>=-10.0-Ttol)t0[i]= t0warm[i];
+        for(int i=0;i<m;i++)if(T[i]<-10.0-Ttol) t0[i] =  10.0 * K + t0dry; // 'cold' snow timescale
         
         // calculate new albedo
@@ -431,5 +450,5 @@
         
         // check if condensation occurs & if it is deposited in solid phase
-        if ( EC > 0 && T[0] < 0) P = P + (EC/dSnow) * 1000;  // add cond to precip [mm]
+        if ( EC > 0.0 + Dtol && T[0] < 0.0-Ttol) P = P + (EC/dSnow) * 1000.0;  // add cond to precip [mm]
         
         a[0] = aSnow - (aSnow - a[0]) * exp(-P/z_snow);
@@ -453,5 +472,5 @@
 	else if (xIsNan(a[0])) _error_("albedo == NAN\n");
 }  /*}}}*/
-void thermo(IssmDouble* pEC, IssmDouble* T, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlwrf, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz,int sid) { /*{{{*/
+void thermo(IssmDouble* pEC, IssmDouble* T, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlwrf, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble dIce, int sid) { /*{{{*/
 
 	/* ENGLACIAL THERMODYNAMICS*/
@@ -495,53 +514,41 @@
 	IssmDouble* sw = NULL;
 	IssmDouble* dT_sw = NULL;
-	IssmDouble* lw = NULL;
 	IssmDouble* T0 = NULL;
 	IssmDouble* Tu = NULL;
 	IssmDouble* Td = NULL;
 
-	IssmDouble  z0;	
-	IssmDouble  dt;
-	IssmDouble max_fdt=0;
-	IssmDouble  Ts=0;
-	IssmDouble  L;
-	IssmDouble  eS;
-	IssmDouble  Ri=0;
-	IssmDouble  coefM;
-	IssmDouble  coefH;
-	IssmDouble An;
-	IssmDouble C;
-	IssmDouble shf;
-	IssmDouble SB;
-	IssmDouble CI; 
-	IssmDouble ds;
-	IssmDouble dAir;
-	IssmDouble TCs;
-	IssmDouble lhf;
-	IssmDouble EC_day;
-	IssmDouble dT_turb;
-	IssmDouble turb;
-	IssmDouble ulw;
-	IssmDouble dT_ulw;
-	IssmDouble dlw;
-	IssmDouble dT_dlw;
+	IssmDouble  z0=0.0;	
+	IssmDouble  dt=0.0;
+	IssmDouble max_fdt=0.0;
+	IssmDouble  Ts=0.0;
+	IssmDouble  L=0.0;
+	IssmDouble  eS=0.0;
+	IssmDouble  Ri=0.0;
+	IssmDouble  coefM=0.0;
+	IssmDouble  coefH=0.0;
+	IssmDouble An=0.0;
+	IssmDouble C=0.0;
+	IssmDouble shf=0.0;
+	IssmDouble ds=0.0;
+	IssmDouble dAir=0.0;
+	IssmDouble TCs=0.0;
+	IssmDouble lhf=0.0;
+	IssmDouble EC_day=0.0;
+	IssmDouble dT_turb=0.0;
+	IssmDouble turb=0.0;
+	IssmDouble ulw=0.0;
+	IssmDouble dT_ulw=0.0;
+	IssmDouble dlw=0.0;
+	IssmDouble dT_dlw=0.0;
 	
 	/*outputs:*/
-	IssmDouble EC;
+	IssmDouble EC=0.0;
 
 	if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   thermal module\n");
 
-	// INITIALIZE
-	CI = 2102;          // heat capacity of snow/ice (J kg-1 k-1)
-	// CA = 1005;                  // heat capacity of air (J kg-1 k-1)
-	// LF = 0.3345E6;              // latent heat of fusion(J kg-1)
-	// LV = 2.495E6;               // latent heat of vaporization(J kg-1)
-	// dIce = 910;                 // density of ice [kg m-3]
-	// dSnow = 300;                // density of snow [kg m-3]
-	SB = 5.67E-8;       // Stefan-Boltzmann constant [W m-2 K-4]
-
 	ds = d[0];      // density of top grid cell
 
 	// calculated air density [kg/m3]
-	dAir = 0.029 * pAir /(8.314 * Ta);
+	dAir = 0.029 * pAir /(R * Ta);
 
 	// thermal capacity of top grid cell [J/k]
@@ -549,5 +556,5 @@
 
 	//initialize Evaporation - Condenstation 
-	EC = 0;
+	EC = 0.0;
 	
 	// check if all SW applied to surface or distributed throught subsurface
@@ -556,10 +563,10 @@
 	// SURFACE ROUGHNESS (Bougamont, 2006)
 	// wind/temperature surface roughness height [m]
-	if (ds < 910 && Ws == 0) z0 = 0.00012;       // 0.12 mm for dry snow
-	else if (ds >= 910) z0 = 0.0032;             // 3.2 mm for ice
+	if (ds < dIce-Dtol && fabs(Ws) < Wtol) z0 = 0.00012;       // 0.12 mm for dry snow
+	else if (ds >= dIce-Dtol) z0 = 0.0032;             // 3.2 mm for ice
 	else z0 = 0.0013;                            // 1.3 mm for wet snow
 
 	// if V = 0 goes to infinity therfore if V = 0 change
-	if(V<.01)V=.01;
+	if(V<0.01+Dtol)V=0.01;
 	
 	// Bulk-transfer coefficient for turbulent fluxes
@@ -574,8 +581,8 @@
 
 	// for snow and firn (density < 910 kg m-3) (Sturn et al, 1997)
-	for(int i=0;i<m;i++) if(d[i]<910) K[i] = 0.138 - 1.01E-3 * d[i] + 3.233E-6 * (pow(d[i],2));
+	for(int i=0;i<m;i++) if(d[i]<dIce-Dtol) K[i] = 0.138 - 1.01E-3 * d[i] + 3.233E-6 * (pow(d[i],2));
 
 	// for ice (density >= 910 kg m-3)
-	for(int i=0;i<m;i++) if(d[i]>=910) K[i] = 9.828 * exp(-5.7E-3*T[i]);
+	for(int i=0;i<m;i++) if(d[i]>=dIce-Dtol) K[i] = 9.828 * exp(-5.7E-3*T[i]);
 
 	// THERMAL DIFFUSION COEFFICIENTS
@@ -632,5 +639,5 @@
 	max_fdt=f[0];
 	for(int i=0;i<45;i++){
-		if (f[i]<dt)if(f[i]>=max_fdt)max_fdt=f[i];
+		if (f[i]<dt-Dtol)if(f[i]>=max_fdt-Dtol)max_fdt=f[i];
 	}
 	dt=max_fdt;
@@ -641,6 +648,6 @@
 	Ap = xNew<IssmDouble>(m);
 	for(int i=0;i<m;i++){
-		Au[i] = pow((dzU[i]/2/KP[i] + dz[i]/2/KU[i]),-1);
-		Ad[i] = pow((dzD[i]/2/KP[i] + dz[i]/2/KD[i]),-1);
+		Au[i] = pow((dzU[i]/2.0/KP[i] + dz[i]/2.0/KU[i]),-1.0);
+		Ad[i] = pow((dzD[i]/2.0/KP[i] + dz[i]/2.0/KD[i]),-1.0);
 		Ap[i] = (d[i]*dz[i]*CI)/dt;
 	}
@@ -653,17 +660,17 @@
 		Nu[i] = Au[i] / Ap[i];
 		Nd[i] = Ad[i] / Ap[i];
-		Np[i]= 1 - Nu[i] - Nd[i];
+		Np[i]= 1.0 - Nu[i] - Nd[i];
 	}
 	
 	// specify boundary conditions: constant flux at bottom
-	Nu[m-1] = 0;
-	Np[m-1] = 1;
+	Nu[m-1] = 0.0;
+	Np[m-1] = 1.0;
 	
 	// zero flux at surface
-	Np[0] = 1 - Nd[0];
+	Np[0] = 1.0 - Nd[0];
 	
 	// Create neighbor arrays for diffusion calculations instead of a tridiagonal matrix
-	Nu[0] = 0;
-	Nd[m-1] = 0;
+	Nu[0] = 0.0;
+	Nd[m-1] = 0.0;
 	
 	/* RADIATIVE FLUXES*/
@@ -704,10 +711,10 @@
 		// calculated.  The estimated enegy balance & melt are significanly
 		// less when Ts is taken as the mean of the x top grid cells.
-		Ts = (T[0] + T[1])/2;
-		Ts = fmin(273.15,Ts);    // don't allow Ts to exceed 273.15 K (0°C)
+		Ts = (T[0] + T[1])/2.0;
+		Ts = fmin(CtoK,Ts);    // don't allow Ts to exceed 273.15 K (0 degC)
 		
 		//TURBULENT HEAT FLUX
     
-		// Monin–Obukhov Stability Correction
+		// Monin-Obukhov Stability Correction
 		// Reference:
 		// Ohmura, A., 1982: Climate and Energy-Balance on the Arctic Tundra.
@@ -715,7 +722,7 @@
 
 		// calculate the Bulk Richardson Number (Ri)
-		Ri = (2*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2.0));
+		Ri = (2.0*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2.0));
 		
-		// calculate Monin-–Obukhov stability factors 'coefM' and 'coefH'
+		// calculate Monin-Obukhov stability factors 'coefM' and 'coefH'
     
 		// do not allow Ri to exceed 0.19
@@ -723,19 +730,19 @@
 
 		// calculate momentum 'coefM' stability factor
-		if (Ri > 0){
+		if (Ri > 0.0+Ttol){
 			// if stable
-			coefM = 1/(1-5.2*Ri);
+			coefM = 1.0/(1.0-5.2*Ri);
 		}
 		else {
-			coefM =pow (1-18*Ri,-0.25);
+			coefM =pow (1.0-18*Ri,-0.25);
 		}
 		
 		// calculate heat/wind 'coef_H' stability factor
-		if (Ri < -0.03) coefH = 1.3 * coefM;
+		if (Ri < -0.03-Ttol) coefH = 1.3 * coefM;
 		else coefH = coefM;
 		
 		//// Sensible Heat
 		// calculate the sensible heat flux [W m-2](Patterson, 1998)
-		shf = C * 1005 * (Ta - Ts);
+		shf = C * CA * (Ta - Ts);
 
 		// adjust using Monin–Obukhov stability theory
@@ -744,6 +751,6 @@
 		//// Latent Heat
 		// determine if snow pack is melting & calcualte surface vapour pressure over ice or liquid water
-		if (Ts >= 273.15){
-			L = 2.495E6;
+		if (Ts >= CtoK-Ttol){
+			L = LV;
 
 			// for an ice surface Murphy and Koop, 2005 [Equation 7]
@@ -751,7 +758,7 @@
 		}
 		else{
-			L = 2.8295E6; // latent heat of sublimation for liquid surface (assume liquid on surface when Ts == 0 deg C)
+			L = LS; // latent heat of sublimation for liquid surface (assume liquid on surface when Ts == 0 deg C)
 			// Wright (1997), US Meteorological Handbook from Murphy and Koop, 2005 Appendix A
-			eS = 611.21 * exp(17.502 * (Ts - 273.15) / (240.97 + Ts - 273.15));
+			eS = 611.21 * exp(17.502 * (Ts - CtoK) / (240.97 + Ts - CtoK));
 		}
 		
@@ -759,10 +766,10 @@
 		lhf = C * L * (eAir - eS) * 0.622 / pAir;
 		
-		// adjust using Monin–Obukhov stability theory (if lhf '+' then there is energy and mass gained at the surface, 
+		// adjust using Monin-Obukhov stability theory (if lhf '+' then there is energy and mass gained at the surface, 
 		// if '-' then there is mass and energy loss at the surface.
 		lhf = lhf / (coefM * coefH);
 
 		//mass loss (-)/acreation(+) due to evaporation/condensation [kg]
-		EC_day = lhf * 86400 / L;
+		EC_day = lhf * dts / L;
 
 		// temperature change due turbulent fluxes
@@ -787,5 +794,5 @@
 		
 		// calculate cumulative evaporation (+)/condensation(-)
-		EC = EC + (EC_day/86400)*dt;
+		EC = EC + (EC_day/dts)*dt;
     
 		/* CHECK FOR ENERGY (E) CONSERVATION [UNITS: J]
@@ -819,5 +826,4 @@
 	xDelete<IssmDouble>(sw);
 	xDelete<IssmDouble>(dT_sw);
-	xDelete<IssmDouble>(lw);
 	xDelete<IssmDouble>(T0);
 	xDelete<IssmDouble>(Tu);
@@ -829,5 +835,5 @@
 
 }  /*}}}*/
-void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, int m, int sid){ /*{{{*/
+void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, IssmDouble dIce, int m, int sid){ /*{{{*/
 
 	// DISTRIBUTES ABSORBED SHORTWAVE RADIATION WITHIN SNOW/ICE
@@ -865,5 +871,5 @@
         
 		// calculate surface shortwave radiation fluxes [W m-2]
-		swf[0] = (1 - as) * dsw;
+		swf[0] = (1.0 - as) * dsw;
 	}
 	else{ // sw radation is absorbed at depth within the glacier
@@ -884,5 +890,5 @@
 			// convert effective radius [mm] to grain size [m]
 			gsz=xNew<IssmDouble>(m);
-			for(int i=0;i<m;i++) gsz[i]= (re[i] * 2) / 1000;
+			for(int i=0;i<m;i++) gsz[i]= (re[i] * 2.0) / 1000.0;
 
 			// Spectral fractions [0.3-0.8um 0.8-1.5um 1.5-2.8um]
@@ -894,6 +900,6 @@
 			B2_cum=xNew<IssmDouble>(m+1);
 			for(int i=0;i<m+1;i++){
-				B1_cum[i]=1;
-				B2_cum[i]=1;
+				B1_cum[i]=1.0;
+				B2_cum[i]=1.0;
 			}
 
@@ -901,7 +907,7 @@
 			// spectral albedos:
 			// 0.3 - 0.8um
-			IssmDouble a0 = fmin(0.98, 1 - 1.58 *pow(gsz[0],0.5));
+			IssmDouble a0 = fmin(0.98, 1.0 - 1.58 *pow(gsz[0],0.5));
 			// 0.8 - 1.5um
-			IssmDouble a1 = fmax(0, 0.95 - 15.4 *pow(gsz[0],0.5));
+			IssmDouble a1 = fmax(0.0, 0.95 - 15.4 *pow(gsz[0],0.5));
 			// 1.5 - 2.8um
 			IssmDouble a2 = fmax(0.127, 0.88 + 346.3*gsz[0] - 32.31*pow(gsz[0],0.5));
@@ -909,7 +915,7 @@
 			// separate net shortwave radiative flux into spectral ranges
 			IssmDouble swfS[3];
-			swfS[0] = (sF[0] * dsw) * (1 - a0);
-			swfS[1] = (sF[1] * dsw) * (1 - a1);
-			swfS[2] = (sF[2] * dsw) * (1 - a2);
+			swfS[0] = (sF[0] * dsw) * (1.0 - a0);
+			swfS[1] = (sF[1] * dsw) * (1.0 - a1);
+			swfS[2] = (sF[2] * dsw) * (1.0 - a2);
 
 			// absorption coeficient for spectral range:
@@ -918,6 +924,6 @@
 			B2 =xNew<IssmDouble>(m);
 			for(int i=0;i<m;i++) h[i]= d[i] /(pow(gsz[i],0.5));
-			for(int i=0;i<m;i++) B1[i] = .0192 * h[i];                 // 0.3 - 0.8um
-			for(int i=0;i<m;i++) B2[i]= .1098 * h[i];                 // 0.8 - 1.5um
+			for(int i=0;i<m;i++) B1[i] = 0.0192 * h[i];                 // 0.3 - 0.8um
+			for(int i=0;i<m;i++) B2[i]= 0.1098 * h[i];                 // 0.8 - 1.5um
 			// B3 = +inf                     // 1.5 - 2.8um
 
@@ -986,16 +992,16 @@
 
 			// calculate surface shortwave radiation fluxes [W m-2]
-			IssmDouble swf_s = SWs * (1 - as) * dsw;
+			IssmDouble swf_s = SWs * (1.0 - as) * dsw;
 
 			// calculate surface shortwave radiation fluxes [W m-2]
-			IssmDouble swf_ss = (1-SWs) * (1 - as) * dsw;
+			IssmDouble swf_ss = (1.0-SWs) * (1.0 - as) * dsw;
 
 			// SW allowed to penetrate into snowpack
-			IssmDouble Bs = 10;    // snow SW extinction coefficient [m-1] (Bassford,2006)
+			IssmDouble Bs = 10.0;    // snow SW extinction coefficient [m-1] (Bassford,2006)
 			IssmDouble Bi = 1.3;   // ice SW extinction coefficient [m-1] (Bassford,2006)
 
 			// calculate extinction coefficient B [m-1] vector
 			B=xNew<IssmDouble>(m);
-			for(int i=0;i<m;i++) B[i] = Bs + (300 - d[i]) * ((Bs - Bi)/(910 - 300));
+			for(int i=0;i<m;i++) B[i] = Bs + (300.0 - d[i]) * ((Bs - Bi)/(dIce - 300.0));
 
 			// cumulative extinction factor
@@ -1004,5 +1010,5 @@
 			for(int i=0;i<m;i++)exp_B[i]=exp(-B[i]*dz[i]);
 
-			B_cum[0]=1;
+			B_cum[0]=1.0;
 			for(int i=0;i<m;i++){
 				IssmDouble cum_B=exp_B[0];
@@ -1032,5 +1038,5 @@
 
 } /*}}}*/ 
-void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, IssmDouble T_air, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, int sid){ /*{{{*/
+void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, IssmDouble T_air, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, IssmDouble dIce, int sid){ /*{{{*/
 
 	// Adds precipitation and deposition to the model grid
@@ -1051,8 +1057,7 @@
 	// MAIN FUNCTION
 	// specify constants
-	const IssmDouble dIce = 910;     // density of ice [kg m-3]
 	const IssmDouble dSnow = 150;    // density of snow [kg m-3]
 	const IssmDouble reNew = 0.1;    // new snow grain size [mm]
-	const IssmDouble gdnNew = 1;     // new snow dendricity 
+	const IssmDouble gdnNew = 1.0;     // new snow dendricity 
 	const IssmDouble gspNew = 0.5;   // new snow sphericity 
 
@@ -1060,5 +1065,7 @@
 	IssmDouble* mInit=NULL;
 	bool        top=true;
-	IssmDouble  mass, massinit, mass_diff;
+	IssmDouble  mass=0;
+	IssmDouble  massinit=0;
+	IssmDouble  mass_diff=0;
 
 	/*output: */
@@ -1071,5 +1078,5 @@
 	IssmDouble* gdn=NULL;
 	IssmDouble* gsp=NULL;
-	int         m;
+	int         m=0;
 
 	if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   accumulation module\n");
@@ -1089,20 +1096,21 @@
 	mInit=xNew<IssmDouble>(m);
 	for(int i=0;i<m;i++) mInit[i]= d[i] * dz[i];
-	massinit=0; for(int i=0;i<m;i++)massinit+=mInit[i];
-
-	if (P > 0){
+	massinit=0.0; 
+	for(int i=0;i<m;i++)massinit+=mInit[i];
+
+	if (P > 0.0+Dtol){
 			
 
-		if (T_air <= 273.15){ // if snow
+		if (T_air <= CtoK+Ttol){ // if snow
 
 			IssmDouble  z_snow = P/dSnow;               // depth of snow
 
 			// if snow depth is greater than specified min dz, new cell created
-			if (z_snow > dzMin){
+			if (z_snow > dzMin+Dtol){
 
 				newcell(&T,T_air,top,m); //new cell T
 				newcell(&dz,z_snow,top,m); //new cell dz
 				newcell(&d,dSnow,top,m); //new cell d
-				newcell(&W,0,top,m); //new cell W
+				newcell(&W,0.0,top,m); //new cell W
 				newcell(&a,aSnow,top,m); //new cell a
 				newcell(&re,reNew,top,m); //new cell grain size
@@ -1137,7 +1145,4 @@
 			  makes the numerics easier.*/
 
-			IssmDouble LF = 0.3345E6;  // latent heat of fusion(J kg-1)
-			IssmDouble CI = 2102;      // specific heat capacity of snow/ice (J kg-1 k-1)
-
 			// grid cell adjust mass
 			mass = mInit[0] + P;
@@ -1151,5 +1156,5 @@
 
 			// if d > the density of ice, d = dIce
-			if (d[0] > dIce){
+			if (d[0] > dIce+Dtol){
 				d[0] = dIce;           // adjust d
 				dz[0] = mass / d[0];    // dz is adjusted to conserve mass
@@ -1163,5 +1168,5 @@
 		
 		#ifndef _HAVE_ADOLC_  //avoid round operation. only check in forward mode.
-		mass_diff = round(mass_diff * 100)/100;
+		mass_diff = round(mass_diff * 100.0)/100.0;
 		if (mass_diff > 0) _error_("mass not conserved in accumulation function");
 		#endif
@@ -1182,5 +1187,5 @@
 	*pm=m;
 } /*}}}*/
-void melt(IssmDouble* pM, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, int sid){ /*{{{*/
+void melt(IssmDouble* pM, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, IssmDouble dIce, int sid){ /*{{{*/
 
 	//// MELT ROUTINE
@@ -1200,5 +1205,5 @@
 	IssmDouble* F=NULL;
 	IssmDouble* flxDn=NULL;
-	IssmDouble  ER=0;
+	IssmDouble  ER=0.0;
 	IssmDouble* EI=NULL;
 	IssmDouble* EW=NULL;
@@ -1206,28 +1211,28 @@
 	int*       D=NULL;
 	
-	IssmDouble sumM;
-	IssmDouble sumER;
-	IssmDouble addE;
-	IssmDouble mSum0;
-	IssmDouble sumE0;
-	IssmDouble mSum1;
-	IssmDouble sumE1;
-	IssmDouble dE;
-	IssmDouble dm;
-	IssmDouble X;
-	IssmDouble Wi;
+	IssmDouble sumM=0.0;
+	IssmDouble sumER=0.0;
+	IssmDouble addE=0.0;
+	IssmDouble mSum0=0.0;
+	IssmDouble sumE0=0.0;
+	IssmDouble mSum1=0.0;
+	IssmDouble sumE1=0.0;
+	IssmDouble dE=0.0;
+	IssmDouble dm=0.0;
+	IssmDouble X=0.0;
+	IssmDouble Wi=0.0;
     
-	IssmDouble Ztot;
-	IssmDouble T_bot;
-	IssmDouble m_bot;
-	IssmDouble dz_bot;
-	IssmDouble d_bot;
-	IssmDouble W_bot;
-	IssmDouble a_bot;
-	IssmDouble re_bot;
-	IssmDouble gdn_bot;
-	IssmDouble gsp_bot;
-	IssmDouble EI_bot;
-	IssmDouble EW_bot;
+	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 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;
     
@@ -1235,14 +1240,14 @@
 	IssmDouble* dzMin2=NULL;
 	IssmDouble zY2=1.025;
-	bool lastCellFlag;
+	bool lastCellFlag = false;
 	int X1=0;
 	int X2=0;
     
-	int        D_size;
+	int        D_size = 0;
 
 	/*outputs:*/
-	IssmDouble  mAdd;
-	IssmDouble dz_add;
-	IssmDouble  Rsum;
+	IssmDouble  mAdd = 0.0;
+	IssmDouble dz_add = 0.0;
+	IssmDouble  Rsum = 0.0;
 	IssmDouble* T=*pT;
 	IssmDouble* d=*pd;
@@ -1254,5 +1259,5 @@
 	IssmDouble* gsp=*pgsp;
 	int         n=*pn;
-	IssmDouble* R=0;
+	IssmDouble* R=NULL;
 	
 	if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   melt module\n");
@@ -1265,11 +1270,4 @@
 	dW=xNew<IssmDouble>(n); 
 
-	// specify constants
-	const IssmDouble CtoK = 273.15;  // clecius to Kelvin conversion
-	const IssmDouble CI = 2102;      // specific heat capacity of snow/ice (J kg-1 k-1)
-	const IssmDouble LF = 0.3345E6;  // latent heat of fusion(J kg-1)
-	const IssmDouble dPHC = 830;     // pore hole close off density[kg m-3]
-	const IssmDouble dIce = 910;     // density of ice [kg m-3]
-
 	// store initial mass [kg] and energy [J]
 	m=xNew<IssmDouble>(n); for(int i=0;i<n;i++) m[i] = dz[i]* d[i];                    // grid cell mass [kg]
@@ -1281,17 +1279,18 @@
 
 	// initialize melt and runoff scalars
-	Rsum = 0;       // runoff [kg]
-	sumM = 0;       // total melt [kg]
-	mAdd = 0;       // mass added/removed to/from base of model [kg]
-	addE = 0;       // energy added/removed to/from base of model [J]
-	dz_add=0;      // thickness of the layer added/removed to/from base of model [m]
+	Rsum = 0.0;       // runoff [kg]
+	sumM = 0.0;       // total melt [kg]
+	mAdd = 0.0;       // mass added/removed to/from base of model [kg]
+	addE = 0.0;       // energy added/removed to/from base of model [J]
+	dz_add=0.0;      // thickness of the layer added/removed to/from base of model [m]
 
 	// calculate temperature excess above 0 deg C
 	exsT=xNewZeroInit<IssmDouble>(n);
-	for(int i=0;i<n;i++) exsT[i]= fmax(0, T[i] - CtoK);        // [K] to [°C]
+	for(int i=0;i<n;i++) exsT[i]= fmax(0.0, T[i] - CtoK);        // [K] to [degC]
 
 	// new grid point center temperature, T [K]
-	for(int i=0;i<n;i++) T[i]-=exsT[i];
-
+	//for(int i=0;i<n;i++) T[i]-=exsT[i];
+	for(int i=0;i<n;i++) T[i]=fmin(T[i],CtoK);
+	
 	// specify irreducible water content saturation [fraction]
 	const IssmDouble Swi = 0.07;                     // assumed constant after Colbeck, 1974
@@ -1299,8 +1298,8 @@
 	//// REFREEZE PORE WATER
 	// check if any pore water
-	if (cellsum(W,n) > 0){
+	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]
-		for(int i=0;i<n;i++) maxF[i] = fmax(0, -((T[i] - CtoK) * m[i] * CI) / LF);
+		for(int i=0;i<n;i++) maxF[i] = fmax(0.0, -((T[i] - CtoK) * m[i] * CI) / LF);
 
 		// freeze pore water and change snow/ice properties
@@ -1319,6 +1318,6 @@
 	exsW=xNew<IssmDouble>(n); 
 	for(int i=0;i<n;i++){
-		Wi= (910 - d[i]) * Swi * (m[i] / d[i]);        // irreducible water content [kg]
-		exsW[i] = fmax(0, W[i] - Wi);                  // water "squeezed" from snow [kg]
+		Wi= (dIce - d[i]) * Swi * (m[i] / d[i]);        // irreducible water content [kg]
+		exsW[i] = fmax(0.0, W[i] - Wi);                  // water "squeezed" from snow [kg]
 	}
 
@@ -1326,5 +1325,5 @@
 
 	// run melt algorithm if there is melt water or excess pore water
-	if ((cellsum(exsT,n) > 0) || (cellsum(exsW,n) > 0)){
+	if ((cellsum(exsT,n) > 0.0 + Ttol) || (cellsum(exsW,n) > 0.0 + Wtol)){
 		// _printf_(""MELT OCCURS");
 		// check to see if thermal energy exceeds energy to melt entire cell
@@ -1332,7 +1331,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, exsT [i]- 159.1342);
-        
-		if (cellsum(surpT,n) > 0 ){
+		surpT=xNew<IssmDouble>(n); for(int i=0;i<n;i++)surpT[i] = fmax(0.0, exsT[i]- 159.1342);
+        
+		if (cellsum(surpT,n) > 0.0 + Ttol ){
 			// _printf_("T Surplus");
 			// calculate surplus energy
@@ -1340,17 +1339,17 @@
             
 			int i = 0;
-			while (cellsum(surpE,n) > 0){
+			while (cellsum(surpE,n) > 0.0 + Ttol){
 				// use surplus energy to increase the temperature of lower cell
 				T[i+1] = surpE[i]/m[i+1]/CI + T[i+1];
                 
-				exsT[i+1] = fmax(0, T[i+1] - CtoK) + exsT[i+1];
+				exsT[i+1] = fmax(0.0, T[i+1] - CtoK) + exsT[i+1];
 				T[i+1] = fmin(CtoK, T[i+1]);
                 
-				surpT[i+1] = fmax(0, exsT[i+1] - 159.1342);
+				surpT[i+1] = fmax(0.0, exsT[i+1] - 159.1342);
 				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;
-				surpE[i] = 0;
+				surpE[i] = 0.0;
 				i = i + 1;
 			}
@@ -1362,5 +1361,5 @@
 
 		// calculate maximum refreeze amount, maxF [kg]
-		for(int i=0;i<n;i++)maxF[i] = fmax(0, -((T[i] - CtoK) * d[i] * dz[i] * CI)/ LF);
+		for(int i=0;i<n;i++)maxF[i] = fmax(0.0, -((T[i] - CtoK) * d[i] * dz[i] * CI)/ LF);
 
 		// initialize refreeze, runoff, flxDn and dW vectors [kg]
@@ -1368,5 +1367,5 @@
 		IssmDouble* R=xNewZeroInit<IssmDouble>(n);
 
-		for(int i=0;i<n;i++)dW[i] = 0;
+		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];
 
@@ -1374,5 +1373,5 @@
 		X = 0;
 		for(int i=n-1;i>=0;i--){
-			if(M[i]>0 || reCast<int,IssmDouble>(exsW[i])){
+			if(M[i]>0.0+Wtol || exsW[i]>0.0+Wtol){
 				X=i;
 				break;
@@ -1386,10 +1385,10 @@
 
 			// break loop if there is no meltwater and if depth is > mw_depth
-			if (inM == 0 && i > X){
+			if (inM == 0.0 && i > X){
 				break;
 			}
 
 			// if reaches impermeable ice layer all liquid water runs off (R)
-			else if (d[i] >= dIce){  // dPHC = pore hole close off [kg m-3]
+			else if (d[i] >= dIce-Dtol){  // dPHC = pore hole close off [kg m-3]
 				// _printf_("ICE LAYER");
 				// no water freezes in this cell
@@ -1398,10 +1397,10 @@
 
 				m[i] = m[i] - M[i];                     // mass after melt
-				Wi = (910-d[i]) * Swi * (m[i]/d[i]);    // irreducible water 
+				Wi = (dIce-d[i]) * Swi * (m[i]/d[i]);    // irreducible water 
 				dW[i] = fmin(inM, Wi - W[i]);            // change in pore water
-				R[i] = fmax(0, inM - dW[i]);             // runoff
+				R[i] = fmax(0.0, inM - dW[i]);             // runoff
 			}
 			// check if no energy to refreeze meltwater
-			else if (maxF[i] == 0){
+			else if (fabs(maxF[i]) < Dtol){
 				// _printf_("REFREEZE == 0");
 				// no water freezes in this cell
@@ -1409,8 +1408,8 @@
 
 				m[i] = m[i] - M[i];                     // mass after melt
-				Wi = (910-d[i]) * Swi * (m[i]/d[i]);    // irreducible water 
+				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, inM-dW[i]);         // meltwater out
-				F[i] = 0;                               // no freeze 
+				flxDn[i+1] = fmax(0.0, inM-dW[i]);         // meltwater out
+				F[i] = 0.0;                               // no freeze 
 			}
 			// some or all meltwater refreezes
@@ -1426,15 +1425,15 @@
 
 				//-----------------------pore water-----------------------------
-				Wi = (910-d[i])* Swi * dz_0;            // irreducible water 
+				Wi = (dIce-d[i])* Swi * dz_0;            // irreducible water 
 				dW[i] = fmin(inM - F1, Wi-W[i]);         // change in pore water
-				if (-dW[i]>W[i] ){
-					dW[i]= W[i];
+				if (-1*dW[i]>W[i]-Wtol ){
+					dW[i]= -1*W[i];
 				}
-				IssmDouble F2 = 0;                                 
-
-				if (dW[i] < 0){                         // excess pore water
+				IssmDouble F2 = 0.0;                                 
+
+				if (dW[i] < 0.0-Wtol){                         // excess pore water
 					dMax = (dIce - d[i])*dz_0;          // maximum refreeze                                             
 					IssmDouble maxF2 = fmin(dMax, maxF[i]-F1);      // maximum refreeze
-					F2 = fmin(-dW[i], maxF2);            // pore water refreeze
+					F2 = fmin(-1*dW[i], maxF2);            // pore water refreeze
 					m[i] = m[i] + F2;                   // mass after refreeze
 					d[i] = m[i]/dz_0;
@@ -1446,10 +1445,10 @@
 
 				// check if an ice layer forms 
-				if (d[i] == dIce){
+				if (fabs(d[i] - dIce) < Dtol){
 					// _printf_("ICE LAYER FORMS");
 					// excess water runs off
 					R[i] = flxDn[i+1];
 					// no water percolates to lower cell
-					flxDn[i+1] = 0;
+					flxDn[i+1] = 0.0;
 				}
 			}
@@ -1458,5 +1457,5 @@
 
 		//// GRID CELL SPACING AND MODEL DEPTH
-		for(int i=0;i<n;i++)if (W[i] < 0) _error_("negative pore water generated in melt equations");
+		for(int i=0;i<n;i++)if (W[i] < 0.0 - Wtol) _error_("negative pore water generated in melt equations");
 		
 		// delete all cells with zero mass
@@ -1465,5 +1464,6 @@
 
 		// delete all cells with zero mass
-		D_size=0; for(int i=0;i<n;i++)if(m[i]!=0)D_size++; D=xNew<int>(D_size);
+		D_size=0; for(int i=0;i<n;i++)if(m[i]!=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++;}
 		
@@ -1504,5 +1504,5 @@
     
 	for (int i=0;i<n;i++){
-		if (Zcum[i]<=zTop){
+		if (Zcum[i]<=zTop+Dtol){
 			dzMin2[i]=dzMin;
 		}
@@ -1515,5 +1515,5 @@
 	X = 0;
 	for(int i=n-1;i>=0;i--){
-		if(dz[i]<dzMin2[i]){
+		if(dz[i]<dzMin2[i]-Dtol){
 			X=i;
 			break;
@@ -1531,5 +1531,5 @@
 
 	for (int i = 0; i<=X;i++){
-		if (dz [i] < dzMin2[i]){                              // merge top cells
+		if (dz[i] < dzMin2[i]-Dtol){                              // merge top cells
 			// _printf_("dz > dzMin * 2')
 			// adjust variables as a linearly weighted function of mass
@@ -1542,5 +1542,5 @@
             
 			// merge with underlying grid cell and delete old cell
-			dz [i+1] = dz[i] + dz[i+1];                 // combine cell depths
+			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
@@ -1572,5 +1572,5 @@
         
 		// merge with underlying grid cell and delete old cell
-  		dz [X1] = dz[X2] + dz[X1];                 // combine cell depths
+  		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
@@ -1582,5 +1582,6 @@
 
 	// delete combined cells
-	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_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++;}
 
@@ -1602,5 +1603,5 @@
 	X=0;
 	for(int i=9;i>=0;i--){
-		if(dz[i]> 2* dzMin){
+		if(dz[i]> 2.0*dzMin+Dtol){
 			X=i;
 			break;
@@ -1608,25 +1609,25 @@
 	}
 	
-	int i=0;
-	while(i<=X){
-		if (dz [i] > dzMin *2){
+	int j=0;
+	while(j<=X){
+		if (dz[j] > dzMin*2.0+Dtol){
 
 				// _printf_("dz > dzMin * 2");
 				// split in two
-				cellsplit(&dz, n, i,.5);
-				cellsplit(&W, n, i,.5);
-				cellsplit(&m, n, i,.5);
-				cellsplit(&T, n, i,1.0);
-				cellsplit(&d, n, i,1.0);
-				cellsplit(&a, n, i,1.0);
-				cellsplit(&EI, n, i,1.0);
-				cellsplit(&EW, n, i,1.0);
-				cellsplit(&re, n, i,1.0);
-				cellsplit(&gdn, n, i,1.0);
-				cellsplit(&gsp, n, i,1.0);
+				cellsplit(&dz, n, j,.5);
+				cellsplit(&W, n, j,.5);
+				cellsplit(&m, n, j,.5);
+				cellsplit(&T, n, j,1.0);
+				cellsplit(&d, n, j,1.0);
+				cellsplit(&a, n, j,1.0);
+				cellsplit(&EI, n, j,1.0);
+				cellsplit(&EW, n, j,1.0);
+				cellsplit(&re, n, j,1.0);
+				cellsplit(&gdn, n, j,1.0);
+				cellsplit(&gsp, n, j,1.0);
 				n++;
 				X=X+1;
 		}
-		else i++;
+		else j++;
 	}
 
@@ -1637,5 +1638,5 @@
 	Ztot = cellsum(dz,n);
     
-	if (Ztot < zMin){
+	if (Ztot < zMin-Dtol){
 		// printf("Total depth < zMin %f \n", Ztot);
 		// mass and energy to be added
@@ -1671,5 +1672,5 @@
 		n=n+1;
 	}
-	else if (Ztot > zMax){
+	else if (Ztot > zMax+Dtol){
 		// printf("Total depth > zMax %f \n", Ztot);
 		// mass and energy loss
@@ -1710,5 +1711,5 @@
 
 	/*checks: */
-	for(int i=0;i<n;i++) if (W[i]<0) _error_("negative pore water generated in melt equations\n");
+	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: */
@@ -1753,5 +1754,5 @@
 
 } /*}}}*/ 
-void densification(IssmDouble* d,IssmDouble* dz, IssmDouble* T, IssmDouble* re, int denIdx, IssmDouble C, IssmDouble dt, IssmDouble Tmean,IssmDouble dIce, int m, int sid){ /*{{{*/
+void densification(IssmDouble* d,IssmDouble* dz, IssmDouble* T, IssmDouble* re, int denIdx, IssmDouble C, IssmDouble dt, IssmDouble Tmean, IssmDouble dIce, int m, int sid){ /*{{{*/
 
 	//// THIS NEEDS TO BE DOUBLE CHECKED AS THERE SEAMS TO BE LITTLE DENSIFICATION IN THE MODEL OUTOUT [MAYBE COMPATION IS COMPNSATED FOR BY TRACES OF SNOW???]
@@ -1795,5 +1796,5 @@
 	//// MAIN FUNCTION
 	// specify constants
-	dt      = dt / 86400;  // convert from [s] to [d]
+	dt      = dt / dts;  // convert from [s] to [d]
 	// R     = 8.314        // gas constant [mol-1 K-1]
 	// Ec    = 60           // activation energy for self-diffusion of water
@@ -1802,5 +1803,7 @@
 
 	/*intermediary: */
-	IssmDouble c0,c1,H;
+	IssmDouble c0=0.0;
+	IssmDouble c1=0.0;
+	IssmDouble H=0.0;
 	
 	if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   densification module\n");
@@ -1815,5 +1818,5 @@
 
 	IssmDouble* obp = xNew<IssmDouble>(m);
-	obp[0]=0;
+	obp[0]=0.0;
 	for(int i=1;i<m;i++)obp[i]=cumdz[i-1]*d[i-1];
 	
@@ -1826,11 +1829,11 @@
 		switch (denIdx){
 			case 1: // Herron and Langway (1980)
-				c0 = (11 * exp(-10160 / (T[i] * 8.314))) * C/1000;
-				c1 = (575 * exp(-21400 / (T[i]* 8.314))) * pow(C/1000,.5);
+				c0 = (11.0 * exp(-10160.0 / (T[i] * 8.314))) * C/1000.0;
+				c1 = (575.0 * exp(-21400.0 / (T[i]* 8.314))) * pow(C/1000.0,.5);
 				break;
 			case 2: // Arthern et al. (2010) [semi-emperical]
 				// common variable
 				// NOTE: Ec=60000, Eg=42400 (i.e. should be in J not kJ)
-				H = exp((-60000./(T[i] * 8.314)) + (42400./(Tmean * 8.314))) * (C * 9.81);
+				H = exp((-60000.0/(T[i] * 8.314)) + (42400.0/(Tmean * 8.314))) * (C * 9.81);
 				c0 = 0.07 * H;
 				c1 = 0.03 * H;
@@ -1840,5 +1843,5 @@
 
 				// common variable
-				H = exp((-60/(T[i] * 8.314))) * obp[i] / pow(re[i]/1000,2.0);
+				H = exp((-60.0/(T[i] * 8.314))) * obp[i] / pow(re[i]/1000.0,2.0);
 				c0 = 9.2e-9 * H;
 				c1 = 3.7e-9 * H;
@@ -1846,5 +1849,5 @@
 
 			case 4: // Li and Zwally (2004)
-				c0 = (C/dIce) * (139.21 - 0.542*Tmean)*8.36*pow(273.15 - T[i],-2.061);
+				c0 = (C/dIce) * (139.21 - 0.542*Tmean)*8.36*pow(CtoK - T[i],-2.061);
 				c1 = c0;
 				break;
@@ -1852,5 +1855,5 @@
 			case 5: // Helsen et al. (2008)
 				// common variable
-				c0 = (C/dIce) * (76.138 - 0.28965*Tmean)*8.36*pow(273.15 - T[i],-2.061);
+				c0 = (C/dIce) * (76.138 - 0.28965*Tmean)*8.36*pow(CtoK - T[i],-2.061);
 				c1 = c0;
 				break;
@@ -1858,11 +1861,11 @@
 
 		// new snow density
-		if(d[i] <= 550) d[i] = d[i] + (c0 * (dIce - d[i]) / 365 * dt);
-		else            d[i] = d[i] + (c1 * (dIce - d[i]) / 365 * dt);
+		if(d[i] <= 550.0+Dtol) d[i] = d[i] + (c0 * (dIce - d[i]) / 365.0 * dt);
+		else            d[i] = d[i] + (c1 * (dIce - d[i]) / 365.0 * dt);
 
 		//disp((num2str(nanmean(c0 .* (dIce - d(idx)) / 365 * dt))))
 
 		// do not allow densities to exceed the density of ice
-		if(d[i]>dIce)d[i]=dIce;
+		if(d[i] > dIce-Dtol) d[i]=dIce;
 
 		// calculate new grid cell length
@@ -1876,5 +1879,5 @@
 
 } /*}}}*/
-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, 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){ /*{{{*/
 
 	//// TURBULENT HEAT FLUX
@@ -1899,35 +1902,33 @@
 	//   Tz: height above ground at which temperature (T) was sampled [m]
 
-	//// FUNCTION INITILIZATION 
-
-	// CA = 1005;                    // heat capacity of air (J kg-1 k-1)
-	// LF = 0.3345E6;                // latent heat of fusion(J kg-1)
-	// LV = 2.495E6;                 // latent heat of vaporization(J kg-1)
-	// dIce = 910;                   // density of ice [kg m-3]
-	
 	/*intermediary:*/
-	IssmDouble d_air;
-	IssmDouble Ri;
-	IssmDouble z0;
-	IssmDouble coef_M,coef_H;
-	IssmDouble An, C;
-	IssmDouble L, eS;
+	IssmDouble d_air=0.0;
+	IssmDouble Ri=0.0;
+	IssmDouble z0=0.0;
+	IssmDouble coef_M=0.0;
+	IssmDouble coef_H=0.0;
+	IssmDouble An=0.0;
+	IssmDouble C=0.0;
+	IssmDouble L=0.0;
+	IssmDouble eS=0.0;
 
 	/*output: */
-	IssmDouble shf, lhf, EC;
+	IssmDouble shf=0.0;
+	IssmDouble lhf=0.0;
+	IssmDouble EC=0.0;
 	
 	if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   turbulentFlux module\n");
 
 	// calculated air density [kg/m3]
-	d_air = 0.029 * pAir /(8.314 * Ta);
+	d_air = 0.029 * pAir /(R * Ta);
 
 	//// Determine Surface Roughness
 	// Bougamont, 2006
 	// wind/temperature surface roughness height [m]
-	if (ds < 910 && Ws == 0) z0 = 0.00012;               // 0.12 mm for dry snow
-	else if (ds >= 910) z0 = 0.0032;                // 3.2 mm for ice 
+	if (ds < dIce-Dtol && fabs(Ws) < Wtol) z0 = 0.00012;               // 0.12 mm for dry snow
+	else if (ds >= dIce-Dtol) z0 = 0.0032;                // 3.2 mm for ice 
 	else z0 = 0.0013;                // 1.3 mm for wet snow
 
-	//// Monin–Obukhov Stability Correction
+	//// Monin-Obukhov Stability Correction
 	// Reference:
 	// Ohmura, A., 1982: Climate and Energy-Balance on the Arctic Tundra.
@@ -1935,17 +1936,17 @@
 
 	// if V = 0 goes to infinity therfore if V = 0 change
-	if(V< .01) V=.01;
+	if(V < 0.01) V=0.01;
 
 	// calculate the Bulk Richardson Number (Ri)
-	Ri = (2*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2));
-
-	// calculate Monin–Obukhov stability factors 'coef_M' and 'coef_H'
+	Ri = (2.0*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2));
+
+	// calculate Monin-Obukhov stability factors 'coef_M' and 'coef_H'
 
 	// do not allow Ri to exceed 0.19
-	if(Ri>.19)Ri= 0.19;
+	if(Ri>0.19-Ttol)Ri= 0.19;
 
 	// calculate momentum 'coef_M' stability factor
-	if (Ri > 0) coef_M = pow(1-5.2*Ri,-1); // if stable
-	else coef_M = pow(1-18*Ri,-0.25);
+	if (Ri > 0.0) coef_M = pow(1.0-5.2*Ri,-1.0); // if stable
+	else coef_M = pow(1.0-18*Ri,-0.25);
 
 	// calculate heat/wind 'coef_H' stability factor
@@ -1959,7 +1960,7 @@
 	//// Sensible Heat
 	// calculate the sensible heat flux [W m-2](Patterson, 1998)
-	shf = C * 1005 * (Ta - Ts);
-
-	// adjust using Monin–Obukhov stability theory
+	shf = C * CA * (Ta - Ts);
+
+	// adjust using Monin-Obukhov stability theory
 	shf = shf / (coef_M * coef_H);
 
@@ -1967,6 +1968,6 @@
 	// determine if snow pack is melting & calcualte surface vapour pressure
 	// over ice or liquid water
-	if (Ts >= 273.15){
-		L = 2.495E6;
+	if (Ts >= CtoK-Ttol){
+		L = LV;
 
 		// for an ice surface Murphy and Koop, 2005 [Equation 7]
@@ -1974,9 +1975,9 @@
 	}
 	else{
-		L = 2.8295E6; // latent heat of sublimation
+		L = LS; // latent heat of sublimation
 		// for liquid surface (assume liquid on surface when Ts == 0 deg C)
 		// Wright (1997), US Meteorological Handbook from Murphy and Koop,
 		// 2005 Apendix A
-		eS = 611.21 * exp(17.502 * (Ts - 273.15) / (240.97 + Ts - 273.15));
+		eS = 611.21 * exp(17.502 * (Ts - CtoK) / (240.97 + Ts - CtoK));
 	}
 
@@ -1984,5 +1985,5 @@
 	lhf = C * L * (eAir - eS) * 0.622 / pAir;
 
-	// adjust using Monin–Obukhov stability theory (if lhf '+' then there is
+	// adjust using Monin-Obukhov stability theory (if lhf '+' then there is
 	// energy and mass gained at the surface, if '-' then there is mass and 
 	// energy loss at the surface. 
@@ -1990,5 +1991,5 @@
 
 	// mass loss (-)/acreation(+) due to evaporation/condensation [kg]
-	EC = lhf * 86400 / L;
+	EC = lhf * dts / L;
 
 	/*assign output poitners: */
Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp	(revision 22248)
+++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp	(revision 22249)
@@ -271,12 +271,12 @@
 	IssmDouble  h = -1.54e-5;
 	IssmDouble  smb,smbref,anomaly,yts,z;
-    
-    /* Get constants */
-    femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
-    /*iomodel->FindConstant(&yts,"md.constants.yts");*/
-    /*this->parameters->FindParam(&yts,ConstantsYtsEnum);*/
-    /*Mathieu original*/
-    /*IssmDouble  smb,smbref,z;*/
-    
+
+	/* Get constants */
+	femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
+	/*iomodel->FindConstant(&yts,"md.constants.yts");*/
+	/*this->parameters->FindParam(&yts,ConstantsYtsEnum);*/
+	/*Mathieu original*/
+	/*IssmDouble  smb,smbref,z;*/
+
 	/*Loop over all the elements of this partition*/
 	for(int i=0;i<femmodel->elements->Size();i++){
@@ -298,16 +298,16 @@
 			anomaly = smblistref[v];
 
-            /* Henning edited acc. to Riannes equations*/
-            /* Set SMB maximum elevation, if dz = 0 -> z_critical = 1675 */
-            z_critical = z_critical + dz;
-            
-            /* Calculate smb acc. to the surface elevation z */
-            if(z<z_critical){
+			/* Henning edited acc. to Riannes equations*/
+			/* Set SMB maximum elevation, if dz = 0 -> z_critical = 1675 */
+			z_critical = z_critical + dz;
+
+			/* Calculate smb acc. to the surface elevation z */
+			if(z<z_critical){
 				smb = a + b*z + c;
 			}
 			else{
-			  smb = (a + b*z)*(f + g*(z-z_critical) + h*(z-z_critical)*(z-z_critical)) + c;
-			}
-            
+				smb = (a + b*z)*(f + g*(z-z_critical) + h*(z-z_critical)*(z-z_critical)) + c;
+			}
+
 			/* Compute smb including anomaly,
 				correct for number of seconds in a year [s/yr]*/
Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h	(revision 22248)
+++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h	(revision 22249)
@@ -25,10 +25,10 @@
 IssmDouble Marbouty(IssmDouble T, IssmDouble d, IssmDouble dT);
 void grainGrowth(IssmDouble* pre, IssmDouble* pgdn, IssmDouble* pgsp, IssmDouble* T,IssmDouble* dz,IssmDouble* d, IssmDouble* W,IssmDouble smb_dt,int m,int aIdx, int sid);
-void albedo(IssmDouble* a,int aIdx, IssmDouble* re, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble* T, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt,int m, int sid);
-void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, int m, int sid);
-void thermo(IssmDouble* pEC, IssmDouble* T, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlw, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, int sid);
-void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, IssmDouble Ta, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, int sid); 
-void melt(IssmDouble* pM, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, int sid);
-void densification(IssmDouble* d,IssmDouble* dz, IssmDouble* T, IssmDouble* re, int denIdx, IssmDouble C, IssmDouble dt, IssmDouble Tmean,IssmDouble rho_ice,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, int sid);
+void albedo(IssmDouble* a,int aIdx, IssmDouble* re, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble* T, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m, int sid);
+void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, IssmDouble dIce, int m, int sid);
+void thermo(IssmDouble* pEC, IssmDouble* T, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlw, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble dIce, int sid);
+void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, IssmDouble Ta, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, IssmDouble dIce, int sid); 
+void melt(IssmDouble* pM, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, IssmDouble dIce, int sid);
+void densification(IssmDouble* d,IssmDouble* dz, IssmDouble* T, IssmDouble* re, int denIdx, 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);
 #endif  /* _SurfaceMassBalancex_H*/
Index: /issm/trunk-jpl/test/NightlyRun/test243.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test243.m	(revision 22248)
+++ /issm/trunk-jpl/test/NightlyRun/test243.m	(revision 22249)
@@ -45,5 +45,5 @@
 %Fields and tolerances to track changes
 field_names      ={'SmbDz','SmbT' ,'SmbD' ,'SmbRe','SmbGdn','SmbGsp','SmbA' ,'SmbEC','SmbMassBalance'};
-field_tolerances ={5e-4,5e-5,0.0006,0.0002,1e-5,0.0003,2e-5,2e-7,1e-7};
+field_tolerances ={1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12};
 
 field_values={...