Index: issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 22902)
+++ issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 22955)
@@ -1454,5 +1454,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;
 
@@ -1478,5 +1478,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;
 
@@ -1502,5 +1502,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;
 
@@ -1526,5 +1526,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;
 
@@ -1550,5 +1550,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;
 
@@ -1574,5 +1574,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;
 
@@ -1619,5 +1619,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;
@@ -1955,5 +1955,5 @@
 								IssmDouble* values    = xNewZeroInit<IssmDouble>(nlines*ncols);
 								IssmDouble* allvalues = xNew<IssmDouble>(nlines*ncols);
-
+								
 								/*Fill-in matrix*/
 								for(int j=0;j<elements->Size();j++){
@@ -1964,5 +1964,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);
@@ -2012,5 +2012,5 @@
 	int domaintype;
 	this->parameters->FindParam(&domaintype,DomainTypeEnum);
-
+	
 	/*1: go throug all elements of this partition and figure out how many
 	 * segments we have (corresopnding to levelset = 0)*/
@@ -2132,11 +2132,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; 
 	}
 
@@ -2269,5 +2269,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;
 		}
@@ -2460,7 +2460,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: */
@@ -2473,5 +2473,5 @@
 	IssmDouble         *surface = NULL;
 	IssmDouble         *bed     = NULL;
-
+			
 	if(VerboseSolution()) _printf0_("   updating vertices positions\n");
 
@@ -2512,5 +2512,5 @@
 
 /*AMR*/
-#if !defined(_HAVE_ADOLC_)
+#if !defined(_HAVE_ADOLC_) 
 void FemModel::ReMesh(void){/*{{{*/
 
@@ -2522,5 +2522,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
@@ -2528,5 +2528,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){
@@ -2541,5 +2541,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);
@@ -2572,5 +2572,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);
@@ -2602,5 +2602,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);
@@ -2619,5 +2619,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
 		}
@@ -2663,5 +2663,5 @@
 	/*Insert bedrock from mismip+ setup*/
 	/*This was used to Misomip project/simulations*/
-
+	
 	if(VerboseSolution())_printf0_("	call Mismip bedrock adjust module\n");
 
@@ -2680,5 +2680,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);
@@ -2693,5 +2693,5 @@
 
 	if(VerboseSolution())_printf0_("	call adjust base and thickness module\n");
-
+	
 	int     numvertices = this->GetElementsWidth();
    IssmDouble rho_water,rho_ice,density,base_float;
@@ -2705,5 +2705,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);
@@ -2717,14 +2717,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;	
 		}
 
@@ -2734,5 +2734,5 @@
 		element->AddInput(BaseEnum,&b[0],P1Enum);
 	}
-
+	
    /*Delete*/
    xDelete<IssmDouble>(phi);
@@ -2762,5 +2762,5 @@
 	Vector<IssmDouble>* input_interpolations	= NULL;
 	IssmDouble* input_interpolations_serial	= NULL;
-   int* pos												= NULL;
+   int* pos												= NULL; 
 	IssmDouble value									= 0;
 
@@ -2782,5 +2782,5 @@
 			P0input_interp = xNew<int>(numP0inputs);
 			P1input_enums  = xNew<int>(numP1inputs);
-			P1input_interp = xNew<int>(numP1inputs);
+			P1input_interp = xNew<int>(numP1inputs);	
 		}
 		numP0inputs = 0;
@@ -2814,6 +2814,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);
@@ -2821,14 +2821,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++){
@@ -2837,6 +2837,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);	
 		}
 	}
@@ -2847,14 +2847,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*/
@@ -2867,5 +2867,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
@@ -2883,5 +2883,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
@@ -2895,5 +2895,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)*/
@@ -2928,15 +2928,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;
@@ -2956,5 +2956,5 @@
 		}
 	}
-
+	
 	/*Cleanup*/
 	xDelete<IssmDouble>(P0inputs);
@@ -2995,5 +2995,5 @@
 
 	if(!this->elements || !this->vertices || !this->results || !this->parameters) return;
-
+	 
 	parameters->FindParam(&step,StepEnum);
 	parameters->FindParam(&time,TimeEnum);
