Index: /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp	(revision 21231)
+++ /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp	(revision 21232)
@@ -58,16 +58,15 @@
 			iomodel->FetchDataToInput(elements,"md.smb.Tz",SmbTzEnum);
 			iomodel->FetchDataToInput(elements,"md.smb.Vz",SmbVzEnum);
-            iomodel->FetchDataToInput(elements,"md.smb.isInitialized",SmbIsInitializedEnum);
-            iomodel->FetchDataToInput(elements,"md.smb.isrestart",SmbIsrestartEnum);
-            iomodel->FetchDataToInput(elements,"md.smb.Dzrst",SmbDzrstEnum);
-            iomodel->FetchDataToInput(elements,"md.smb.Drst",SmbDrstEnum);
-            iomodel->FetchDataToInput(elements,"md.smb.Rerst",SmbRerstEnum);
-            iomodel->FetchDataToInput(elements,"md.smb.Gdnrst",SmbGdnrstEnum);
-            iomodel->FetchDataToInput(elements,"md.smb.Gsprst",SmbGsprstEnum);
-            iomodel->FetchDataToInput(elements,"md.smb.ECrst",SmbECrstEnum);
-            iomodel->FetchDataToInput(elements,"md.smb.Wrst",SmbWrstEnum);
-            iomodel->FetchDataToInput(elements,"md.smb.Arst",SmbArstEnum);
-            iomodel->FetchDataToInput(elements,"md.smb.Trst",SmbTrstEnum);
-            iomodel->FetchDataToInput(elements,"md.smb.Sizerst",SmbSizerstEnum);
+			InputUpdateFromConstantx(elements,0.,SmbIsInitializedEnum);
+			iomodel->FetchDataToInput(elements,"md.smb.Dzini",SmbDziniEnum);
+			iomodel->FetchDataToInput(elements,"md.smb.Dini",SmbDiniEnum);
+			iomodel->FetchDataToInput(elements,"md.smb.Reini",SmbReiniEnum);
+			iomodel->FetchDataToInput(elements,"md.smb.Gdnini",SmbGdniniEnum);
+			iomodel->FetchDataToInput(elements,"md.smb.Gspini",SmbGspiniEnum);
+			iomodel->FetchDataToInput(elements,"md.smb.ECini",SmbECiniEnum);
+			iomodel->FetchDataToInput(elements,"md.smb.Wini",SmbWiniEnum);
+			iomodel->FetchDataToInput(elements,"md.smb.Aini",SmbAiniEnum);
+			iomodel->FetchDataToInput(elements,"md.smb.Tini",SmbTiniEnum);
+			iomodel->FetchDataToInput(elements,"md.smb.Sizeini",SmbSizeiniEnum);
 			break;
 		case SMBpddEnum:
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 21231)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 21232)
@@ -1475,4 +1475,5 @@
 				name==SmbRefreezeEnum ||
 				name==SmbEvaporationEnum ||
+				name==SmbIsInitializedEnum ||
 				name==BasalforcingsGroundediceMeltingRateEnum ||
 				name==BasalforcingsFloatingiceMeltingRateEnum ||
@@ -2246,7 +2247,6 @@
 
 	/*Intermediary variables: {{{*/
-    bool isinitialized;
-    bool isrestart;
-    IssmDouble zTop,dzTop,zMax,zMin,zY,dzMin;
+	IssmDouble isinitialized;
+	IssmDouble zTop,dzTop,zMax,zMin,zY,dzMin;
 	IssmDouble Tmean; 
 	IssmDouble C; 
@@ -2266,9 +2266,10 @@
 	IssmDouble lhf, shf, dayEC;
 	IssmDouble initMass;
-    IssmDouble sumR, sumM, sumEC, sumP, sumW,sumMassAdd;
-    IssmDouble sumdz_add;
-    IssmDouble sumMass,dMass;
+	IssmDouble sumR, sumM, sumEC, sumP, sumW,sumMassAdd;
+	IssmDouble sumdz_add;
+	IssmDouble sumMass,dMass;
 	bool isgraingrowth,isalbedo,isshortwave,isthermal,isaccumulation,ismelt,isdensification,isturbulentflux;
 	IssmDouble init_scaling;
+
 	/*}}}*/
 	/*Output variables:{{{ */
@@ -2287,15 +2288,15 @@
 	IssmDouble  R; 
 	IssmDouble  mAdd;
-    IssmDouble  dz_add;
+	IssmDouble  dz_add;
     
-    IssmDouble* dzrst=NULL;
-    IssmDouble* drst = NULL;
-    IssmDouble* rerst = NULL;
-    IssmDouble* gdnrst = NULL;
-    IssmDouble* gsprst = NULL;
-    IssmDouble* Wrst = NULL;
-    IssmDouble* arst = NULL;
-    IssmDouble* Trst = NULL;
-
+	IssmDouble* dzini=NULL;
+	IssmDouble* dini = NULL;
+	IssmDouble* reini = NULL;
+	IssmDouble* gdnini = NULL;
+	IssmDouble* gspini = NULL;
+	IssmDouble* Wini = NULL;
+	IssmDouble* aini = NULL;
+	IssmDouble* Tini = NULL;
+    
 	int         m;
 	int         count=0;
@@ -2313,5 +2314,5 @@
 	parameters->FindParam(&time,TimeEnum);                        /*transient core time at which we run the smb core*/
 	parameters->FindParam(&dt,TimesteppingTimeStepEnum);          /*transient core time step*/
-    parameters->FindParam(&yts,ConstantsYtsEnum);
+	parameters->FindParam(&yts,ConstantsYtsEnum);
 	parameters->FindParam(&smb_dt,SmbDtEnum);                     /*time period for the smb solution,  usually smaller than the glaciological dt*/
 	parameters->FindParam(&aIdx,SmbAIdxEnum);
@@ -2331,5 +2332,5 @@
 	parameters->FindParam(&isturbulentflux,SmbIsturbulentfluxEnum);
 	parameters->FindParam(&init_scaling,SmbInitDensityScalingEnum);
-
+    
 	/*}}}*/
 	/*Retrieve inputs: {{{*/
@@ -2351,7 +2352,5 @@
 	Input* eAir_input=this->GetInput(SmbEAirEnum); _assert_(eAir_input);
 	Input* pAir_input=this->GetInput(SmbPAirEnum); _assert_(pAir_input);
-    Input* isinitialized_input=this->GetInput(SmbIsInitializedEnum);  _assert_(isinitialized_input);
-    Input* isrestart_input=this->GetInput(SmbIsrestartEnum);  _assert_(isrestart_input);
-    
+	Input* isinitialized_input=this->GetInput(SmbIsInitializedEnum); _assert_(isinitialized_input);
 	/*Retrieve input values:*/
 	Gauss* gauss=this->NewGauss(1); gauss->GaussPoint(0);
@@ -2367,100 +2366,100 @@
 	Tz_input->GetInputValue(&Tz,gauss);
 	Vz_input->GetInputValue(&Vz,gauss);
