Index: /issm/trunk-jpl/jenkins/macosx_pine-island
===================================================================
--- /issm/trunk-jpl/jenkins/macosx_pine-island	(revision 23045)
+++ /issm/trunk-jpl/jenkins/macosx_pine-island	(revision 23046)
@@ -5,7 +5,7 @@
 
 #MATLAB path
-MATLAB_PATH="/Applications/MATLAB_R2017b.app/"
+MATLAB_PATH="/Applications/MATLAB_R2015b.app/"
 
-#ISSM CONFIGURATION 
+#ISSM CONFIGURATION
 ISSM_CONFIG='--prefix=$ISSM_DIR \
 	--with-matlab-dir=$MATLAB_PATH \
@@ -54,5 +54,5 @@
 #by Matlab and runme.m
 #ex: "'id',[101 102 103]"
-##                           FS                     
+##                           FS
 PYTHON_NROPTIONS=""
 MATLAB_NROPTIONS="'exclude',[701,702,703,435,IdFromString('Dakota')]"
Index: /issm/trunk-jpl/jenkins/macosx_pine-island_static
===================================================================
--- /issm/trunk-jpl/jenkins/macosx_pine-island_static	(revision 23045)
+++ /issm/trunk-jpl/jenkins/macosx_pine-island_static	(revision 23046)
@@ -5,7 +5,7 @@
 
 #MATLAB path
-MATLAB_PATH="/Applications/MATLAB_R2017b.app/"
+MATLAB_PATH="/Applications/MATLAB_R2015b.app/"
 
-#ISSM CONFIGURATION 
+#ISSM CONFIGURATION
 ISSM_CONFIG='--prefix=$ISSM_DIR \
 	--disable-static \
@@ -23,5 +23,5 @@
 	--with-m1qn3-dir=$ISSM_DIR/externalpackages/m1qn3/install \
 	--with-math77-dir=$ISSM_DIR/externalpackages/math77/install \
-	--with-fortran-lib="/usr/local/gfortran/lib/libgfortran.a /usr/local/gfortran/lib/libquadmath.a /usr/local/gfortran/lib/gcc/x86_64-apple-darwin15/6.1.0/libgcc.a" \
+	--with-fortran-lib="/usr/local/gfortran/lib/libgfortran.a /usr/local/gfortran/lib/libquadmath.a /usr/local/gfortran/lib/gcc/x86_64-apple-darwin14/5.2.0/libgcc.a" \
 	--with-numthreads=4'
 
@@ -61,5 +61,5 @@
 #by Matlab and runme.m
 #ex: "'id',[101 102 103]"