@@ -3013,5 +3013,5 @@
 	this->results->AddResult(new GenericExternalResult<IssmDouble*>(this->results->Size()+1,MeshYEnum,
 					y,numberofvertices,1,step,time));
-
+	
 	/*Cleanup*/
 	xDelete<IssmDouble>(x);
@@ -3074,5 +3074,5 @@
    /*Assemble*/
 	vmasklevelset->Assemble();
-
+	
 	/*Serialize and set output*/
 	(*pmasklevelset)=vmasklevelset->ToMPISerial();
@@ -3088,5 +3088,5 @@
 
 	/*newelementslist is in Matlab indexing*/
-
+	
 	/*Creating connectivity table*/
 	int* connectivity=NULL;
@@ -3099,10 +3099,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;
@@ -3115,6 +3115,6 @@
 			newvertex->connectivity=connectivity[i];
 			newvertex->clone=false;//itapopo check this
-			vertices->AddObject(newvertex);
-		}
+			vertices->AddObject(newvertex);	
+		} 
 	}
 
@@ -3143,5 +3143,5 @@
 			}
 			else newtria->element_type_list=NULL;
-
+			
 			/*Element hook*/
 			int matpar_id=newnumberofelements+1; //retrieve material parameter id (last pointer in femodel->materials)
@@ -3149,5 +3149,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;
@@ -3161,6 +3161,6 @@
 			/*Clean up*/
 			xDelete<int>(vertex_ids);
-			elements->AddObject(newtria);
-		}
+			elements->AddObject(newtria);	
+		} 
 	}
 
@@ -3172,12 +3172,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	    
 }
 /*}}}*/
@@ -3186,8 +3186,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;
@@ -3195,12 +3195,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;
@@ -3216,5 +3216,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]);
@@ -3231,5 +3231,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();
@@ -3237,5 +3237,5 @@
 	IssmDouble* x			= NULL;
 	IssmDouble* y			= NULL;
-	IssmDouble* z			= NULL;
+	IssmDouble* z			= NULL;	
 	int* elementslist 	= NULL;
 	int* elem_vertices	= NULL;
@@ -3246,5 +3246,5 @@
 	/*Get vertices coordinates*/
 	VertexCoordinatesx(&x,&y,&z,femmodel_vertices,false) ;
-
+	
 	/*Get element vertices*/
 	elem_vertices				= xNew<int>(elementswidth);
@@ -3261,5 +3261,5 @@
     	vid3->SetValue(element->sid,elem_vertices[2],INS_VAL);
    }
-
+		
 	/*Assemble*/
    vid1->Assemble();
@@ -3271,5 +3271,5 @@
    id2 = vid2->ToMPISerial();
 	id3 = vid3->ToMPISerial();
-
+	
 	/*Construct elements list*/
 	elementslist=xNew<int>(numberofelements*elementswidth);
@@ -3280,5 +3280,5 @@
 		elementslist[elementswidth*i+2] = reCast<int>(id3[i])+1; //InterpMesh wants Matlab indexing
 	}
-
+	
 	/*Assign pointers*/
 	*px				= x;
@@ -3301,5 +3301,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)
@@ -3308,9 +3308,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
@@ -3318,8 +3318,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();
@@ -3334,5 +3334,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));
@@ -3341,6 +3341,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;
@@ -3348,5 +3348,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*/
@@ -3359,7 +3359,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
@@ -3385,5 +3385,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();
@@ -3393,5 +3393,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");
@@ -3400,12 +3400,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
@@ -3423,5 +3423,5 @@
 								spc,numberofvertices,numberofcols,
 								newx,newy,newnumberofvertices,NULL);
-
+	
 	/*Now, insert the interpolated constraints in the data set (constraints)*/
 	count=0;
@@ -3440,5 +3440,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++;
@@ -3497,11 +3497,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++){
@@ -3523,5 +3523,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);
@@ -3533,12 +3533,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
 			}
@@ -3552,16 +3552,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);
@@ -3573,8 +3573,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++){
@@ -3621,5 +3621,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];
@@ -3646,5 +3646,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*/
 