-    isinitialized_input->GetInputValue(&isinitialized);
-    isrestart_input->GetInputValue(&isrestart);
+	isinitialized_input->GetInputValue(&isinitialized);
 	/*}}}*/
 
-    /*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){
+	/*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(VerboseSmb() && this->Sid()==0)_printf0_("Snow properties initialized w DEFAULT values\n");
-        GembgridInitialize(&dz, &m, zTop, dzTop, zMax, zY);
         //if(this->Sid()==1) for(int i=0;i<m;i++)_printf_("z[" << i << "]=" <<
         //dz[i] << "\n");
-        /*initialize profile variables:*/
-        d = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)d[i]=rho_ice*init_scaling; //ice density scaled by a factor
-        re = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)re[i]=2.5;         //set grain size to old snow [mm]
-        gdn = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)gdn[i]=0;         //set grain dentricity to old snow
-        gsp = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)gsp[i]=0;         //set grain sphericity to old snow
-        EC = 0;                                                                //surface evaporation (-) condensation (+) [kg m-2]
-        W = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)W[i]=0;             //set water content to zero [kg m-2]
-        a = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)a[i]=aSnow;         //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]
         
+        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];
-        
-        /*Flag the initialization:*/
-        this->AddInput(new BoolInput(SmbIsInitializedEnum,false));
-    }
-    else if(isrestart){ //Retrieve the snow properties from previous run
-        if(VerboseSmb() && this->Sid()==0)_printf0_("Snow properties initialized w RESTART values\n");
-        
-        /*Recover inputs: */
-        DoubleArrayInput* dzrst_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDzrstEnum)); _assert_(dzrst_input);
-        DoubleArrayInput* drst_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDrstEnum));_assert_(drst_input);
-        DoubleArrayInput* rerst_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbRerstEnum));_assert_(rerst_input);
-        DoubleArrayInput* gdnrst_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGdnrstEnum));_assert_(gdnrst_input);
-        DoubleArrayInput* gsprst_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGsprstEnum));_assert_(gsprst_input);
-        DoubleInput* ECrst_input= dynamic_cast<DoubleInput*>(this->GetInput(SmbECrstEnum));_assert_(ECrst_input);
-        DoubleArrayInput* Wrst_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbWrstEnum));_assert_(Wrst_input);
-        DoubleArrayInput* arst_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbArstEnum));_assert_(arst_input);
-        DoubleArrayInput* Trst_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbTrstEnum));_assert_(Trst_input);
-        
-        /*Recover arrays: */
-        dzrst_input->GetValues(&dzrst,&m);
-        drst_input->GetValues(&drst,&m);
-        rerst_input->GetValues(&rerst,&m);
-        gdnrst_input->GetValues(&gdnrst,&m);
-        gsprst_input->GetValues(&gsprst,&m);
-        ECrst_input->GetInputValue(&EC); 
-        Wrst_input->GetValues(&Wrst,&m);
-        arst_input->GetValues(&arst,&m);
-        Trst_input->GetValues(&Trst,&m);
-        
-        /*Retrive the correct value of m (without the zeroes at the end)*/
-        Input* Sizerst_input=this->GetInput(SmbSizerstEnum); _assert_(Sizerst_input);
-        Sizerst_input->GetInputValue(&m);
-        
-        dz = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)dz[i]=dzrst[i];
-        d = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)d[i]=drst[i];
-        re = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)re[i]=rerst[i];
-        gdn = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)gdn[i]=gdnrst[i];
-        gsp = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)gsp[i]=gsprst[i];
-        W = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)W[i]=Wrst[i];
-        a = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)a[i]=arst[i];
-        T = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)T[i]=Trst[i];
-        
-        //fixed lower temperatuer bounday condition - T is fixed
-        T_bottom=T[m-1];
-        
-        /*Flag the initialization:*/
-         this->AddInput(new BoolInput(SmbIsrestartEnum,false));
-    }
-    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 temperatuer bounday condition - T is fixed
-		T_bottom=T[m-1];
 
     } /*}}}*/
@@ -2469,7 +2468,7 @@
 	initMass=0; for(int i=0;i<m;i++) initMass += dz[i]*d[i] + W[i];
     
-    // initialize cumulative variables
-    sumR = 0; sumM = 0; sumEC = 0; sumP = 0; sumMassAdd = 0;
-    sumdz_add=0;
+        // initialize cumulative variables
+	sumR = 0; sumM = 0; sumEC = 0; sumP = 0; sumMassAdd = 0;
+	sumdz_add=0;
 
 	//before starting loop, realize that the transient core runs this smb_core at time = time +deltaT. 
@@ -2502,35 +2501,35 @@
 					
 		/*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,m,this->Sid());
 		
 		/*Calculate net shortwave [W m-2]*/
-        netSW = cellsum(swf,m);
+		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,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];
+		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());
-
-		/*Calculate water production, M [kg m-2] resulting from snow/ice temperature exceeding 273.15 deg K 
+		if(isaccumulation)accumulation(&T, &dz, &d, &W, &a, &re, &gdn, &gsp, &m, Ta, P, dzMin, aSnow,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,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());
+		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*/
-        ulw = 5.67E-8 * pow(T[0],4.0);
+		ulw = 5.67E-8 * pow(T[0],4.0);
 
 		/*Calculate net longwave [W m-2]*/
-        netLW = dlw - ulw;
+		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,this->Sid());
 		
 		/*Verbose some resuls in debug mode: {{{*/
@@ -2550,25 +2549,25 @@
 		
 		/*Sum component mass changes [kg m-2]*/
-        sumMassAdd = mAdd + sumMassAdd;
-        sumM = M + sumM;
-        sumR = R + sumR;
-        sumW = cellsum(W,m);
-        sumP = P +  sumP;
-        sumEC = sumEC + EC;  // evap (-)/cond(+)
-        sumdz_add=dz_add+sumdz_add; //*CL*
+		sumMassAdd = mAdd + sumMassAdd;
+		sumM = M + sumM;
+		sumR = R + sumR;
+		sumW = cellsum(W,m);
+		sumP = P +  sumP;
+		sumEC = sumEC + EC;  // evap (-)/cond(+)
+		sumdz_add=dz_add+sumdz_add;
 
 		/*Calculate total system mass:*/
-        sumMass=0; for(int i=0;i<m;i++) sumMass += dz[i]*d[i];
-
-        #ifndef _HAVE_ADOLC_ //we want to avoid the round operation at all cost. Not differentiable.
-        dMass = sumMass + sumR + sumW - sumP - sumEC - initMass - sumMassAdd;
+		sumMass=0; for(int i=0;i<m;i++) sumMass += dz[i]*d[i];
+
+		#ifndef _HAVE_ADOLC_ //we want to avoid the round operation at all cost. Not differentiable.
+		dMass = sumMass + sumR + sumW - sumP - sumEC - initMass - sumMassAdd;
 		dMass = round(dMass * 100.0)/100.0;
 
 		/*Check mass conservation:*/