-##                           bamg mesh   FS                     
+##                           bamg mesh   FS
 #PYTHON_NROPTIONS=""
 #MATLAB_NROPTIONS="'exclude',[119,243,514,701,702,703,435,IdFromString('Dakota')]"
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 23045)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 23046)
@@ -1464,5 +1464,5 @@
 	/*Figure out maximum across the cluster: */
 	ISSM_MPI_Reduce(&maxvy,&node_maxvy,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
-	ISSM_MPI_Bcast(&node_maxvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
+	ISSM_MPI_Bcast(&node_maxvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
 	maxvy=node_maxvy;
 
@@ -1488,5 +1488,5 @@
 	/*Figure out maximum across the cluster: */
 	ISSM_MPI_Reduce(&maxvz,&node_maxvz,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
-	ISSM_MPI_Bcast(&node_maxvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
+	ISSM_MPI_Bcast(&node_maxvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
 	maxvz=node_maxvz;
 
@@ -1512,5 +1512,5 @@
 	/*Figure out minimum across the cluster: */
 	ISSM_MPI_Reduce(&minvel,&node_minvel,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
-	ISSM_MPI_Bcast(&node_minvel,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
+	ISSM_MPI_Bcast(&node_minvel,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
 	minvel=node_minvel;
 
@@ -1536,5 +1536,5 @@
 	/*Figure out minimum across the cluster: */
 	ISSM_MPI_Reduce(&minvx,&node_minvx,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
-	ISSM_MPI_Bcast(&node_minvx,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
+	ISSM_MPI_Bcast(&node_minvx,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
 	minvx=node_minvx;
 
@@ -1560,5 +1560,5 @@
 	/*Figure out minimum across the cluster: */
 	ISSM_MPI_Reduce(&minvy,&node_minvy,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
-	ISSM_MPI_Bcast(&node_minvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
+	ISSM_MPI_Bcast(&node_minvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
 	minvy=node_minvy;
 
@@ -1584,5 +1584,5 @@
 	/*Figure out minimum across the cluster: */
 	ISSM_MPI_Reduce(&minvz,&node_minvz,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
-	ISSM_MPI_Bcast(&node_minvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
+	ISSM_MPI_Bcast(&node_minvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
 	minvz=node_minvz;
 
@@ -1629,5 +1629,5 @@
 			omega_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss);
 
-			/*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ 
+			/*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
 			//J+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight;
 			J+=weight*1/2*pow(dp[0]*dp[0]+dp[1]*dp[1],2)*Jdet*gauss->weight;
@@ -1965,5 +1965,5 @@
 								IssmDouble* values    = xNewZeroInit<IssmDouble>(nlines*ncols);
 								IssmDouble* allvalues = xNew<IssmDouble>(nlines*ncols);
-								
+
 								/*Fill-in matrix*/
 								for(int j=0;j<elements->Size();j++){
@@ -1974,5 +1974,5 @@
 								ISSM_MPI_Allreduce((void*)values,(void*)allvalues,ncols*nlines,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
 								xDelete<IssmDouble>(values);
-								
+
 								if(save_results)results->AddResult(new GenericExternalResult<IssmDouble*>(results->Size()+1,output_enum,allvalues,nlines,ncols,step,time));
 								xDelete<IssmDouble>(allvalues);
@@ -2075,11 +2075,11 @@
 		case VelEnum:                            this->ElementResponsex(responses,VelEnum); break;
 		case FrictionCoefficientEnum:            NodalValuex(responses, FrictionCoefficientEnum,elements,nodes, vertices, loads, materials, parameters); break;
-		default: 
+		default:
 			if(response_descriptor_enum>=Outputdefinition1Enum && response_descriptor_enum <=Outputdefinition100Enum){
 				IssmDouble double_result = OutputDefinitionsResponsex(this,response_descriptor_enum);
 				*responses=double_result;
 			}
-			else _error_("response descriptor \"" << EnumToStringx(response_descriptor_enum) << "\" not supported yet!"); 
-			break; 
+			else _error_("response descriptor \"" << EnumToStringx(response_descriptor_enum) << "\" not supported yet!");
+			break;
 	}
 
@@ -2212,5 +2212,5 @@
 			thickness_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss);
 
-			/*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ 
+			/*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
 			J+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight;
 		}
@@ -2403,7 +2403,7 @@
 	analysis->UpdateConstraints(this);
 	delete analysis;
-	
+
 	/*Second, constraints might be time dependent: */
-	SpcNodesx(nodes,constraints,parameters,analysis_type); 
+	SpcNodesx(nodes,constraints,parameters,analysis_type);
 
 	/*Now, update degrees of freedoms: */
@@ -2416,5 +2416,5 @@
 	IssmDouble         *surface = NULL;
 	IssmDouble         *bed     = NULL;
-			
+
 	if(VerboseSolution()) _printf0_("   updating vertices positions\n");
 
@@ -2455,5 +2455,5 @@
 
 /*AMR*/
-#if !defined(_HAVE_ADOLC_) 
+#if !defined(_HAVE_ADOLC_)
 void FemModel::ReMesh(void){/*{{{*/
 
@@ -2465,5 +2465,5 @@
 	int newnumberofvertices	= -1;
 	int newnumberofelements = -1;
-	bool* my_elements			= NULL; 
+	bool* my_elements			= NULL;
 	int* my_vertices			= NULL;
 	int elementswidth       = this->GetElementsWidth();//just tria elements in this version
@@ -2471,5 +2471,5 @@
 	bool isgroundingline;
 
-	/*Branch to specific amr depending on requested method*/	
+	/*Branch to specific amr depending on requested method*/
 	this->parameters->FindParam(&amrtype,AmrTypeEnum);
 	switch(amrtype){
@@ -2484,5 +2484,5 @@
 		default: _error_("not implemented yet");
 	}
-	
+
 	/*Partitioning the new mesh. Maybe ElementsAndVerticesPartitioning.cpp could be modified to set this without iomodel.*/
 	this->ElementsAndVerticesPartitioning(newnumberofvertices,newnumberofelements,elementswidth,newelementslist,&my_elements,&my_vertices);
@@ -2515,5 +2515,5 @@
 
 		/*As the domain is 2D, it is not necessary to create nodes for this analysis*/
-		if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue;    
+		if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue;
 
 		this->CreateNodes(newnumberofvertices,my_vertices,nodecounter,analysis_enum,new_nodes);
@@ -2545,5 +2545,5 @@
 		//SetCurrentConfiguration(analysis_type);
 
-		this->analysis_counter=i;	
+		this->analysis_counter=i;
 		/*Now, plug analysis_counter and analysis_type inside the parameters: */
 		this->parameters->SetParam(this->analysis_counter,AnalysisCounterEnum);
@@ -2562,5 +2562,5 @@
 
 		ConfigureObjectsx(new_elements,this->loads,new_nodes,new_vertices,new_materials,this->parameters);
-		if(i==0){ 
+		if(i==0){
 			VerticesDofx(new_vertices,this->parameters); //only call once, we only have one set of vertices
 		}
@@ -2606,5 +2606,5 @@
 	/*Insert bedrock from mismip+ setup*/
 	/*This was used to Misomip project/simulations*/
-	
+
 	if(VerboseSolution())_printf0_("	call Mismip bedrock adjust module\n");
 
@@ -2623,5 +2623,5 @@
 			by		= 500./(1.+exp((-2./4000.)*(y-80000./2.-24000.)))+500./(1.+exp((2./4000.)*(y-80000./2.+24000.)));
 			r[i]	= max(bx+by,-720.);
-		}	
+		}
 		/*insert new bedrock*/
 		element->AddInput(BedEnum,&r[0],P1Enum);
@@ -2636,5 +2636,5 @@
 
 	if(VerboseSolution())_printf0_("	call adjust base and thickness module\n");
-	
+
 	int     numvertices = this->GetElementsWidth();
    IssmDouble rho_water,rho_ice,density,base_float;
@@ -2648,5 +2648,5 @@
 	for(int el=0;el<this->elements->Size();el++){
 		Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(el));
-	
+
 		element->GetInputListOnVertices(&s[0],SurfaceEnum);
 		element->GetInputListOnVertices(&r[0],BedEnum);
@@ -2660,14 +2660,14 @@
 			base_float = rho_ice*s[i]/(rho_ice-rho_water);
 			if(r[i]>base_float){
-				b[i] = r[i];			
-			} 
+				b[i] = r[i];
+			}
 			else {
-				b[i] = base_float;	
-			} 
+				b[i] = base_float;
+			}
 
 			if(fabs(sl[i])>0) _error_("Sea level value "<<sl[i]<<" not supported!");
 			/*update thickness and mask grounded ice level set*/
 			h[i]	  = s[i]-b[i];
-			phi[i]  = h[i]+r[i]/density;	
+			phi[i]  = h[i]+r[i]/density;
 		}
 
@@ -2677,5 +2677,5 @@
 		element->AddInput(BaseEnum,&b[0],P1Enum);
 	}
-	
+
    /*Delete*/
    xDelete<IssmDouble>(phi);
@@ -2705,5 +2705,5 @@
 	Vector<IssmDouble>* input_interpolations	= NULL;
 	IssmDouble* input_interpolations_serial	= NULL;
-   int* pos												= NULL; 
+   int* pos												= NULL;
 	IssmDouble value									= 0;
 
@@ -2725,5 +2725,5 @@
 			P0input_interp = xNew<int>(numP0inputs);
 			P1input_enums  = xNew<int>(numP1inputs);
-			P1input_interp = xNew<int>(numP1inputs);	
+			P1input_interp = xNew<int>(numP1inputs);
 		}
 		numP0inputs = 0;
@@ -2757,6 +2757,6 @@
 		}
 	}
-	
-	/*Get P0 and P1 inputs over the elements*/	
+
+	/*Get P0 and P1 inputs over the elements*/
 	pos		= xNew<int>(elementswidth);
 	vP0inputs= new Vector<IssmDouble>(numberofelements*numP0inputs);
@@ -2764,14 +2764,14 @@
 	for(int i=0;i<this->elements->Size();i++){
 		Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
-		
+
 		/*Get P0 inputs*/
 		for(int j=0;j<numP0inputs;j++){
-			TriaInput* input=xDynamicCast<TriaInput*>(element->GetInput(P0input_enums[j]));		
+			TriaInput* input=xDynamicCast<TriaInput*>(element->GetInput(P0input_enums[j]));
 			input->GetInputAverage(&value);
 			pos[0]=element->Sid()*numP0inputs+j;
-			/*Insert input in the vector*/	
+			/*Insert input in the vector*/
 			vP0inputs->SetValues(1,pos,&value,INS_VAL);
 		}
-		
+
 		/*Get P1 inputs*/
 		for(int j=0;j<numP1inputs;j++){
@@ -2780,6 +2780,6 @@
 			pos[1]=element->vertices[1]->Sid()*numP1inputs+j;
 			pos[2]=element->vertices[2]->Sid()*numP1inputs+j;
-			/*Insert input in the vector*/	
-			vP1inputs->SetValues(elementswidth,pos,input->values,INS_VAL);	
+			/*Insert input in the vector*/
+			vP1inputs->SetValues(elementswidth,pos,input->values,INS_VAL);
 		}
 	}
@@ -2790,14 +2790,14 @@
 	P0inputs=vP0inputs->ToMPISerial();
 	P1inputs=vP1inputs->ToMPISerial();
-	
+
 	/*Assign pointers*/
-	*pnumP0inputs		= numP0inputs; 
-	*pP0inputs			= P0inputs; 
+	*pnumP0inputs		= numP0inputs;
+	*pP0inputs			= P0inputs;
 	*pP0input_enums	= P0input_enums;
 	*pP0input_interp	= P0input_interp;
-	*pnumP1inputs		= numP1inputs; 
-	*pP1inputs			= P1inputs; 
+	*pnumP1inputs		= numP1inputs;
+	*pP1inputs			= P1inputs;
 	*pP1input_enums	= P1input_enums;
-	*pP1input_interp	= P1input_interp;	
+	*pP1input_interp	= P1input_interp;
 
 	/*Cleanup*/
@@ -2810,5 +2810,5 @@
 /*}}}*/
 void FemModel::InterpolateInputs(Vertices* newfemmodel_vertices,Elements* newfemmodel_elements){/*{{{*/
-	
+
 	int numberofelements			= this->elements->NumberOfElements();	//global, entire old mesh
 	int newnumberofelements		= newfemmodel_elements->Size();			//just on the new partition
@@ -2826,5 +2826,5 @@
 	int* P1input_enums  			= NULL;
 	int* P1input_interp 			= NULL;
-	IssmDouble* values			= NULL;	
+	IssmDouble* values			= NULL;
    IssmDouble* vector      	= NULL;
 	IssmDouble* x					= NULL;//global, entire old mesh
@@ -2838,5 +2838,5 @@
 	IssmDouble* newyc				= NULL;//just on the new partition
 	int* newelementslist			= NULL;//just on the new partition
-	int* sidtoindex				= NULL;//global vertices sid to partition index 
+	int* sidtoindex				= NULL;//global vertices sid to partition index
 
 	/*Get old P0 and P1  inputs (entire mesh)*/
@@ -2871,15 +2871,15 @@
 
 	/*Insert P0 and P1 inputs into the new elements (just on the new partition)*/
-	values=xNew<IssmDouble>(elementswidth);	
+	values=xNew<IssmDouble>(elementswidth);
 	for(int i=0;i<newfemmodel_elements->Size();i++){//just on the new partition
 		Element* element=xDynamicCast<Element*>(newfemmodel_elements->GetObjectByOffset(i));
 		/*newP0inputs is just on the new partition*/
 		for(int j=0;j<numP0inputs;j++){
-			switch(P0input_interp[j]){	
+			switch(P0input_interp[j]){
 				case P0Enum:
 				case DoubleInputEnum:
 					element->AddInput(new DoubleInput(P0input_enums[j],newP0inputs[i*numP0inputs+j]));
 					break;
-				case IntInputEnum: 
+				case IntInputEnum:
 					element->AddInput(new IntInput(P0input_enums[j],reCast<int>(newP0inputs[i*numP0inputs+j])));
 					break;
@@ -2899,5 +2899,5 @@
 		}
 	}
-	
+
 	/*Cleanup*/
 	xDelete<IssmDouble>(P0inputs);
@@ -2938,5 +2938,5 @@
 
 	if(!this->elements || !this->vertices || !this->results || !this->parameters) return;
-	 
+
 	parameters->FindParam(&step,StepEnum);
 	parameters->FindParam(&time,TimeEnum);
@@ -2956,5 +2956,5 @@
 	this->results->AddResult(new GenericExternalResult<IssmDouble*>(this->results->Size()+1,MeshYEnum,
 					y,numberofvertices,1,step,time));
-	
+
 	/*Cleanup*/
 	xDelete<IssmDouble>(x);
@@ -3017,5 +3017,5 @@
    /*Assemble*/
 	vmasklevelset->Assemble();
-	
+
 	/*Serialize and set output*/
 	(*pmasklevelset)=vmasklevelset->ToMPISerial();
@@ -3031,5 +3031,5 @@
 
 	/*newelementslist is in Matlab indexing*/
-	
+
 	/*Creating connectivity table*/
 	int* connectivity=NULL;
@@ -3042,10 +3042,10 @@
 			connectivity[vertexid-1]+=1;//Matlab to C indexing
 		}
-	}	
+	}
 
 	/*Create vertex and insert in vertices*/
 	for(int i=0;i<newnumberofvertices;i++){
 		if(my_vertices[i]){
-			Vertex *newvertex=new Vertex();	
+			Vertex *newvertex=new Vertex();
 			newvertex->id=i+1;
 			newvertex->sid=i;
@@ -3058,6 +3058,6 @@
 			newvertex->connectivity=connectivity[i];
 			newvertex->clone=false;//itapopo check this
-			vertices->AddObject(newvertex);	
-		} 
+			vertices->AddObject(newvertex);
+		}
 	}
 
@@ -3086,5 +3086,5 @@
 			}
 			else newtria->element_type_list=NULL;
-			
+
 			/*Element hook*/
 			int matpar_id=newnumberofelements+1; //retrieve material parameter id (last pointer in femodel->materials)
@@ -3092,5 +3092,5 @@
 			/*retrieve vertices ids*/
 			int* vertex_ids=xNew<int>(elementswidth);
-			for(int j=0;j<elementswidth;j++)	vertex_ids[j]=reCast<int>(newelementslist[elementswidth*i+j]);//this Hook wants Matlab indexing	
+			for(int j=0;j<elementswidth;j++)	vertex_ids[j]=reCast<int>(newelementslist[elementswidth*i+j]);//this Hook wants Matlab indexing
 			/*Setting the hooks*/
 			newtria->numanalyses =this->nummodels;
@@ -3104,6 +3104,6 @@
 			/*Clean up*/
 			xDelete<int>(vertex_ids);
-			elements->AddObject(newtria);	
-		} 
+			elements->AddObject(newtria);
+		}
 	}
 
@@ -3115,12 +3115,12 @@
 	for(int i=0;i<newnumberofelements;i++){
 		if(my_elements[i]){
-			materials->AddObject(new Matice(i+1,i,MaticeEnum));	
-		} 
-	}
-	
+			materials->AddObject(new Matice(i+1,i,MaticeEnum));
+		}
+	}
+
 	/*Add new constant material property to materials, at the end: */
 	Matpar *newmatpar=static_cast<Matpar*>(this->materials->GetObjectByOffset(this->materials->Size()-1)->copy());
 	newmatpar->SetMid(newnumberofelements+1);
-	materials->AddObject(newmatpar);//put it at the end of the materials	    
+	materials->AddObject(newmatpar);//put it at the end of the materials
 }
 /*}}}*/
@@ -3129,8 +3129,8 @@
 	int lid=0;
 	for(int j=0;j<newnumberofvertices;j++){
-		if(my_vertices[j]){				
-			
-			Node* newnode=new Node();	
-			
+		if(my_vertices[j]){
+
+			Node* newnode=new Node();
+
 			/*id: */
 			newnode->id=nodecounter+j+1;
@@ -3138,12 +3138,12 @@
 			newnode->lid=lid++;
 			newnode->analysis_enum=analysis_enum;
-			
+
 			/*Initialize coord_system: Identity matrix by default*/
 			for(int k=0;k<3;k++) for(int l=0;l<3;l++) newnode->coord_system[k][l]=0.0;
 			for(int k=0;k<3;k++) newnode->coord_system[k][k]=1.0;
-			
+
 			/*indexing:*/
 			newnode->indexingupdate=true;
-			
+
 			Analysis* analysis=EnumToAnalysis(analysis_enum);
 			int *doftypes=NULL;
@@ -3159,5 +3159,5 @@
 			/*Stressbalance Horiz*/
 			if(analysis_enum==StressbalanceAnalysisEnum){
-				// itapopo this code is rarely used. 
+				// itapopo this code is rarely used.
 				/*Coordinate system provided, convert to coord_system matrix*/
 				//XZvectorsToCoordinateSystem(&this->coord_system[0][0],&iomodel->Data(StressbalanceReferentialEnum)[j*6]);
@@ -3174,5 +3174,5 @@
 	if(!femmodel_vertices) _error_("GetMesh: vertices are NULL.");
 	if(!femmodel_elements) _error_("GetMesh: elements are NULL.");
-	
+
 	int numberofvertices = femmodel_vertices->NumberOfVertices();
 	int numberofelements = femmodel_elements->NumberOfElements();
@@ -3180,5 +3180,5 @@
 	IssmDouble* x			= NULL;
 	IssmDouble* y			= NULL;
-	IssmDouble* z			= NULL;	
+	IssmDouble* z			= NULL;
 	int* elementslist 	= NULL;
 	int* elem_vertices	= NULL;
@@ -3189,5 +3189,5 @@
 	/*Get vertices coordinates*/
 	VertexCoordinatesx(&x,&y,&z,femmodel_vertices,false) ;
-	
+
 	/*Get element vertices*/
 	elem_vertices				= xNew<int>(elementswidth);
@@ -3204,5 +3204,5 @@
     	vid3->SetValue(element->sid,elem_vertices[2],INS_VAL);
    }
-		
+
 	/*Assemble*/
    vid1->Assemble();
@@ -3214,5 +3214,5 @@
    id2 = vid2->ToMPISerial();
 	id3 = vid3->ToMPISerial();
-	
+
 	/*Construct elements list*/
 	elementslist=xNew<int>(numberofelements*elementswidth);
@@ -3223,5 +3223,5 @@
 		elementslist[elementswidth*i+2] = reCast<int>(id3[i])+1; //InterpMesh wants Matlab indexing
 	}
-	
+
 	/*Assign pointers*/
 	*px				= x;
@@ -3244,5 +3244,5 @@
 	if(!femmodel_vertices) _error_("GetMesh: vertices are NULL.");
 	if(!femmodel_elements) _error_("GetMesh: elements are NULL.");
-	
+
 	int numberofvertices			= femmodel_vertices->Size();	//number of vertices of this partition
 	int numbertotalofvertices	= femmodel_vertices->NumberOfVertices();	//number total of vertices (entire mesh)
@@ -3251,9 +3251,9 @@
 	IssmDouble* x					= NULL;
 	IssmDouble* y					= NULL;
-	IssmDouble* z					= NULL;	
+	IssmDouble* z					= NULL;
 	int* elementslist				= NULL;
 	int* sidtoindex				= NULL;
 	int* elem_vertices			= NULL;
-	
+
 	/*Get vertices coordinates of this partition*/
 	sidtoindex	= xNewZeroInit<int>(numbertotalofvertices);//entire mesh, all vertices
@@ -3261,8 +3261,8 @@
 	y				= xNew<IssmDouble>(numberofvertices);//just this partitio;
 	z				= xNew<IssmDouble>(numberofvertices);//just this partitio;
-	
+
 	/*Go through in this partition (vertices)*/
 	for(int i=0;i<numberofvertices;i++){//just this partition
-		Vertex* vertex=(Vertex*)femmodel_vertices->GetObjectByOffset(i);	
+		Vertex* vertex=(Vertex*)femmodel_vertices->GetObjectByOffset(i);
 		/*Attention: no spherical coordinates*/
 		x[i]=vertex->GetX();
@@ -3277,5 +3277,5 @@
 	elementslist = xNew<int>(numberofelements*elementswidth);
 	if(numberofelements*elementswidth<0) _error_("numberofelements negative.");
-	
+
 	for(int i=0;i<numberofelements;i++){//just this partition
     	Element* element=xDynamicCast<Element*>(femmodel_elements->GetObjectByOffset(i));
@@ -3284,6 +3284,6 @@
 		elementslist[elementswidth*i+1] = sidtoindex[elem_vertices[1]]+1; //InterpMesh wants Matlab indexing
 		elementslist[elementswidth*i+2] = sidtoindex[elem_vertices[2]]+1; //InterpMesh wants Matlab indexing
-	}	
-		
+	}
+
 	/*Assign pointers*/
 	*px				= x;
@@ -3291,5 +3291,5 @@
 	*pz				= z;
 	*pelementslist = elementslist; //Matlab indexing. InterMesh uses this type.
-	*psidtoindex	= sidtoindex;  //it is ncessary to insert inputs 
+	*psidtoindex	= sidtoindex;  //it is ncessary to insert inputs
 
 	/*Cleanup*/
@@ -3302,7 +3302,7 @@
 	/*OTHERS CONSTRAINTS MUST BE IMPLEMENTED*/
 	if(analysis_enum!=StressbalanceAnalysisEnum) return;
-	
+
 	int numberofnodes_analysistype= this->nodes->NumberOfNodes(analysis_enum);
-	int dofpernode						= 2;														//vx and vy 
+	int dofpernode						= 2;														//vx and vy
 	int numberofcols					= dofpernode*2;										//to keep dofs and flags in the vspc vector
 	int numberofvertices				= this->vertices->NumberOfVertices();			//global, entire old mesh
@@ -3328,5 +3328,5 @@
 	newy=xNew<IssmDouble>(newnumberofvertices);//just the new partition
 	for(int i=0;i<newnumberofvertices;i++){//just the new partition
-		Vertex* vertex=(Vertex*)newfemmodel_vertices->GetObjectByOffset(i);	
+		Vertex* vertex=(Vertex*)newfemmodel_vertices->GetObjectByOffset(i);
 		/*Attention: no spherical coordinates*/
 		newx[i]=vertex->GetX();
@@ -3336,5 +3336,5 @@
 	/*Get spcvx and spcvy of old mesh*/
 	for(int i=0;i<this->constraints->Size();i++){
-		
+
 		Constraint* constraint=(Constraint*)constraints->GetObjectByOffset(i);
 		if(!constraint->InAnalysis(analysis_enum)) _error_("AMR create constraints for "<<EnumToStringx(analysis_enum)<<" not supported yet!\n");
@@ -3343,12 +3343,12 @@
 		int dof					= spcstatic->GetDof();
 		int node					= spcstatic->GetNodeId();
-		IssmDouble spcvalue	= spcstatic->GetValue(); 
+		IssmDouble spcvalue	= spcstatic->GetValue();
 		int nodeindex			= node-1;
-		
+
 		/*vx and vx flag insertion*/
 		if(dof==0) {//vx
 			vspc->SetValue(nodeindex*numberofcols,spcvalue,INS_VAL);    //vx
 			vspc->SetValue(nodeindex*numberofcols+dofpernode,1,INS_VAL);//vxflag
-		} 
+		}
 		/*vy and vy flag insertion*/
 		if(dof==1){//vy
@@ -3366,5 +3366,5 @@
 								spc,numberofvertices,numberofcols,
 								newx,newy,newnumberofvertices,NULL);
-	
+
 	/*Now, insert the interpolated constraints in the data set (constraints)*/
 	count=0;
@@ -3383,5 +3383,5 @@
 		/*spcvy*/
 		if(!xIsNan<IssmDouble>(newspc[i*numberofcols+1]) && newspc[i*numberofcols+dofpernode+1]>(1-eps) ){
-			newfemmodel_constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+vertex->Sid()+1,1,newspc[i*numberofcols+1],analysis_enum)); 
+			newfemmodel_constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+vertex->Sid()+1,1,newspc[i*numberofcols+1],analysis_enum));
 			//add count'th spc, on node i+1, setting dof 1 to vx.
 			count++;
@@ -3440,11 +3440,11 @@
 	bool *my_elements = NULL;
 	int *my_vertices  = NULL;
-	
-	_assert_(newnumberofvertices>0); 
-	_assert_(newnumberofelements>0); 
+
+	_assert_(newnumberofvertices>0);
+	_assert_(newnumberofelements>0);
 	epart=xNew<int>(newnumberofelements);
 	npart=xNew<int>(newnumberofvertices);
    index=xNew<int>(elementswidth*newnumberofelements);
-   
+
 	for (int i=0;i<newnumberofelements;i++){
    	for (int j=0;j<elementswidth;j++){
@@ -3466,5 +3466,5 @@
 		for (int i=0;i<newnumberofvertices;i++) npart[i]=0;
 	}
-	else _error_("At least one processor is required");	    
+	else _error_("At least one processor is required");
 
 	my_vertices=xNew<int>(newnumberofvertices);
@@ -3476,12 +3476,12 @@
 	for(int i=0;i<newnumberofelements;i++){
 		/*!All elements have been partitioned above, only deal with elements for this cpu: */
-		if(my_rank==epart[i]){ 
+		if(my_rank==epart[i]){
 			my_elements[i]=true;
-			/*Now that we are here, we can also start building the list of vertices belonging to this cpu partition: we use 
-			 *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) 
-			 into the vertices coordinates. If we start plugging 1 into my_vertices for each index[n][i] (i=0:2), then my_vertices 
+			/*Now that we are here, we can also start building the list of vertices belonging to this cpu partition: we use
+			 *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing)
+			 into the vertices coordinates. If we start plugging 1 into my_vertices for each index[n][i] (i=0:2), then my_vertices
 			 will hold which vertices belong to this partition*/
 			for(int j=0;j<elementswidth;j++){
-				_assert_(newelementslist[elementswidth*i+j]-1<newnumberofvertices);//newelementslist is in Matlab indexing 
+				_assert_(newelementslist[elementswidth*i+j]-1<newnumberofvertices);//newelementslist is in Matlab indexing
 				my_vertices[newelementslist[elementswidth*i+j]-1]=1;//newelementslist is in Matlab indexing
 			}
@@ -3495,16 +3495,16 @@
 	/*Free ressources:*/
 	xDelete<int>(epart);
-	xDelete<int>(npart);	    
+	xDelete<int>(npart);
 	xDelete<int>(index);
 }
 /*}}}*/
 void FemModel::SmoothedDeviatoricStressTensor(IssmDouble** ptauxx,IssmDouble** ptauyy,IssmDouble** ptauxy){/*{{{*/
-	
+
 	int elementswidth							= this->GetElementsWidth();//just 2D mesh, tria elements
    int numberofvertices						= this->vertices->NumberOfVertices();
    IssmDouble weight 						= 0.;
-	IssmDouble*	tauxx							= NULL; 
-	IssmDouble*	tauyy							= NULL; 
-	IssmDouble*	tauxy							= NULL; 
+	IssmDouble*	tauxx							= NULL;
+	IssmDouble*	tauyy							= NULL;
+	IssmDouble*	tauxy							= NULL;
    IssmDouble* totalweight 				= NULL;
 	IssmDouble* deviatoricstressxx 		= xNew<IssmDouble>(elementswidth);
@@ -3516,8 +3516,8 @@
    Vector<IssmDouble>* vectauxy			= new Vector<IssmDouble>(numberofvertices);
    Vector<IssmDouble>* vectotalweight	= new Vector<IssmDouble>(numberofvertices);
-	
+
 	/*Update the Deviatoric Stress tensor over the elements*/
 	this->DeviatoricStressx();
-	
+
    /*Calculate the Smoothed Deviatoric Stress tensor*/
 	for(int i=0;i<this->elements->Size();i++){
@@ -3564,5 +3564,5 @@
 	/*Divide for the total weight*/
 	for(int i=0;i<numberofvertices;i++){
-		_assert_(totalweight[i]>0);	
+		_assert_(totalweight[i]>0);
 		tauxx[i] = tauxx[i]/totalweight[i];
 		tauyy[i] = tauyy[i]/totalweight[i];
@@ -3589,5 +3589,5 @@
 void FemModel::ZZErrorEstimator(IssmDouble** pelementerror){/*{{{*/
 
-	/*Compute the Zienkiewicz and Zhu (ZZ) error estimator for the deviatoric stress tensor. 
+	/*Compute the Zienkiewicz and Zhu (ZZ) error estimator for the deviatoric stress tensor.
 	 * Ref.: Zienkiewicz and Zhu, A Simple Error Estimator and Adaptive Procedure for Practical Engineering Analysis, Int. J. Numer. Meth. Eng, 1987*/
 
@@ -3609,5 +3609,5 @@
 	/*Get smoothed deviatoric stress tensor*/
 	this->SmoothedDeviatoricStressTensor(&smoothedtauxx,&smoothedtauyy,&smoothedtauxy);
-	
+
 	/*Integrate the error over elements*/
    for(int i=0;i<this->elements->Size();i++){
@@ -3617,5 +3617,5 @@
       element->GetInputListOnVertices(tauxy,DeviatoricStressxyEnum);
       element->GetVerticesSidList(elem_vertices);
-		
+
 		/*Integrate*/
 		element->GetVerticesCoordinates(&xyz_list);
@@ -3632,7 +3632,7 @@
 				ftxy+=(tauxy[n]-smoothedtauxy[elem_vertices[n]])*basis[n];
 			}
-			error+=Jdet*gauss->weight*( pow(ftxx,2)+pow(ftyy,2)+pow(ftxy,2) ); //e^2 
-		}
-		/*Set the error in the global vector*/	
+			error+=Jdet*gauss->weight*( pow(ftxx,2)+pow(ftyy,2)+pow(ftxy,2) ); //e^2
+		}
+		/*Set the error in the global vector*/
       sid=element->Sid();
 		error = sqrt(error);//sqrt(e^2)
@@ -3648,5 +3648,5 @@
    /*Serialize and set output*/
    (*pelementerror)=velementerror->ToMPISerial();
-	
+
 	/*Cleanup*/
 	xDelete<IssmDouble>(smoothedtauxx);
@@ -3692,5 +3692,5 @@
       Tria* triaelement = xDynamicCast<Tria*>(element);
       weight            = triaelement->GetArea();//the tria area is a choice for the weight
-      
+
 		/*dH/dx*/
       vecdHdx->SetValue(elem_vertices[0],weight*GradH[0],ADD_VAL);
@@ -3760,5 +3760,5 @@
    /*Get smoothed deviatoric stress tensor*/
    this->SmoothedGradThickness(&smoothed_dHdx,&smoothed_dHdy);
-   
+
 	/*Integrate the error over elements*/
    for(int i=0;i<this->elements->Size();i++){
@@ -3846,5 +3846,5 @@
 	IssmDouble* x					= NULL;
 	IssmDouble* y					= NULL;
-	IssmDouble* z					= NULL;	
+	IssmDouble* z					= NULL;
 	IssmDouble* xyz_list			= NULL;
 	IssmDouble x1,y1,x2,y2,x3,y3;
@@ -3855,5 +3855,5 @@
       //element->GetVerticesSidList(elem_vertices);
       int sid = element->Sid();
-		element->GetVerticesCoordinates(&xyz_list); 
+		element->GetVerticesCoordinates(&xyz_list);
 		x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1];
 		x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1];
@@ -3888,9 +3888,9 @@
 		_error_("level set type not implemented yet!");
 	}
-	
+
 	/*Outputs*/
 	IssmDouble* zerolevelset_points			= NULL;
 	int npoints										= 0;
-	
+
 	/*Intermediaries*/
  	int elementswidth                   	= this->GetElementsWidth();
@@ -3905,5 +3905,5 @@
 	int count,sid;
 	IssmDouble xc,yc,x1,y1,x2,y2,x3,y3;
-	
+
 	/*Use the element center coordinate if level set is zero (grounding line or ice front), otherwise set NAN*/
    for(int i=0;i<this->elements->Size();i++){
@@ -3912,13 +3912,13 @@
 		element->GetVerticesSidList(elem_vertices);
 		sid= element->Sid();
-		element->GetVerticesCoordinates(&xyz_list); 
+		element->GetVerticesCoordinates(&xyz_list);
 		x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1];
 		x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1];
 		x3 = xyz_list[3*2+0];y3 = xyz_list[3*2+1];
 		xc	= NAN;
-		yc	= NAN;	
+		yc	= NAN;
      	Tria* tria 	= xDynamicCast<Tria*>(element);
 		if(tria->IsIceInElement()){/*verify if there is ice in the element*/
-			if(levelset[0]*levelset[1]<0. || levelset[0]*levelset[2]<0. ||	
+			if(levelset[0]*levelset[1]<0. || levelset[0]*levelset[2]<0. ||
 				abs(levelset[0]*levelset[1])<DBL_EPSILON || abs(levelset[0]*levelset[2])<DBL_EPSILON) {
 				xc=(x1+x2+x3)/3.;
@@ -3950,5 +3950,5 @@
 		}
 	}
-	
+
 	/*Assign outputs*/
 	numberofpoints				= npoints;
@@ -3990,5 +3990,5 @@
 	responses_pointer=d_responses;
 
-	//watch out, we have more d_numresponses than numresponsedescriptors, because the responses have been expanded if they were scaled. 
+	//watch out, we have more d_numresponses than numresponsedescriptors, because the responses have been expanded if they were scaled.
 	//because we don't know the d_responses descriptors (the scaled ones) we can't key off them, so we will key off the responses_descriptors: */
 
@@ -4083,8 +4083,8 @@
 
 	int         ns,nsmax;
-	
+
 	/*Go through elements, and add contribution from each element to the deflection vector wg:*/
 	ns = elements->Size();
-	
+
 	/*Figure out max of ns: */
 	ISSM_MPI_Reduce(&ns,&nsmax,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
@@ -4105,5 +4105,5 @@
 		}
 	}
-	
+
 	/*One last time: */
 	pUp->Assemble();
@@ -4124,9 +4124,9 @@
 
 	int         ns,nsmax;
-	
+
 	/*Go through elements, and add contribution from each element to the deflection vector wg:*/
 	ns = elements->Size();
-	
-	/*First, figure out the surface area of Earth: */ 
+
+	/*First, figure out the surface area of Earth: */
 	for(int i=0;i<ns;i++){
 		Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
@@ -4152,5 +4152,5 @@
 		}
 	}
-	
+
 	/*One last time: */
 	pUp->Assemble();
@@ -4214,5 +4214,5 @@
 	}
 	if(VerboseConvergence())_printf0_("\n");
-		
+
 	/*One last time: */
 	pRSLgi->Assemble();
@@ -4232,10 +4232,10 @@
 	/*serialized vectors:*/
 	IssmDouble* RSLg_old=NULL;
-	
+
 	IssmDouble  eartharea=0;
 	IssmDouble  eartharea_cpu=0;
 
 	int         ns,nsmax;
-	
+
 	/*Serialize vectors from previous iteration:*/
 	RSLg_old=pRSLg_old->ToMPISerial();
@@ -4249,5 +4249,5 @@
 		eartharea_cpu += element->GetAreaSpherical();
 	}
-	
+
 	ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
 	ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
@@ -4267,5 +4267,5 @@
 	}
 	if(verboseconvolution)if(VerboseConvergence())_printf0_("\n");