@@ -3666,5 +3666,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++){
@@ -3674,5 +3674,5 @@
       element->GetInputListOnVertices(tauxy,DeviatoricStressxyEnum);
       element->GetVerticesSidList(elem_vertices);
-
+		
 		/*Integrate*/
 		element->GetVerticesCoordinates(&xyz_list);
@@ -3689,7 +3689,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)
@@ -3705,5 +3705,5 @@
    /*Serialize and set output*/
    (*pelementerror)=velementerror->ToMPISerial();
-
+	
 	/*Cleanup*/
 	xDelete<IssmDouble>(smoothedtauxx);
@@ -3749,5 +3749,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);
@@ -3817,5 +3817,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++){
@@ -3903,5 +3903,5 @@
 	IssmDouble* x					= NULL;
 	IssmDouble* y					= NULL;
-	IssmDouble* z					= NULL;
+	IssmDouble* z					= NULL;	
 	IssmDouble* xyz_list			= NULL;
 	IssmDouble x1,y1,x2,y2,x3,y3;
@@ -3912,5 +3912,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];
@@ -3945,9 +3945,9 @@
 		_error_("level set type not implemented yet!");
 	}
-
+	
 	/*Outputs*/
 	IssmDouble* zerolevelset_points			= NULL;
 	int npoints										= 0;
-
+	
 	/*Intermediaries*/
  	int elementswidth                   	= this->GetElementsWidth();
@@ -3962,5 +3962,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++){
@@ -3969,13 +3969,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.;
@@ -4007,5 +4007,5 @@
 		}
 	}
-
+	
 	/*Assign outputs*/
 	numberofpoints				= npoints;
@@ -4047,5 +4047,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: */
 
@@ -4140,8 +4140,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());
@@ -4162,5 +4162,5 @@
 		}
 	}
-
+	
 	/*One last time: */
 	pUp->Assemble();
@@ -4181,9 +4181,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));
@@ -4209,5 +4209,5 @@
 		}
 	}
-
+	
 	/*One last time: */
 	pUp->Assemble();
@@ -4226,5 +4226,5 @@
 #endif
 #ifdef _HAVE_SEALEVELRISE_
-void FemModel::SealevelriseEustatic(Vector<IssmDouble>* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius) { /*{{{*/
+void FemModel::SealevelriseEustatic(Vector<IssmDouble>* pRSLgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop) { /*{{{*/
 
 	/*serialized vectors:*/
@@ -4262,16 +4262,16 @@
 		if(i<ns){
 
-			if(VerboseConvergence())if(i%100==0)_printf0_("\r" << "      convolution progress: " << (double)i/(double)ns*100 << "%  ");
+			if(VerboseConvergence())if(i%100==0)_printf0_("\r" << "              convolution progress: " << (double)i/(double)ns*100 << "%  ");
 
 			Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-			element->SealevelriseEustatic(pSgi,&eustatic_cpu_e,latitude,longitude,radius,oceanarea,eartharea);
+			element->SealevelriseEustatic(pRSLgi,&eustatic_cpu_e,latitude,longitude,radius,oceanarea,eartharea);
 			eustatic_cpu+=eustatic_cpu_e;
 		}
-		if(i%100==0)pSgi->Assemble();
+		if(i%loop==0)pRSLgi->Assemble();
 	}
 	if(VerboseConvergence())_printf0_("\n");
-
+		
 	/*One last time: */
-	pSgi->Assemble();
+	pRSLgi->Assemble();
 
 	/*Sum all eustatic components from all cpus:*/
@@ -4285,16 +4285,16 @@
 }
 /*}}}*/