-        if (dMass != 0.0) _printf_("total system mass not conserved in MB function");
+		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");
+		if (T[m-1]!=T_bottom) _printf_("T(end)~=T_bottom" << "\n");
 		
 		/*Free ressources: */
@@ -2579,5 +2578,5 @@
 	} //for (t=time;t<time+dt;t=t+smb_dt)
 
-    /*Save generated inputs: */
+	/*Save generated inputs: */
 	this->AddInput(new DoubleArrayInput(SmbDzEnum,dz,m));
 	this->AddInput(new DoubleArrayInput(SmbDEnum,d,m));
@@ -2589,11 +2588,11 @@
 	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(SmbDz_addEnum,sumdz_add/yts));
-    this->AddInput(new DoubleInput(SmbM_addEnum,sumMassAdd/yts));
-
-    /*Free allocations:{{{*/
+	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(SmbDz_addEnum,sumdz_add/yts));
+	this->AddInput(new DoubleInput(SmbM_addEnum,sumMassAdd/yts));
+
+	/*Free allocations:{{{*/
 	xDelete<IssmDouble>(dz);
 	xDelete<IssmDouble>(d);
Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 21231)
+++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 21232)
@@ -12,6 +12,6 @@
 
 	for(int i=0;i<femmodel->elements->Size();i++){
-		Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
-		element->SmbGemb();
+        Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+        element->SmbGemb();
 	}
 
@@ -67,5 +67,5 @@
 		gpBottom++;
 	}
-   //initialize bottom vectors
+	//initialize bottom vectors
 	dzB = xNewZeroInit<IssmDouble>(gpBottom);
 	gp0 = dzTop;
@@ -98,5 +98,5 @@
 
 	// calculates grain growth according to Fig. 9 of Marbouty, 1980
-	// ------NO GRAIN GROWTH FOR d > 400 kg m-3 because (H is set to zero)------
+	// ------NO GRAIN GROWTH FOR d > 400 kg m-3 because H is set to zero------
 	// ---------------this is a major limitation of the model-------------------
 
@@ -354,8 +354,8 @@
 
 	if(aIdx==1){ 
-		//function of effective grain radius
+        //function of effective grain radius
         
         //convert effective radius to specific surface area [cm2 g-1]
-		IssmDouble S = 3.0 / (.091 * re[0]);
+        IssmDouble S = 3.0 / (.091 * re[0]);
         
         //determine broadband albedo
@@ -364,8 +364,8 @@
 	else if(aIdx==2){
 		
-		// Spectral fractions  (Lefebre et al., 2003)
+        // Spectral fractions  (Lefebre et al., 2003)
         // [0.3-0.8um 0.8-1.5um 1.5-2.8um]
         
-		IssmDouble sF[3] = {0.606, 0.301, 0.093};
+        IssmDouble sF[3] = {0.606, 0.301, 0.093};
         
         // convert effective radius to grain size in meters
@@ -386,5 +386,5 @@
 	else if(aIdx==3){
 		
-		// a as a function of density
+        // a as a function of density
         
         // calculate albedo
@@ -393,5 +393,5 @@
 	else if(aIdx==4){
 		
-		// exponential time decay & wetness
+        // exponential time decay & wetness
         
         // change in albedo with time:
@@ -402,8 +402,8 @@
         
         // initialize variables
-		IssmDouble* t0=xNew<IssmDouble>(m);
-		IssmDouble* T=xNew<IssmDouble>(m);
-		IssmDouble* t0warm=xNew<IssmDouble>(m);
-		IssmDouble* d_a=xNew<IssmDouble>(m);
+        IssmDouble* t0=xNew<IssmDouble>(m);
+        IssmDouble* T=xNew<IssmDouble>(m);
+        IssmDouble* t0warm=xNew<IssmDouble>(m);
+        IssmDouble* d_a=xNew<IssmDouble>(m);
         
         // specify constants
@@ -418,7 +418,7 @@
         
         // 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++) t0warm[i]= fabs(T[i]) * K + t0dry; //// 'warm' snow timescale
+        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++) 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
@@ -440,9 +440,9 @@
         // a_surf = a_wet - (a_wet - a_surf) * exp(-W_surf/W0);
 
-		/*Free ressources:*/
-		xDelete<IssmDouble>(t0);
-		xDelete<IssmDouble>(T);
-		xDelete<IssmDouble>(t0warm);
-		xDelete<IssmDouble>(d_a);
+        /*Free ressources:*/
+        xDelete<IssmDouble>(t0);
+        xDelete<IssmDouble>(T);
+        xDelete<IssmDouble>(t0warm);
+        xDelete<IssmDouble>(d_a);
 
 	}
@@ -1219,22 +1219,22 @@
 	IssmDouble Wi;
     
-    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;
-    bool        top=false;
+	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;
+	bool        top=false;
     
-    IssmDouble* Zcum=NULL;
-    IssmDouble* dzMin2=NULL;
-    IssmDouble zY2=1.025;
-    bool lastCellFlag;
-    int X1=0;
-    int X2=0;
+	IssmDouble* Zcum=NULL;
+	IssmDouble* dzMin2=NULL;
+	IssmDouble zY2=1.025;
+	bool lastCellFlag;
+	int X1=0;
+	int X2=0;
     
 	int        D_size;
@@ -1242,5 +1242,5 @@
 	/*outputs:*/
 	IssmDouble  mAdd;
-    IssmDouble dz_add;
+	IssmDouble dz_add;
 	IssmDouble  Rsum;
 	IssmDouble* T=*pT;
@@ -1284,5 +1284,5 @@
 	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]
+	dz_add=0;      // thickness of the layer added/removed to/from base of model [m]
 
 	// calculate temperature excess above 0 deg C
@@ -1322,37 +1322,37 @@
 	}
 