-	
+
 	/*Free ressources:*/
 	xDelete<IssmDouble>(RSLg_old);
@@ -4355,6 +4355,6 @@
 	/*Free ressources:*/
 	xDelete<IssmDouble>(RSLg_old);
-	
-	
+
+
 }
 /*}}}*/
@@ -4363,10 +4363,10 @@
 	/*serialized vectors:*/
 	IssmDouble* RSLg=NULL;
-	
+
 	IssmDouble  eartharea=0;
 	IssmDouble  eartharea_cpu=0;
 
 	int         ns,nsmax;
-	
+
 	/*Serialize vectors from previous iteration:*/
 	RSLg=pRSLg->ToMPISerial();
@@ -4374,5 +4374,5 @@
 	/*Go through elements, and add contribution from each element to the deflection vector wg:*/
 	ns = elements->Size();
-	
+
 	/*First, figure out the area of the ocean, which is needed to compute the eustatic component: */
 	for(int i=0;i<ns;i++){
@@ -4402,5 +4402,5 @@
 		}
 	}
-	
+
 	/*One last time: */
 	pUp->Assemble();
@@ -4436,5 +4436,5 @@
 	ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
 	ISSM_MPI_Bcast(&oceanarea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-	
+
 	ISSM_MPI_Reduce (&oceanvalue_cpu,&oceanvalue,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
 	ISSM_MPI_Bcast(&oceanvalue,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
@@ -4442,5 +4442,5 @@
 	/*Free ressources:*/
 	xDelete<IssmDouble>(RSLg_serial);
-	
+
 	return oceanvalue/oceanarea;
 }
@@ -4458,5 +4458,5 @@
 	int*                eplzigzag_counter =	NULL;
 	int                 eplflip_lock;
-	
+
 	HydrologyDCEfficientAnalysis* effanalysis =  new HydrologyDCEfficientAnalysis();
 	HydrologyDCInefficientAnalysis* inefanalysis =  new HydrologyDCInefficientAnalysis();
@@ -4465,8 +4465,8 @@
 	mask=new Vector<IssmDouble>(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
 	recurence=new Vector<IssmDouble>(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
-	this->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum); 
-	this->parameters->FindParam(&eplflip_lock,HydrologydcEplflipLockEnum); 
+	this->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum);
+	this->parameters->FindParam(&eplflip_lock,HydrologydcEplflipLockEnum);
 	GetVectorFromInputsx(&old_active,this,HydrologydcMaskEplactiveNodeEnum,NodeSIdEnum);
-	
+
 	for (int i=0;i<elements->Size();i++){
 		Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
@@ -4486,6 +4486,6 @@
 	/*Assemble and serialize*/
 	mask->Assemble();
-	serial_mask=mask->ToMPISerial();	
-	
+	serial_mask=mask->ToMPISerial();
+
 	xDelete<int>(eplzigzag_counter);
 	xDelete<IssmDouble>(serial_rec);
@@ -4529,5 +4529,5 @@
 	int sum_counter;
 	ISSM_MPI_Reduce(&counter,&sum_counter,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-	ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm());                
+	ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm());
 	counter=sum_counter;
 	*pEplcount = counter;
@@ -4650,5 +4650,5 @@
 	int sum_counter;
 	ISSM_MPI_Reduce(&counter,&sum_counter,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-	ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm());                
+	ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm());
 	counter=sum_counter;
 	*pL2count = counter;
@@ -4735,5 +4735,5 @@
 }
 /*}}}*/