-void FemModel::SealevelriseNonEustatic(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution){/*{{{*/
+void FemModel::SealevelriseNonEustatic(Vector<IssmDouble>* pRSLgo, Vector<IssmDouble>* pRSLg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution,int loop){/*{{{*/
 
 	/*serialized vectors:*/
-	IssmDouble* Sg_old=NULL;
-
+	IssmDouble* RSLg_old=NULL;
+	
 	IssmDouble  eartharea=0;
 	IssmDouble  eartharea_cpu=0;
 
 	int         ns,nsmax;
-
+	
 	/*Serialize vectors from previous iteration:*/
-	Sg_old=pSg_old->ToMPISerial();
+	RSLg_old=pRSLg_old->ToMPISerial();
 
 	/*Go through elements, and add contribution from each element to the deflection vector wg:*/
@@ -4306,5 +4306,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());
@@ -4317,20 +4317,20 @@
 	for(int i=0;i<nsmax;i++){
 		if(i<ns){
-			if(verboseconvolution)if(VerboseConvergence())if(i%100==0)_printf_("\r" << "      convolution progress: " << (double)i/(double)ns*100 << "%   ");
+			if(verboseconvolution)if(VerboseConvergence())if(i%100==0)_printf0_("\r" << "              convolution progress: " << (double)i/(double)ns*100 << "%   ");
 			Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-			element->SealevelriseNonEustatic(pSgo,Sg_old,latitude,longitude,radius,eartharea);
-		}
-		if(i%100==0)pSgo->Assemble();
-	}
-	if(verboseconvolution)if(VerboseConvergence())_printf_("\n");
-
+			element->SealevelriseNonEustatic(pRSLgo,RSLg_old,latitude,longitude,radius,eartharea);
+		}
+		if(i%loop==0)pRSLgo->Assemble();
+	}
+	if(verboseconvolution)if(VerboseConvergence())_printf0_("\n");
+	
 	/*Free ressources:*/
-	xDelete<IssmDouble>(Sg_old);
-}
-/*}}}*/
-void FemModel::SealevelriseRotationalFeedback(Vector<IssmDouble>* pSgo_rot, Vector<IssmDouble>* pSg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/
+	xDelete<IssmDouble>(RSLg_old);
+}
+/*}}}*/
+void FemModel::SealevelriseRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/
 
 	/*serialized vectors:*/
-	IssmDouble* Sg_old=NULL;
+	IssmDouble* RSLg_old=NULL;
 	IssmDouble  eartharea=0;
 	IssmDouble  eartharea_cpu=0;
@@ -4341,5 +4341,5 @@
 
 	/*Serialize vectors from previous iteration:*/
-	Sg_old=pSg_old->ToMPISerial();
+	RSLg_old=pRSLg_old->ToMPISerial();
 
 	/*First, figure out the area of the ocean, which is needed to compute the eustatic component: */
@@ -4355,5 +4355,5 @@
 	for(int i=0;i<elements->Size();i++){
 		Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-		element->SealevelriseMomentOfInertia(&moi_list[0],Sg_old,eartharea);
+		element->SealevelriseMomentOfInertia(&moi_list[0],RSLg_old,eartharea);
 		moi_list_cpu[0] += moi_list[0];
 		moi_list_cpu[1] += moi_list[1];
@@ -4399,36 +4399,37 @@
 						(-m3/6.0 + 0.5*m3*cos(2.0*lati) - 0.5*sin(2.*lati)*(m1*cos(longi)+m2*sin(longi)));
 
-		pSgo_rot->SetValue(sid,value,INS_VAL); //INS_VAL ensures that you don't add several times
+		pRSLgo_rot->SetValue(sid,value,INS_VAL); //INS_VAL ensures that you don't add several times
 	}
 
 	/*Assemble mesh velocity*/
-	pSgo_rot->Assemble();
+	pRSLgo_rot->Assemble();
 
 	/*Assign output pointers:*/
-	*pIxz=moi_list[0];
-	*pIyz=moi_list[1];
-	*pIzz=moi_list[2];
+	if(pIxz)*pIxz=moi_list[0];
+	if(pIyz)*pIyz=moi_list[1];
+	if(pIzz)*pIzz=moi_list[2];
 
 	/*Free ressources:*/