-    //// MELT, PERCOLATION AND REFREEZE
-    
-    // run melt algorithm if there is melt water or excess pore water
-    if ((cellsum(exsT,n) > 0) || (cellsum(exsW,n) > 0)){
-        // _printf_(""MELT OCCURS");
-        // check to see if thermal energy exceeds energy to melt entire cell
-        // if so redistribute temperature to lower cells (temperature surplus)
-        // (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 ){
-            // _printf_("T Surplus");
-            // calculate surplus energy
-            surpE=xNew<IssmDouble>(n); for(int i=0;i<n;i++)surpE[i] = surpT[i] * CI * m[i];//*CL*
+	//// MELT, PERCOLATION AND REFREEZE
+
+	// run melt algorithm if there is melt water or excess pore water
+	if ((cellsum(exsT,n) > 0) || (cellsum(exsW,n) > 0)){
+		// _printf_(""MELT OCCURS");
+		// check to see if thermal energy exceeds energy to melt entire cell
+		// if so redistribute temperature to lower cells (temperature surplus)
+		// (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 ){
+			// _printf_("T Surplus");
+			// calculate surplus energy
+			surpE=xNew<IssmDouble>(n); for(int i=0;i<n;i++)surpE[i] = surpT[i] * CI * m[i];
             
-            int i = 0;
-            while (cellsum(surpE,n) > 0){
-                // use surplus energy to increase the temperature of lower cell
-                T[i+1] = surpE[i]/m[i+1]/CI + T[i+1];//*CL*
+			int i = 0;
+			while (cellsum(surpE,n) > 0){
+				// 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];
-                T[i+1] = fmin(CtoK, T[i+1]);
+				exsT[i+1] = fmax(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);
-                surpE[i+1] = surpT[i+1] * CI * m[i+1];//*CL*
+				surpT[i+1] = fmax(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;   
-                i = i + 1;
-            }
-        }
+				// adjust current cell properties (again 159.1342 is the max T)
+				exsT[i] = 159.1342;
+				surpE[i] = 0;
+				i = i + 1;
+			}
+		}
 
 		// convert temperature excess to melt [kg]
@@ -1364,6 +1364,6 @@
 
 		// initialize refreeze, runoff, flxDn and dW vectors [kg]
-		IssmDouble* F = xNewZeroInit<IssmDouble>(n); 
-		IssmDouble* R=xNewZeroInit<IssmDouble>(n); 
+ 		IssmDouble* F = xNewZeroInit<IssmDouble>(n);
+		IssmDouble* R=xNewZeroInit<IssmDouble>(n);
 
 		for(int i=0;i<n;i++)dW[i] = 0;
@@ -1390,5 +1390,5 @@
 
 			// 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){  // dPHC = pore hole close off [kg m-3]
 				// _printf_("ICE LAYER");
 				// no water freezes in this cell
@@ -1401,5 +1401,5 @@
 				R[i] = fmax(0, inM - dW[i]);             // runoff
 			}