-#ifdef _HAVE_JAVASCRIPT_ 
+#ifdef _HAVE_JAVASCRIPT_
 FemModel::FemModel(IssmDouble* buffer, int buffersize, char* toolkits, char* solution, char* modelname,ISSM_MPI_Comm incomm, bool trace){ /*{{{*/
 	/*configuration: */
@@ -4750,10 +4750,10 @@
 	/*From command line arguments, retrieve different filenames needed to create the FemModel: */
 	solution_type=StringToEnumx(solution);
-	
+
 	/*Create femmodel from input files: */
 	profiler->Start(MPROCESSOR);
 	this->InitFromBuffers((char*)buffer,buffersize,toolkits, solution_type,trace,NULL);
 	profiler->Stop(MPROCESSOR);
-	
+
 	/*Save communicator in the parameters dataset: */
 	this->parameters->AddObject(new GenericParam<ISSM_MPI_Comm>(incomm,FemModelCommEnum));
@@ -4770,5 +4770,5 @@
 	size_t* poutputbuffersize;
 
-	
+
 	/*Before we delete the profiler, report statistics for this run: */
 	profiler->Stop(TOTAL);  //final tagging
@@ -4783,5 +4783,5 @@
 				);
 	_printf0_("\n");
-	
+
 	/*Before we close the output file, recover the buffer and size:*/
 	outputbufferparam = xDynamicCast<GenericParam<char**>*>(this->parameters->FindParamObject(OutputBufferPointerEnum));