-	xDelete<IssmDouble>(Sg_old);
-
-}
-/*}}}*/
-void FemModel::SealevelriseGeodetic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){/*{{{*/
+	xDelete<IssmDouble>(RSLg_old);
+	
+	
+}
+/*}}}*/
+void FemModel::SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pRSLg, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz){/*{{{*/
 
 	/*serialized vectors:*/
-	IssmDouble* Sg=NULL;
-
+	IssmDouble* RSLg=NULL;
+	
 	IssmDouble  eartharea=0;
 	IssmDouble  eartharea_cpu=0;
 
 	int         ns,nsmax;
-
+	
 	/*Serialize vectors from previous iteration:*/
-	Sg=pSg->ToMPISerial();
+	RSLg=pRSLg->ToMPISerial();
 
 	/*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++){
@@ -4446,37 +4447,37 @@
 	for(int i=0;i<nsmax;i++){
 		if(i<ns){
+			if(VerboseConvergence())if(i%100==0)_printf0_("\r" << "              convolution progress: " << (double)i/(double)ns*100 << "%  ");
 			Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-			element->SealevelriseGeodetic(pUp,pNorth,pEast,Sg,latitude,longitude,radius,xx,yy,zz,eartharea);
-		}
-		if(i%100==0){
+			element->SealevelriseGeodetic(pUp,pNorth,pEast,RSLg,latitude,longitude,radius,xx,yy,zz,eartharea,horiz);
+		}
+		if(i%loop==0){
 			pUp->Assemble();
-			pNorth->Assemble();
-			pEast->Assemble();
-		}
-	}
-
+			if (horiz){
+				pNorth->Assemble();
+				pEast->Assemble();
+			}
+		}
+	}
+	
 	/*One last time: */
 	pUp->Assemble();
-	pNorth->Assemble();
-	pEast->Assemble();
+	if(horiz){
+		pNorth->Assemble();
+		pEast->Assemble();
+	}
+	if(VerboseConvergence())_printf0_("\n");
 
 	/*Free ressources:*/
-	xDelete<IssmDouble>(Sg);
-	xDelete<IssmDouble>(latitude);
-	xDelete<IssmDouble>(longitude);
-	xDelete<IssmDouble>(radius);
-	xDelete<IssmDouble>(xx);
-	xDelete<IssmDouble>(yy);
-	xDelete<IssmDouble>(zz);
-}
-/*}}}*/
-IssmDouble FemModel::SealevelriseOceanAverage(Vector<IssmDouble>* Sg) { /*{{{*/
-
-	IssmDouble* Sg_serial=NULL;
+	xDelete<IssmDouble>(RSLg);
+}
+/*}}}*/
+IssmDouble FemModel::SealevelriseOceanAverage(Vector<IssmDouble>* RSLg) { /*{{{*/
+
+	IssmDouble* RSLg_serial=NULL;
 	IssmDouble  oceanvalue,oceanvalue_cpu;
 	IssmDouble  oceanarea,oceanarea_cpu;
 
 	/*Serialize vectors from previous iteration:*/
-	Sg_serial=Sg->ToMPISerial();
+	RSLg_serial=RSLg->ToMPISerial();
 
 	/*Initialize:*/
@@ -4488,15 +4489,15 @@
 		Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
 		oceanarea_cpu += element->OceanArea();
-		oceanvalue_cpu += element->OceanAverage(Sg_serial);
+		oceanvalue_cpu += element->OceanAverage(RSLg_serial);
 	}
 	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());
 
 	/*Free ressources:*/
-	xDelete<IssmDouble>(Sg_serial);
-
+	xDelete<IssmDouble>(RSLg_serial);
+	
 	return oceanvalue/oceanarea;
 }
@@ -4514,5 +4515,5 @@
 	int*                eplzigzag_counter =	NULL;
 	int                 eplflip_lock;
-
+	
 	HydrologyDCEfficientAnalysis* effanalysis =  new HydrologyDCEfficientAnalysis();
 	HydrologyDCInefficientAnalysis* inefanalysis =  new HydrologyDCInefficientAnalysis();
@@ -4521,8 +4522,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));
@@ -4542,6 +4543,6 @@
 	/*Assemble and serialize*/
 	mask->Assemble();
-	serial_mask=mask->ToMPISerial();
-
+	serial_mask=mask->ToMPISerial();	
+	
 	xDelete<int>(eplzigzag_counter);
 	xDelete<IssmDouble>(serial_rec);
@@ -4585,83 +4586,8 @@
 	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;
 	if(VerboseSolution()) _printf0_("   Number of active nodes in EPL layer: "<< counter <<"\n");