-			// check if no energy to refreeze meltwater     
+			// check if no energy to refreeze meltwater
 			else if (maxF[i] == 0){
 				// _printf_("REFREEZE == 0");
@@ -1432,5 +1432,5 @@
 				IssmDouble F2 = 0;                                 
 
-				if (dW[i] < 0){                            // excess pore water
+				if (dW[i] < 0){                         // excess pore water
 					dMax = (dIce - d[i])*dz_0;          // maximum refreeze                                             
 					IssmDouble maxF2 = fmin(dMax, maxF[i]-F1);      // maximum refreeze
@@ -1464,5 +1464,5 @@
 
 		// 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++;}
 		
@@ -1484,5 +1484,5 @@
 		for(int i=0;i<n;i++)dz[i] = m[i] / d[i];
 
-		//calculate Rsum: 
+		//calculate Rsum:
 		Rsum=cellsum(R,n);
 
@@ -1492,27 +1492,27 @@
 	}
     
-    //Merging of cells as they are burried under snow.
-    Zcum=xNew<IssmDouble>(n);
-    dzMin2=xNew<IssmDouble>(n);
+	//Merging of cells as they are burried under snow.
+	Zcum=xNew<IssmDouble>(n);
+	dzMin2=xNew<IssmDouble>(n);
     
-    Zcum[0]=dz[0]; // Compute a cumulative depth vector
+	Zcum[0]=dz[0]; // Compute a cumulative depth vector
     
-    for (int i=1;i<n;i++){
-        Zcum[i]=Zcum[i-1]+dz[i];
-    }
+	for (int i=1;i<n;i++){
+		Zcum[i]=Zcum[i-1]+dz[i];
+	}
     
-    for (int i=0;i<n;i++){
-        if (Zcum[i]<=zTop){
-            dzMin2[i]=dzMin;
-        }
-        else{
-            dzMin2[i]=zY2*dzMin2[i-1];
-        }
-    }
+	for (int i=0;i<n;i++){
+		if (Zcum[i]<=zTop){
+			dzMin2[i]=dzMin;
+		}
+		else{
+			dzMin2[i]=zY2*dzMin2[i-1];
+		}
+	}
 
 	// check if depth is too small
 	X = 0;
 	for(int i=n-1;i>=0;i--){
-         if(dz[i]<dzMin2[i]){
+		if(dz[i]<dzMin2[i]){
 			X=i;
 			break;
@@ -1520,63 +1520,63 @@
 	}
     
-    //Last cell has to be treated separately if has to be merged (no underlying cell so need to merge with overlying cell)
-    if(X==n-1){
-        lastCellFlag = true;
-        X = X-1;
-    }
-    else{
-        lastCellFlag = false;
-    }
-
-    for (int i = 0; i<=X;i++){
-        if (dz [i] < dzMin2[i]){                               // merge top cells
-            //                                                                          _printf_("dz > dzMin * 2')
-            // adjust variables as a linearly weighted function of mass
-            IssmDouble m_new = m[i] + m[i+1];
-            T[i+1] = (T[i]*m[i] + T[i+1]*m[i+1]) / m_new;
-            a[i+1] = (a[i]*m[i] + a[i+1]*m[i+1]) / m_new;
-            re[i+1] = (re[i]*m[i] + re[i+1]*m[i+1]) / m_new;
-            gdn[i+1] = (gdn[i]*m[i] + gdn[i+1]*m[i+1]) / m_new;
-            gsp[i+1] = (gsp[i]*m[i] + gsp[i+1]*m[i+1]) / m_new;
+	//Last cell has to be treated separately if has to be merged (no underlying cell so need to merge with overlying cell)
+	if(X==n-1){
+		lastCellFlag = true;
+		X = X-1;
+	}
+	else{
+		lastCellFlag = false;
+	}
+
+	for (int i = 0; i<=X;i++){
+		if (dz [i] < dzMin2[i]){                              // merge top cells
+			// _printf_("dz > dzMin * 2')
+			// adjust variables as a linearly weighted function of mass
+			IssmDouble m_new = m[i] + m[i+1];
+			T[i+1] = (T[i]*m[i] + T[i+1]*m[i+1]) / m_new;
+			a[i+1] = (a[i]*m[i] + a[i+1]*m[i+1]) / m_new;
+			re[i+1] = (re[i]*m[i] + re[i+1]*m[i+1]) / m_new;
+			gdn[i+1] = (gdn[i]*m[i] + gdn[i+1]*m[i+1]) / m_new;
+			gsp[i+1] = (gsp[i]*m[i] + gsp[i+1]*m[i+1]) / m_new;
             
-            // merge with underlying grid cell and delete old cell
-            dz [i+1] = dz[i] + dz[i+1];                 // combine cell depths
-            d[i+1] = m_new / dz[i+1];                   // combine top densities
-            W[i+1] = W[i+1] + W[i];                     // combine liquid water
-            m[i+1] = m_new;                             // combine top masses
+			// merge with underlying grid cell and delete old cell
+			dz [i+1] = dz[i] + dz[i+1];                 // combine cell depths
+			d[i+1] = m_new / dz[i+1];                   // combine top densities
+			W[i+1] = W[i+1] + W[i];                     // combine liquid water
+			m[i+1] = m_new;                             // combine top masses
             
-            // set cell to 99999 for deletion
-            m[i] = 99999;
-        }
-    }
-
-    //If last cell has to be merged
-    if(lastCellFlag){
-        //find closest cell to merge with
-        for(int i=n-2;i>=0;i--){
-            if(m[i]!=99999){
-                X2=i;
-                X1=n-1;
-                break;
-            }
-        }
-        
-        // adjust variables as a linearly weighted function of mass
-        IssmDouble m_new = m[X2] + m[X1];
-        T[X1] = (T[X2]*m[X2] + T[X1]*m[X1]) / m_new;
-        a[X1] = (a[X2]*m[X2] + a[X1]*m[X1]) / m_new;
-        re[X1] = (re[X2]*m[X2] + re[X1]*m[X1]) / m_new;
-        gdn[X1] = (gdn[X2]*m[X2] + gdn[X1]*m[X1]) / m_new;
-        gsp[X1] = (gsp[X2]*m[X2] + gsp[X1]*m[X1]) / m_new;
-        
-        // merge with underlying grid cell and delete old cell
-        dz [X1] = dz[X2] + dz[X1];                 // combine cell depths
-        d[X1] = m_new / dz[X1];                   // combine top densities
-        W[X1] = W[X1] + W[X2];                     // combine liquid water
-        m[X1] = m_new;                             // combine top masses
-        
-        // set cell to 99999 for deletion
-        m[X2] = 99999;
-    }
+			// set cell to 99999 for deletion
+			m[i] = 99999;
+		}
+	}
+
+	//If last cell has to be merged
+	if(lastCellFlag){
+         //find closest cell to merge with
+		for(int i=n-2;i>=0;i--){
+			if(m[i]!=99999){
+				X2=i;
+				X1=n-1;
+				break;
+			}
+		}
+        
+		// adjust variables as a linearly weighted function of mass
+		IssmDouble m_new = m[X2] + m[X1];
+		T[X1] = (T[X2]*m[X2] + T[X1]*m[X1]) / m_new;
+		a[X1] = (a[X2]*m[X2] + a[X1]*m[X1]) / m_new;
+		re[X1] = (re[X2]*m[X2] + re[X1]*m[X1]) / m_new;
+		gdn[X1] = (gdn[X2]*m[X2] + gdn[X1]*m[X1]) / m_new;
+		gsp[X1] = (gsp[X2]*m[X2] + gsp[X1]*m[X1]) / m_new;
+        
+		// merge with underlying grid cell and delete old cell
+  		dz [X1] = dz[X2] + dz[X1];                 // combine cell depths
+		d[X1] = m_new / dz[X1];                   // combine top densities
+ 		W[X1] = W[X1] + W[X2];                     // combine liquid water
+		m[X1] = m_new;                             // combine top masses
+        
+		// set cell to 99999 for deletion
+		m[X2] = 99999;
+	}
 
 	// delete combined cells
@@ -1597,5 +1597,5 @@
 	n=D_size;
 	xDelete<int>(D);
-
+    
 	// check if any of the top 10 cell depths are too large
 	X=0;
@@ -1611,84 +1611,84 @@
 		if (dz [i] > dzMin *2){
 
-			// _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);
-			n++;
-			X=X+1;
+				// _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);
+				n++;
+				X=X+1;
 		}
 		else i++;
 	}
 
-    //// CORRECT FOR TOTAL MODEL DEPTH
-    // WORKS FINE BUT HAS BEEN DISABLED FOR CONVIENCE OF MODEL OUTPUT
-    // INTERPRETATION
-    // // calculate total model depth
-    Ztot = cellsum(dz,n);
+	//// CORRECT FOR TOTAL MODEL DEPTH
+	// WORKS FINE BUT HAS BEEN DISABLED FOR CONVIENCE OF MODEL OUTPUT
+	// INTERPRETATION
+	// // calculate total model depth
+	Ztot = cellsum(dz,n);
     
-    if (Ztot < zMin){
-//        printf("Total depth < zMin %f \n", Ztot);
-        // mass and energy to be added
-        mAdd = m[n-1]+W[n-1];
-        addE = T[n-1]*m[n-1]*CI;
-        
-        // add a grid cell of the same size and temperature to the bottom
-        dz_bot=dz[n-1];
-        T_bot=T[n-1];
-        W_bot=W[n-1];
-        m_bot=m[n-1];
-        d_bot=d[n-1];
-        a_bot=a[n-1];
-        re_bot=re[n-1];
-        gdn_bot=gdn[n-1];
-        gsp_bot=gsp[n-1];
-        
-        dz_add=dz_bot;
-        
-        newcell(&dz,dz_bot,top,n);
-        newcell(&T,T_bot,top,n);
-        newcell(&W,W_bot,top,n);
-        newcell(&m,m_bot,top,n);
-        newcell(&d,d_bot,top,n);
-        newcell(&a,a_bot,top,n);
-        newcell(&re,re_bot,top,n);
-        newcell(&gdn,gdn_bot,top,n);
-        newcell(&gsp,gsp_bot,top,n);
-        n=n+1;
-    }
-    else if (Ztot > zMax){
-//        printf("Total depth > zMax %f \n", Ztot);
-        // mass and energy loss
-        mAdd = -(m[n-1]+W[n-1]);
-        addE = -(T[n-1]*m[n-1]*CI);
-        dz_add=-(dz[n-1]);
-        
-        // add a grid cell of the same size and temperature to the bottom
-        D_size=n-1;
-        D=xNew<int>(D_size);
-        
-        for(int i=0;i<n-1;i++) D[i]=i;
-        celldelete(&dz,n,D,D_size);
-        celldelete(&T,n,D,D_size);
-        celldelete(&W,n,D,D_size);
-        celldelete(&m,n,D,D_size);
-        celldelete(&d,n,D,D_size);
-        celldelete(&a,n,D,D_size);
-        celldelete(&re,n,D,D_size);
-        celldelete(&gdn,n,D,D_size);
-        celldelete(&gsp,n,D,D_size);
-        n=D_size;
-        xDelete<int>(D);
-        
-    }
+	if (Ztot < zMin){
+		// printf("Total depth < zMin %f \n", Ztot);
+		// mass and energy to be added
+		mAdd = m[n-1]+W[n-1];
+		addE = T[n-1]*m[n-1]*CI;
+        
+		// add a grid cell of the same size and temperature to the bottom
+		dz_bot=dz[n-1];
+		T_bot=T[n-1];
+		W_bot=W[n-1];
+		m_bot=m[n-1];
+		d_bot=d[n-1];
+		a_bot=a[n-1];
+		re_bot=re[n-1];
+		gdn_bot=gdn[n-1];
+		gsp_bot=gsp[n-1];
+        
+		dz_add=dz_bot;
+        
+		newcell(&dz,dz_bot,top,n);
+		newcell(&T,T_bot,top,n);
+		newcell(&W,W_bot,top,n);
+		newcell(&m,m_bot,top,n);
+		newcell(&d,d_bot,top,n);
+		newcell(&a,a_bot,top,n);
+		newcell(&re,re_bot,top,n);
+		newcell(&gdn,gdn_bot,top,n);
+		newcell(&gsp,gsp_bot,top,n);
+		n=n+1;
+	}
+	else if (Ztot > zMax){
+		// printf("Total depth > zMax %f \n", Ztot);
+		// mass and energy loss
+		mAdd = -(m[n-1]+W[n-1]);
+		addE = -(T[n-1]*m[n-1]*CI);
+		dz_add=-(dz[n-1]);
+        
+		// add a grid cell of the same size and temperature to the bottom
+		D_size=n-1;
+		D=xNew<int>(D_size);
+        
+		for(int i=0;i<n-1;i++) D[i]=i;
+		celldelete(&dz,n,D,D_size);
+		celldelete(&T,n,D,D_size);
+		celldelete(&W,n,D,D_size);
+		celldelete(&m,n,D,D_size);
+		celldelete(&d,n,D,D_size);
+		celldelete(&a,n,D,D_size);
+		celldelete(&re,n,D,D_size);
+		celldelete(&gdn,n,D,D_size);
+		celldelete(&gsp,n,D,D_size);
+		n=D_size;
+		xDelete<int>(D);
+        
+	}
 
 	//// CHECK FOR MASS AND ENERGY CONSERVATION
@@ -1709,5 +1709,5 @@
 	dm = round(mSum0 - mSum1 + mAdd);
 	dE = round(sumE0 - sumE1 - sumER +  addE);
-	if (dm !=0  || dE !=0) _error_("mass or energy are not conserved in melt equations\n" 
+	if (dm !=0  || dE !=0) _error_("mass or energy are not conserved in melt equations\n"
 			<< "dm: " << dm << " dE: " << dE << "\n");
 	#endif
@@ -1726,6 +1726,6 @@
 	if(D)xDelete<int>(D);
 	if(M)xDelete<IssmDouble>(M);
-    xDelete<IssmDouble>(Zcum);
-    xDelete<IssmDouble>(dzMin2);
+ 	xDelete<IssmDouble>(Zcum);
+	xDelete<IssmDouble>(dzMin2);
     
 	/*Assign output pointers:*/
@@ -1733,5 +1733,5 @@
 	*pR=Rsum;
 	*pmAdd=mAdd;
-    *pdz_add=dz_add;
+	*pdz_add=dz_add;
 
 	*pT=T;
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 21231)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 21232)
@@ -166,5 +166,5 @@
 	HydrologysommersEnum,
 	HydrologyHeadEnum,
-   HydrologyHeadOldEnum,
+	HydrologyHeadOldEnum,
 	HydrologyGapHeightEnum,
 	HydrologyBumpSpacingEnum,
@@ -361,15 +361,14 @@
 	SmbRequestedOutputsEnum,
 	SmbIsInitializedEnum,
-    SmbIsrestartEnum,
-    SmbDzrstEnum,
-    SmbDrstEnum,
-    SmbRerstEnum,
-    SmbGdnrstEnum,
-    SmbGsprstEnum,
-    SmbECrstEnum,
-    SmbWrstEnum,
-    SmbArstEnum,
-    SmbTrstEnum,
-    SmbSizerstEnum,
+	SmbDziniEnum,
+	SmbDiniEnum,
+	SmbReiniEnum,
+	SmbGdniniEnum,
+	SmbGspiniEnum,
+	SmbECiniEnum,
+	SmbWiniEnum,
+	SmbAiniEnum,
+	SmbTiniEnum,
+	SmbSizeiniEnum,
 	/*SMBforcing*/
 	SMBforcingEnum,
@@ -423,6 +422,6 @@
 	SmbIsdensificationEnum,
 	SmbIsturbulentfluxEnum,
-    SmbDz_addEnum,
-    SmbM_addEnum,
+	SmbDz_addEnum,
+	SmbM_addEnum,
 	/*SMBpdd*/
 	SMBpddEnum,	
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 21231)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 21232)
@@ -365,15 +365,14 @@
 		case SmbRequestedOutputsEnum : return "SmbRequestedOutputs";
 		case SmbIsInitializedEnum : return "SmbIsInitialized";