@@ -4823,6 +4823,6 @@
 
 	/*Open output file once for all and add output file descriptor to parameters*/
-	output_fid=open_memstream(&outputbuffer,&outputsize); 
-	if(output_fid==NULL||output_fid==0)_error_("could not initialize output stream");
+	output_fid=open_memstream(&outputbuffer,&outputsize);
+	if(output_fid==NULL)_error_("could not initialize output stream");
 	this->parameters->SetParam(output_fid,OutputFilePointerEnum);
 	this->parameters->AddObject(new GenericParam<char**>(&outputbuffer,OutputBufferPointerEnum));
@@ -4865,5 +4865,5 @@
 		/*Initialize hmaxvertices with NAN*/
 		hmaxvertices_serial=xNew<IssmDouble>(numberofvertices);
-		for(int i=0;i<numberofvertices;i++) hmaxvertices_serial[i]=NAN; 
+		for(int i=0;i<numberofvertices;i++) hmaxvertices_serial[i]=NAN;
 		/*Fill hmaxvertices*/
 		if(this->amrbamg->groundingline_distance>0)		this->GethmaxVerticesFromZeroLevelSetDistance(hmaxvertices_serial,MaskGroundediceLevelsetEnum);
@@ -4872,5 +4872,5 @@
 		if(this->amrbamg->deviatoricerror_threshold>0)	this->GethmaxVerticesFromEstimators(hmaxvertices_serial,DeviatoricStressErrorEstimatorEnum);
 	}