-
-	/*Update dof indexings*/
-	this->UpdateConstraintsx();
-
-}
-/*}}}*/
-void FemModel::HydrologyIDSupdateDomainx(IssmDouble* pIDScount){ /*{{{*/
-
-	bool                isthermal;
-	Vector<IssmDouble>* mask				= NULL;
-	Vector<IssmDouble>* active				= NULL;
-	IssmDouble*         serial_mask	= NULL;
-	IssmDouble*         serial_active	= NULL;
-
-	HydrologyDCInefficientAnalysis* inefanalysis =  new HydrologyDCInefficientAnalysis();
-	parameters->FindParam(&isthermal,TransientIsthermalEnum);
-
-	/*When solving a thermal model we update the thawed nodes*/
-	if(isthermal){
-		/*Step 1: update mask, the mask correspond to thawed nodes (that have a meltingrate)*/
-		mask=new Vector<IssmDouble>(this->nodes->NumberOfNodes(HydrologyDCInefficientAnalysisEnum));
-
-		for (int i=0;i<elements->Size();i++){
-			Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-			inefanalysis->HydrologyIDSGetMask(mask,element);
-		}
-		/*Assemble and serialize*/
-		mask->Assemble();
-		serial_mask=mask->ToMPISerial();
-		delete mask;
-	}
-	/*for other cases we just grab the mask from the initialisation value*/
-	else{
-		GetVectorFromInputsx(&serial_mask,this,HydrologydcMaskThawedNodeEnum,NodeSIdEnum);
-	}
-	/*Update Mask and elementize*/
-	InputUpdateFromVectorx(this,serial_mask,HydrologydcMaskThawedNodeEnum,NodeSIdEnum);
-	xDelete<IssmDouble>(serial_mask);
-	inefanalysis->ElementizeIdsMask(this);
-
-	/*get node mask coherent with element mask*/
-	active=new Vector<IssmDouble>(nodes->NumberOfNodes(HydrologyDCInefficientAnalysisEnum));
-	for (int i=0;i<elements->Size();i++){
-		Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-		inefanalysis->HydrologyIdsGetActive(active,element);
-	}
-
-	/*Assemble and serialize*/
-	active->Assemble();
-	serial_active=active->ToMPISerial();
-	delete active;
-
-	/*Update node activation accordingly*/
-	int counter =0;
-	for (int i=0;i<nodes->Size();i++){
-		Node* node=xDynamicCast<Node*>(nodes->GetObjectByOffset(i));
-		if(node->InAnalysis(HydrologyDCInefficientAnalysisEnum)){
-			if(serial_active[node->Sid()]==1.){
-				node->Activate();
-				if(!node->IsClone()) counter++;
-			}
-			else{
-				node->Deactivate();
-			}
-		}
-	}
-
-	xDelete<IssmDouble>(serial_active);
-	delete inefanalysis;
-	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());
-	counter=sum_counter;
-	*pIDScount = counter;
-	if(VerboseSolution()) _printf0_("   Number of active nodes in IDS layer: "<< counter <<"\n");
 
 	/*Update dof indexings*/
@@ -4706,5 +4632,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;
@@ -4791,5 +4717,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: */
@@ -4806,10 +4732,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));
@@ -4826,5 +4752,5 @@
 	size_t* poutputbuffersize;
 
-
+	
 	/*Before we delete the profiler, report statistics for this run: */
 	profiler->Stop(TOTAL);  //final tagging
@@ -4839,5 +4765,5 @@
 				);
 	_printf0_("\n");
-
+	
 	/*Before we close the output file, recover the buffer and size:*/
 	outputbufferparam = xDynamicCast<GenericParam<char**>*>(this->parameters->FindParamObject(OutputBufferPointerEnum));
@@ -4879,5 +4805,5 @@
 
 	/*Open output file once for all and add output file descriptor to parameters*/
-	output_fid=open_memstream(&outputbuffer,&outputsize);
+	output_fid=open_memstream(&outputbuffer,&outputsize); 
 	if(output_fid==NULL)_error_("could not initialize output stream");
 	this->parameters->SetParam(output_fid,OutputFilePointerEnum);