-		case SmbIsrestartEnum : return "SmbIsrestart";
-		case SmbDzrstEnum : return "SmbDzrst";
-		case SmbDrstEnum : return "SmbDrst";
-		case SmbRerstEnum : return "SmbRerst";
-		case SmbGdnrstEnum : return "SmbGdnrst";
-		case SmbGsprstEnum : return "SmbGsprst";
-		case SmbECrstEnum : return "SmbECrst";
-		case SmbWrstEnum : return "SmbWrst";
-		case SmbArstEnum : return "SmbArst";
-		case SmbTrstEnum : return "SmbTrst";
-		case SmbSizerstEnum : return "SmbSizerst";
+		case SmbDziniEnum : return "SmbDzini";
+		case SmbDiniEnum : return "SmbDini";
+		case SmbReiniEnum : return "SmbReini";
+		case SmbGdniniEnum : return "SmbGdnini";
+		case SmbGspiniEnum : return "SmbGspini";
+		case SmbECiniEnum : return "SmbECini";
+		case SmbWiniEnum : return "SmbWini";
+		case SmbAiniEnum : return "SmbAini";
+		case SmbTiniEnum : return "SmbTini";
+		case SmbSizeiniEnum : return "SmbSizeini";
 		case SMBforcingEnum : return "SMBforcing";
 		case SmbMassBalanceEnum : return "SmbMassBalance";
Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 21231)
+++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 21232)
@@ -371,21 +371,20 @@
 	      else if (strcmp(name,"SmbRequestedOutputs")==0) return SmbRequestedOutputsEnum;
 	      else if (strcmp(name,"SmbIsInitialized")==0) return SmbIsInitializedEnum;
-	      else if (strcmp(name,"SmbIsrestart")==0) return SmbIsrestartEnum;
-	      else if (strcmp(name,"SmbDzrst")==0) return SmbDzrstEnum;
-	      else if (strcmp(name,"SmbDrst")==0) return SmbDrstEnum;
-	      else if (strcmp(name,"SmbRerst")==0) return SmbRerstEnum;
-	      else if (strcmp(name,"SmbGdnrst")==0) return SmbGdnrstEnum;
-	      else if (strcmp(name,"SmbGsprst")==0) return SmbGsprstEnum;
-	      else if (strcmp(name,"SmbECrst")==0) return SmbECrstEnum;
-	      else if (strcmp(name,"SmbWrst")==0) return SmbWrstEnum;
-	      else if (strcmp(name,"SmbArst")==0) return SmbArstEnum;
-	      else if (strcmp(name,"SmbTrst")==0) return SmbTrstEnum;
-	      else if (strcmp(name,"SmbSizerst")==0) return SmbSizerstEnum;
+	      else if (strcmp(name,"SmbDzini")==0) return SmbDziniEnum;
+	      else if (strcmp(name,"SmbDini")==0) return SmbDiniEnum;
+	      else if (strcmp(name,"SmbReini")==0) return SmbReiniEnum;
+	      else if (strcmp(name,"SmbGdnini")==0) return SmbGdniniEnum;
+	      else if (strcmp(name,"SmbGspini")==0) return SmbGspiniEnum;
+	      else if (strcmp(name,"SmbECini")==0) return SmbECiniEnum;
+	      else if (strcmp(name,"SmbWini")==0) return SmbWiniEnum;
+	      else if (strcmp(name,"SmbAini")==0) return SmbAiniEnum;
+	      else if (strcmp(name,"SmbTini")==0) return SmbTiniEnum;
+	      else if (strcmp(name,"SmbSizeini")==0) return SmbSizeiniEnum;
 	      else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum;