-	
+
 	if(my_rank==0){
 		this->amrbamg->ExecuteRefinementBamg(vector_serial,hmaxvertices_serial,&newnumberofvertices,&newnumberofelements,&newx,&newy,&newz,&newelementslist);
@@ -4892,8 +4892,8 @@
 		newelementslist=xNew<int>(newnumberofelements*this->GetElementsWidth());
 	}
-	ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());	
-	ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());	
-	ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());	
-	ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm());	
+	ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+	ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+	ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+	ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm());
 
 	/*Assign output pointers*/
@@ -4928,5 +4928,5 @@
 	/*Create bamg data structures for bamg*/
 	this->amrbamg = new AmrBamg();
-	
+
 	/*Get amr parameters*/
 	this->parameters->FindParam(&hmin,AmrHminEnum);
@@ -4952,5 +4952,5 @@
 
 	/*Re-create original mesh and put it in bamg structure (only cpu 0)*/
-	if(my_rank==0){ 
+	if(my_rank==0){
 		this->amrbamg->Initialize(elements,x,y,numberofvertices,numberofelements);
 	}
@@ -4966,5 +4966,5 @@
 
 	if(!hmaxvertices) _error_("hmaxvertices is NULL!\n");
-	
+
 	/*Intermediaries*/
 	int numberofvertices			 = this->vertices->NumberOfVertices();
@@ -4973,5 +4973,5 @@
 
 	switch(levelset_type){
-		case MaskGroundediceLevelsetEnum: 
+		case MaskGroundediceLevelsetEnum:
 			threshold	= this->amrbamg->groundingline_distance;
 			resolution	= this->amrbamg->groundingline_resolution;
@@ -4987,5 +4987,5 @@
 	this->GetVerticeDistanceToZeroLevelSet(&verticedistance,levelset_type);
 	if(!verticedistance) _error_("verticedistance is NULL!\n");
-	
+
 	/*Fill hmaxVertices*/
 	for(int i=0;i<numberofvertices;i++){
@@ -5001,5 +5001,5 @@
 /*}}}*/
 void FemModel::GethmaxVerticesFromEstimators(IssmDouble* hmaxvertices,int errorestimator_type){/*{{{*/
-   
+
 	if(!hmaxvertices) _error_("hmaxvertices is NULL!\n");
 
@@ -5009,5 +5009,5 @@
 	int numberofvertices						= this->vertices->NumberOfVertices();
 	IssmDouble* maxlength					= xNew<IssmDouble>(numberofelements);
-	IssmDouble* error_vertices				= xNewZeroInit<IssmDouble>(numberofvertices);	
+	IssmDouble* error_vertices				= xNewZeroInit<IssmDouble>(numberofvertices);
 	IssmDouble* error_elements				= NULL;
 	IssmDouble* x								= NULL;
@@ -5022,5 +5022,5 @@
 	/*Fill variables*/
 	switch(errorestimator_type){
-		case ThicknessErrorEstimatorEnum: 
+		case ThicknessErrorEstimatorEnum:
 			threshold		= this->amrbamg->thicknesserror_threshold;
 			groupthreshold	= this->amrbamg->thicknesserror_groupthreshold;
@@ -5047,10 +5047,10 @@
       	case ThicknessErrorEstimatorEnum:			this->amrbamg->thicknesserror_maximum 	= maxerror;break;
       	case DeviatoricStressErrorEstimatorEnum: 	this->amrbamg->deviatoricerror_maximum = maxerror;break;
-   	}	
+   	}
 	}
 
 	/*Get mesh*/
 	this->GetMesh(this->vertices,this->elements,&x,&y,&z,&index);
-	
+
 	/*Fill error_vertices (this is the sum of all elements connected to the vertex)*/
 	for(int i=0;i<numberofelements;i++){
@@ -5066,5 +5066,5 @@
 		error_vertices[v2]+=error_elements[i];
 		error_vertices[v3]+=error_elements[i];
-	}	
+	}
 
 	/*Fill hmaxvertices with the criteria*/
@@ -5080,5 +5080,5 @@
 			}
 		}