@@ -4887,4 +4813,5 @@
 }/*}}}*/
 #endif
+
 
 #if defined(_HAVE_BAMG_) && !defined(_HAVE_ADOLC_)
@@ -4920,5 +4847,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);
@@ -4927,5 +4854,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);
@@ -4947,8 +4874,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*/
@@ -4983,5 +4910,5 @@
 	/*Create bamg data structures for bamg*/
 	this->amrbamg = new AmrBamg();
-
+	
 	/*Get amr parameters*/
 	this->parameters->FindParam(&hmin,AmrHminEnum);
@@ -5007,5 +4934,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);
 	}
@@ -5021,5 +4948,5 @@
 
 	if(!hmaxvertices) _error_("hmaxvertices is NULL!\n");
-
+	
 	/*Intermediaries*/
 	int numberofvertices			 = this->vertices->NumberOfVertices();
@@ -5028,5 +4955,5 @@
 
 	switch(levelset_type){
-		case MaskGroundediceLevelsetEnum:
+		case MaskGroundediceLevelsetEnum: 
 			threshold	= this->amrbamg->groundingline_distance;
 			resolution	= this->amrbamg->groundingline_resolution;
@@ -5042,5 +4969,5 @@
 	this->GetVerticeDistanceToZeroLevelSet(&verticedistance,levelset_type);
 	if(!verticedistance) _error_("verticedistance is NULL!\n");
-
+	
 	/*Fill hmaxVertices*/
 	for(int i=0;i<numberofvertices;i++){
@@ -5056,5 +4983,5 @@
 /*}}}*/
 void FemModel::GethmaxVerticesFromEstimators(IssmDouble* hmaxvertices,int errorestimator_type){/*{{{*/
-
+   
 	if(!hmaxvertices) _error_("hmaxvertices is NULL!\n");
 
@@ -5064,5 +4991,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;
@@ -5077,5 +5004,5 @@
 	/*Fill variables*/
 	switch(errorestimator_type){
-		case ThicknessErrorEstimatorEnum:
+		case ThicknessErrorEstimatorEnum: 
 			threshold		= this->amrbamg->thicknesserror_threshold;
 			groupthreshold	= this->amrbamg->thicknesserror_groupthreshold;
@@ -5102,10 +5029,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++){
@@ -5121,5 +5048,5 @@
 		error_vertices[v2]+=error_elements[i];
 		error_vertices[v3]+=error_elements[i];
-	}
+	}	
 
 	/*Fill hmaxvertices with the criteria*/
@@ -5135,5 +5062,5 @@
 			}
 		}
-		/*Now, fill the hmaxvertices if requested*/
+		/*Now, fill the hmaxvertices if requested*/	  
 		if(refine){
 			for(int j=0;j<elementswidth;j++){
@@ -5165,5 +5092,5 @@
 	/*Output*/
 	IssmDouble* verticedistance;
-
+	
 	/*Intermediaries*/
    int numberofvertices       = this->vertices->NumberOfVertices();
@@ -5177,6 +5104,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);
 
@@ -5187,7 +5114,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*/
@@ -5223,10 +5150,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.");
@@ -5242,10 +5169,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;
@@ -5263,5 +5190,5 @@
 /*}}}*/
 void FemModel::InitializeAdaptiveRefinementNeopz(void){/*{{{*/
-
+	
 	/*Define variables*/
 	int my_rank										= IssmComm::GetRank();
@@ -5272,5 +5199,5 @@
 	IssmDouble* z									= NULL;
 	int* elements									= NULL;
-
+	
 	/*Initialize field as NULL for now*/
 	this->amr = NULL;
@@ -5280,10 +5207,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);
@@ -5298,5 +5225,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);
 	}
@@ -5319,5 +5246,5 @@
 	/*Output*/
 	IssmDouble* elementdistance;
-
+	
 	/*Intermediaries*/
    int numberofelements							= this->elements->NumberOfElements();
@@ -5328,6 +5255,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);
 
@@ -5335,5 +5262,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];
@@ -5345,9 +5272,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*/