+	      else if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum;
          else stage=4;
    }
    if(stage==4){
-	      if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum;
-	      else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum;
+	      if (strcmp(name,"SMBgemb")==0) return SMBgembEnum;
 	      else if (strcmp(name,"SmbInitDensityScaling")==0) return SmbInitDensityScalingEnum;
 	      else if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
@@ -506,9 +505,9 @@
 	      else if (strcmp(name,"VxAverage")==0) return VxAverageEnum;
 	      else if (strcmp(name,"Vx")==0) return VxEnum;
+	      else if (strcmp(name,"VxPicard")==0) return VxPicardEnum;
          else stage=5;
    }
    if(stage==5){
-	      if (strcmp(name,"VxPicard")==0) return VxPicardEnum;
-	      else if (strcmp(name,"VyAverage")==0) return VyAverageEnum;
+	      if (strcmp(name,"VyAverage")==0) return VyAverageEnum;
 	      else if (strcmp(name,"Vy")==0) return VyEnum;
 	      else if (strcmp(name,"VyPicard")==0) return VyPicardEnum;
@@ -629,9 +628,9 @@
 	      else if (strcmp(name,"Outputdefinition46")==0) return Outputdefinition46Enum;
 	      else if (strcmp(name,"Outputdefinition47")==0) return Outputdefinition47Enum;
+	      else if (strcmp(name,"Outputdefinition48")==0) return Outputdefinition48Enum;
          else stage=6;
    }
    if(stage==6){
-	      if (strcmp(name,"Outputdefinition48")==0) return Outputdefinition48Enum;
-	      else if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum;
+	      if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum;
 	      else if (strcmp(name,"Outputdefinition50")==0) return Outputdefinition50Enum;
 	      else if (strcmp(name,"Outputdefinition51")==0) return Outputdefinition51Enum;
@@ -752,9 +751,9 @@
 	      else if (strcmp(name,"LliboutryDuval")==0) return LliboutryDuvalEnum;
 	      else if (strcmp(name,"ExtrapolationVariable")==0) return ExtrapolationVariableEnum;
+	      else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum;
          else stage=7;
    }
    if(stage==7){
-	      if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum;
-	      else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum;
+	      if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum;
 	      else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum;
 	      else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum;
@@ -875,9 +874,9 @@
 	      else if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum;
 	      else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum;
+	      else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum;
          else stage=8;
    }
    if(stage==8){
-	      if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum;
-	      else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum;
+	      if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum;
 	      else if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum;
 	      else if (strcmp(name,"DepthAverageAnalysis")==0) return DepthAverageAnalysisEnum;
Index: /issm/trunk-jpl/src/m/classes/SMBgemb.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/SMBgemb.m	(revision 21231)
+++ /issm/trunk-jpl/src/m/classes/SMBgemb.m	(revision 21232)
@@ -21,7 +21,5 @@
 		ismelt;
 		isdensification;
-		isturbulentflux;
-        isInitialized;
-        isrestart;    
+		isturbulentflux;    
 
 		%inputs: 
@@ -39,15 +37,15 @@
 		Vz    = NaN; %height above ground at which wind (V) eas sampled [m]
         
-        % Initialization of snow properties w restart
-        Dzrst = NaN; %cell depth (m)
-        Drst = NaN; %snow density (kg m-3)
-        Rerst = NaN; %effective grain size (mm)
-        Gdnrst = NaN; %grain dendricity (0-1)
-        Gsprst = NaN; %grain sphericity (0-1)
-        ECrst = NaN; %evaporation/condensation (kg m-2)
-        Wrst = NaN; %Water content (kg m-2)
-        Arst = NaN; %albedo (0-1)
-        Trst = NaN; %snow temperature (K)
-        Sizerst = NaN; %Number of layers
+		% Initialization of snow properties
+		Dzini = NaN; %cell depth (m)
+		Dini = NaN; %snow density (kg m-3)
+		Reini = NaN; %effective grain size (mm)
+		Gdnini = NaN; %grain dendricity (0-1)
+		Gspini = NaN; %grain sphericity (0-1)
+		ECini = NaN; %evaporation/condensation (kg m-2)
+		Wini = NaN; %Water content (kg m-2)
+		Aini = NaN; %albedo (0-1)
+		Tini = NaN; %snow temperature (K)
+		Sizeini = NaN; %Number of layers
 
 		%settings: 
@@ -131,6 +129,4 @@
 		self.isdensification=1;
 		self.isturbulentflux=1;
-        self.isInitialized = true*ones(mesh.numberofelements,1);
-        self.isrestart = false*ones(mesh.numberofelements,1);
 	
 		self.aIdx = 1;
@@ -142,5 +138,5 @@
 		self.InitDensityScaling = 1.0;
 		
-        self.zMax=250*ones(mesh.numberofelements,1);
+		self.zMax=250*ones(mesh.numberofelements,1);
 		self.zMin=130*ones(mesh.numberofelements,1);
 		self.zY = 1.10*ones(mesh.numberofelements,1);
@@ -155,14 +151,16 @@
 		self.K = 7;
         
-        self.Dzrst=0.05*ones(mesh.numberofelements,1);
-        self.Drst=910.0*ones(mesh.numberofelements,1); 
-        self.Rerst=2.5*ones(mesh.numberofelements,1);
-        self.Gdnrst=0.0*ones(mesh.numberofelements,1);
-        self.Gsprst=0.0*ones(mesh.numberofelements,1);
-        self.ECrst=0.0*ones(mesh.numberofelements,1);
-        self.Wrst=0.0*ones(mesh.numberofelements,1);
-        self.Arst=self.aSnow*ones(mesh.numberofelements,1); 
-        self.Trst=273.15*ones(mesh.numberofelements,1); %*CL* Default value must be Tmean
-        self.Sizerst=200*ones(mesh.numberofelements,1);
+		self.Dzini=0.05*ones(mesh.numberofelements,2);
+		self.Dini=910.0*ones(mesh.numberofelements,2); 
+		self.Reini=2.5*ones(mesh.numberofelements,2);
+		self.Gdnini=0.0*ones(mesh.numberofelements,2);
+		self.Gspini=0.0*ones(mesh.numberofelements,2);
+		self.ECini=0.0*ones(mesh.numberofelements,1);
+		self.Wini=0.0*ones(mesh.numberofelements,2);
+		self.Aini=self.aSnow*ones(mesh.numberofelements,2); 
+		self.Tini=273.15*ones(mesh.numberofelements,2); 
+% 		/!\ Default value of Tini must be equal to Tmean but don't know Tmean yet (computed when atmospheric forcings are interpolated on mesh). 
+% 		If initialization without restart, this value will be overwritten when snow parameters are retrieved in Element.cpp 
+		self.Sizeini=2*ones(mesh.numberofelements,1);
 
 		end % }}}