-		/*Now, fill the hmaxvertices if requested*/	  
+		/*Now, fill the hmaxvertices if requested*/
 		if(refine){
 			for(int j=0;j<elementswidth;j++){
@@ -5110,5 +5110,5 @@
 	/*Output*/
 	IssmDouble* verticedistance;
-	
+
 	/*Intermediaries*/
    int numberofvertices       = this->vertices->NumberOfVertices();
@@ -5122,6 +5122,6 @@
 	/*Get vertices coordinates*/
 	VertexCoordinatesx(&x,&y,&z,this->vertices,false) ;
-	
-	/*Get points which level set is zero (center of elements with zero level set)*/	
+
+	/*Get points which level set is zero (center of elements with zero level set)*/
 	this->GetZeroLevelSetPoints(&levelset_points,numberofpoints,levelset_type);
 
@@ -5132,7 +5132,7 @@
 		for(int j=0;j<numberofpoints;j++){
 			distance=sqrt((x[i]-levelset_points[2*j])*(x[i]-levelset_points[2*j])+(y[i]-levelset_points[2*j+1])*(y[i]-levelset_points[2*j+1]));
-			verticedistance[i]=min(distance,verticedistance[i]);		
-		}
-	}	
+			verticedistance[i]=min(distance,verticedistance[i]);
+		}
+	}
 
 	/*Assign the pointer*/
@@ -5168,10 +5168,10 @@
 	if(this->amr->groundingline_distance>0)		this->GetElementDistanceToZeroLevelSet(&gl_distance,MaskGroundediceLevelsetEnum);
    if(this->amr->icefront_distance>0)				this->GetElementDistanceToZeroLevelSet(&if_distance,MaskIceLevelsetEnum);
-   if(this->amr->thicknesserror_threshold>0)		this->ThicknessZZErrorEstimator(&thicknesserror);	
-	if(this->amr->deviatoricerror_threshold>0)	this->ZZErrorEstimator(&deviatoricerror);	
-	
+   if(this->amr->thicknesserror_threshold>0)		this->ThicknessZZErrorEstimator(&thicknesserror);
+	if(this->amr->deviatoricerror_threshold>0)	this->ZZErrorEstimator(&deviatoricerror);
+
 	if(my_rank==0){
 		this->amr->ExecuteRefinement(gl_distance,if_distance,deviatoricerror,thicknesserror,
-												&newnumberofvertices,&newnumberofelements,&newx,&newy,&newelementslist); 
+												&newnumberofvertices,&newnumberofelements,&newx,&newy,&newelementslist);
 		newz=xNewZeroInit<IssmDouble>(newnumberofvertices);
 		if(newnumberofvertices<=0 || newnumberofelements<=0) _error_("Error in the ReMeshNeopz.");
@@ -5187,10 +5187,10 @@
 		newelementslist=xNew<int>(newnumberofelements*this->GetElementsWidth());
 	}
-	ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());	
-	ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());	
-	ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());	
-	ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm());	
-
-	/*Assign the pointers*/	
+	ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+	ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+	ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+	ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm());
+
+	/*Assign the pointers*/
 	(*pnewelementslist) 	= newelementslist; //Matlab indexing
 	(*pnewx)				  	= newx;
@@ -5208,5 +5208,5 @@
 /*}}}*/
 void FemModel::InitializeAdaptiveRefinementNeopz(void){/*{{{*/
-	
+
 	/*Define variables*/
 	int my_rank										= IssmComm::GetRank();
@@ -5217,5 +5217,5 @@
 	IssmDouble* z									= NULL;
 	int* elements									= NULL;
-	
+
 	/*Initialize field as NULL for now*/
 	this->amr = NULL;
@@ -5225,10 +5225,10 @@
 	this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elements);
 
-	/*Create initial mesh (coarse mesh) in neopz data structure*/ 
+	/*Create initial mesh (coarse mesh) in neopz data structure*/
 	/*Just CPU #0 should keep AMR object*/
    /*Initialize refinement pattern*/
 	this->SetRefPatterns();
 	this->amr = new AdaptiveMeshRefinement();
-	this->amr->refinement_type=1;//1 is refpattern; 0 is uniform (faster) 
+	this->amr->refinement_type=1;//1 is refpattern; 0 is uniform (faster)
 	/*Get amr parameters*/
 	this->parameters->FindParam(&this->amr->level_max,AmrLevelMaxEnum);