@@ -178,6 +176,4 @@
 			md = checkfield(md,'fieldname','smb.isdensification','values',[0 1]);
 			md = checkfield(md,'fieldname','smb.isturbulentflux','values',[0 1]);
-            md = checkfield(md,'fieldname','smb.isInitialized','values',[false true]);
-            md = checkfield(md,'fieldname','smb.isrestart','values',[false true]);
 
 			md = checkfield(md,'fieldname','smb.Ta','timeseries',1,'NaN',1,'Inf',1,'>',273-100,'<',273+100); %-100/100 celsius min/max value
@@ -222,9 +218,4 @@
 			end
 			md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1);
-            
-            % Check that isinitialization and isrestart have different values
-            if any(self.isInitialized==self.isrestart),
-                error('isInitialized and isrestart have the same value');
-            end
 
 		end % }}}
@@ -266,15 +257,15 @@
 									'4: exponential time decay & wetness [Bougamont & Bamber, 2005]'});
                                 
-            %snow properties init
-            fielddisplay(self,'Dzrst','Initial cel depth when restart [m]');
-            fielddisplay(self,'Drst','Initial snow density when restart [kg m-3]');
-            fielddisplay(self,'Rerst','Initial grain size when restart [mm]');
-            fielddisplay(self,'Gdnrst','Initial grain dendricity when restart [-]');
-            fielddisplay(self,'Gsprst','Initial grain sphericity when restart [-]');
-            fielddisplay(self,'ECrst','Initial evaporation/condensation when restart [kg m-2]');
-            fielddisplay(self,'Wrst','Initial snow water content when restart [kg m-2]');
-            fielddisplay(self,'Arst','Initial albedo when restart [-]');
-            fielddisplay(self,'Trst','Initial snow temperature when restart [K]');
-            fielddisplay(self,'Sizerst','Initial number of layers when restart [K]');
+			%snow properties init
+			fielddisplay(self,'Dzini','Initial cel depth when restart [m]');
+			fielddisplay(self,'Dini','Initial snow density when restart [kg m-3]');
+			fielddisplay(self,'Reini','Initial grain size when restart [mm]');
+			fielddisplay(self,'Gdnini','Initial grain dendricity when restart [-]');
+			fielddisplay(self,'Gspini','Initial grain sphericity when restart [-]');
+			fielddisplay(self,'ECini','Initial evaporation/condensation when restart [kg m-2]');
+			fielddisplay(self,'Wini','Initial snow water content when restart [kg m-2]');
+			fielddisplay(self,'Aini','Initial albedo when restart [-]');
+			fielddisplay(self,'Tini','Initial snow temperature when restart [K]');
+			fielddisplay(self,'Sizeini','Initial number of layers when restart [K]');
             
             
@@ -341,8 +332,4 @@
 			WriteData(fid,prefix,'object',self,'class','smb','fieldname','denIdx','format','Integer');
 			WriteData(fid,prefix,'object',self,'class','smb','fieldname','InitDensityScaling','format','Double');
-            
-            WriteData(fid,prefix,'object',self,'class','smb','fieldname','isInitialized','format','BooleanMat','mattype',2);
-            WriteData(fid,prefix,'object',self,'class','smb','fieldname','isrestart','format','BooleanMat','mattype',2); 
-
 			WriteData(fid,prefix,'object',self,'class','smb','fieldname','outputFreq','format','Double');
 			WriteData(fid,prefix,'object',self,'class','smb','fieldname','aSnow','format','Double');
@@ -353,15 +340,15 @@
 			WriteData(fid,prefix,'object',self,'class','smb','fieldname','K','format','Double');
             
-            %snow properties init
-            WriteData(fid,prefix,'object',self,'class','smb','fieldname','Dzrst','format','DoubleMat','mattype',3);
-            WriteData(fid,prefix,'object',self,'class','smb','fieldname','Drst','format','DoubleMat','mattype',3);
-            WriteData(fid,prefix,'object',self,'class','smb','fieldname','Rerst','format','DoubleMat','mattype',3);
-            WriteData(fid,prefix,'object',self,'class','smb','fieldname','Gdnrst','format','DoubleMat','mattype',3);
-            WriteData(fid,prefix,'object',self,'class','smb','fieldname','Gsprst','format','DoubleMat','mattype',3);
-            WriteData(fid,prefix,'object',self,'class','smb','fieldname','ECrst','format','DoubleMat','mattype',2);
-            WriteData(fid,prefix,'object',self,'class','smb','fieldname','Wrst','format','DoubleMat','mattype',3);
-            WriteData(fid,prefix,'object',self,'class','smb','fieldname','Arst','format','DoubleMat','mattype',3);
-            WriteData(fid,prefix,'object',self,'class','smb','fieldname','Trst','format','DoubleMat','mattype',3);
-            WriteData(fid,prefix,'object',self,'class','smb','fieldname','Sizerst','format','IntMat','mattype',2);
+			%snow properties init
+			WriteData(fid,prefix,'object',self,'class','smb','fieldname','Dzini','format','DoubleMat','mattype',3);
+			WriteData(fid,prefix,'object',self,'class','smb','fieldname','Dini','format','DoubleMat','mattype',3);
+			WriteData(fid,prefix,'object',self,'class','smb','fieldname','Reini','format','DoubleMat','mattype',3);
+			WriteData(fid,prefix,'object',self,'class','smb','fieldname','Gdnini','format','DoubleMat','mattype',3);
+			WriteData(fid,prefix,'object',self,'class','smb','fieldname','Gspini','format','DoubleMat','mattype',3);
+			WriteData(fid,prefix,'object',self,'class','smb','fieldname','ECini','format','DoubleMat','mattype',2);
+			WriteData(fid,prefix,'object',self,'class','smb','fieldname','Wini','format','DoubleMat','mattype',3);
+			WriteData(fid,prefix,'object',self,'class','smb','fieldname','Aini','format','DoubleMat','mattype',3);
+			WriteData(fid,prefix,'object',self,'class','smb','fieldname','Tini','format','DoubleMat','mattype',3);
+			WriteData(fid,prefix,'object',self,'class','smb','fieldname','Sizeini','format','IntMat','mattype',2);
 
 			%figure out dt from forcings: 
@@ -372,8 +359,8 @@
 			WriteData(fid,prefix,'data',dt,'name','md.smb.dt','format','Double','scale',yts);
             
-            % Check if smb_dt goes evenly into transient core time step
-            if (mod(md.timestepping.time_step,dt) >= 1e-10)
+			% Check if smb_dt goes evenly into transient core time step
+			if (mod(md.timestepping.time_step,dt) >= 1e-10)
                 error('smb_dt/dt = %f. The number of SMB time steps in one transient core time step has to be an an integer',md.timestepping.time_step/dt);
-            end
+			end
 			
 			%process requested outputs