@@ -5243,5 +5243,5 @@
 	this->parameters->FindParam(&this->amr->deviatoricerror_groupthreshold,AmrDeviatoricErrorGroupThresholdEnum);
 	this->parameters->FindParam(&this->amr->deviatoricerror_maximum,AmrDeviatoricErrorMaximumEnum);
-	if(my_rank==0){ 
+	if(my_rank==0){
 		this->amr->CreateInitialMesh(numberofvertices,numberofelements,x,y,elements);
 	}
@@ -5264,5 +5264,5 @@
 	/*Output*/
 	IssmDouble* elementdistance;
-	
+
 	/*Intermediaries*/
    int numberofelements							= this->elements->NumberOfElements();
@@ -5273,6 +5273,6 @@
 	IssmDouble xc,yc,x1,y1,x2,y2,x3,y3;
 	int numberofpoints;
-	
-	/*Get points which level set is zero (center of elements with zero level set, levelset_points is serial)*/	
+
+	/*Get points which level set is zero (center of elements with zero level set, levelset_points is serial)*/
 	this->GetZeroLevelSetPoints(&levelset_points,numberofpoints,levelset_type);
 
@@ -5280,5 +5280,5 @@
       Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
       int sid = element->Sid();
-		element->GetVerticesCoordinates(&xyz_list); 
+		element->GetVerticesCoordinates(&xyz_list);
 		x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1];
 		x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1];
@@ -5290,9 +5290,9 @@
 		for(int j=0;j<numberofpoints;j++){
 			distance =sqrt((xc-levelset_points[2*j])*(xc-levelset_points[2*j])+(yc-levelset_points[2*j+1])*(yc-levelset_points[2*j+1]));
-			mindistance=min(distance,mindistance);		
+			mindistance=min(distance,mindistance);
 		}
 		velementdistance->SetValue(sid,mindistance,INS_VAL);
 		xDelete<IssmDouble>(xyz_list);
-	}	
+	}
 
    /*Assemble*/
Index: /issm/trunk-jpl/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp	(revision 23045)
+++ /issm/trunk-jpl/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp	(revision 23046)
@@ -1,5 +1,5 @@
 /*!\file:  InterpFromGridToMeshx.cpp
  * \brief  "c" core code for interpolating values from a structured grid.
- */ 
+ */
 
 #ifdef HAVE_CONFIG_H
@@ -24,7 +24,4 @@
 	double* y=NULL;
 	int     i;
-
-	// TEST
-	_printf_("\r      N: "<<N<<" M:"<<M<<"\n");
 
 	/*Some checks on arguments: */
@@ -151,5 +148,5 @@
 			 *    |         |     |
 			 * y1 x---------+-----x Q21
-			 *    x1                 x2       
+			 *    x1                 x2
 			 *
 			 */
@@ -225,5 +222,5 @@
 	/*split the rectangle in 2 triangle and
 	 * use Lagrange P1 interpolation
-	 *       
+	 *
 	 *   +3----+2,3' Q12----Q22
 	 *   |    /|     |    /|
@@ -282,5 +279,5 @@
 	_assert_(x<=x2 && x>=x1 && y<=y2 && y>=y1);
 
-	return 
+	return
 	  +Q11*(x2-x)*(y2-y)/((x2-x1)*(y2-y1))
 	  +Q21*(x-x1)*(y2-y)/((x2-x1)*(y2-y1))
@@ -301,5 +298,5 @@
 	 *    |        |         |
 	 * y1 x--------x---------x Q21
-	 *    x1       xm        x2 
+	 *    x1       xm        x2
 	 *
 	 */
Index: /issm/trunk-jpl/src/wrappers/InterpFromGridToMesh/InterpFromGridToMesh.js
===================================================================
--- /issm/trunk-jpl/src/wrappers/InterpFromGridToMesh/InterpFromGridToMesh.js	(revision 23045)
+++ /issm/trunk-jpl/src/wrappers/InterpFromGridToMesh/InterpFromGridToMesh.js	(revision 23046)
@@ -6,5 +6,5 @@
  * Usage:
  *	var data_mesh=InterpFromGridToMesh(xIn,yIn,dataIn,xMeshIn,yMeshIn,defaultValue,interpType);\
- * 
+ *
  *	xIn,yIn						: coordinates of matrix data. (x and y must be in increasing order)
  *	dataIn						: matrix holding the data to be interpolated onto the mesh
@@ -51,5 +51,5 @@
 	var yMesh 			= {};
 	//}}}
-	
+
 
 	/*
@@ -57,5 +57,5 @@
 	*/
 	//{{{
-	
+
 	/*
 		Input
@@ -68,5 +68,5 @@
 	dxHeap.set(new Uint8Array(dx.buffer));
 	x 			= dxHeap.byteOffset;
-	
+
 	dy 			= new Float64Array(yIn);
 	ny 			= dy.length * dy.BYTES_PER_ELEMENT;
@@ -75,5 +75,5 @@
 	dyHeap.set(new Uint8Array(dy.buffer));
 	y 			= dyHeap.byteOffset;
-	
+
 	ddata 		= new Float64Array(dataIn);
 	ndata 		= ddata.length * ddata.BYTES_PER_ELEMENT;
@@ -82,5 +82,5 @@
 	ddataHeap.set(new Uint8Array(ddata.buffer));
 	data 		= ddataHeap.byteOffset;
-	
+
 	dxMesh 		= new Float64Array(xMeshIn);
 	nxMesh 		= dxMesh.length * dxMesh.BYTES_PER_ELEMENT;
@@ -89,5 +89,5 @@
 	dxMeshHeap.set(new Uint8Array(dxMesh.buffer));
 	xMesh 		= dxMeshHeap.byteOffset;
-	
+
 	dyMesh 		= new Float64Array(yMeshIn);
 	nyMesh 		= dyMesh.length * dyMesh.BYTES_PER_ELEMENT;
@@ -96,10 +96,9 @@
 	dyMeshHeap.set(new Uint8Array(dyMesh.buffer));
 	yMesh 		= dyMeshHeap.byteOffset;
-	
+
 	nods 		= xMeshIn.length;
 	meshNumRows	= xMeshIn.length;
-	console.log('InterpFromGridToMesh.js: meshNumRows => ' + meshNumRows);
-	
-	
+
+
 	/*
 		Retrieve interpolation type
@@ -112,5 +111,5 @@
 	}
 	//}}}
-	
+
 	/*
 		Output
@@ -118,8 +117,8 @@
 	pdataMesh = Module._malloc(4);
 	//}}}
-	
+
 	//}}}
-	
-	
+
+
 	/*
 		Declare InterpFromGridToMesh module
@@ -127,6 +126,6 @@
 	//{{{
 	InterpFromGridToMeshModule = Module.cwrap(
-		'InterpFromGridToMeshModule', 
-		'number', 
+		'InterpFromGridToMeshModule',
+		'number',
 		[
 			'number', // output : pdataMesh
@@ -145,6 +144,6 @@
 	);
 	//}}}
-	
-	
+
+
 	/*
 		Call InterpFromGridToMesh module
@@ -152,11 +151,11 @@
 	//{{{
 	InterpFromGridToMeshModule(
-		pdataMesh, 
-		x, 
-		y, 
-		data, 
-		xMesh, 
-		yMesh, 
-		defaultValue, 
+		pdataMesh,
+		x,
+		y,
+		data,
+		xMesh,
+		yMesh,
+		defaultValue,
 		nods,
 		dataNumRowsIn,
@@ -166,6 +165,6 @@
 	);
 	//}}}
-	
-	
+
+
 	/*
 		Dynamic copying from heap
@@ -175,6 +174,6 @@
 	dataMesh	= Module.HEAPF64.slice(dataMeshPtr / 8, dataMeshPtr / 8 + nods);
 	//}}}
-	
-	
+
+
 	/*
 		Free resources
@@ -183,6 +182,6 @@
 	Module._free(pdataMesh);
 	//}}}
-	
-	
+
+
 	return dataMesh;
 }