Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 21381)
+++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 21382)
@@ -108,5 +108,5 @@
 
 	bool dakota_analysis,ismovingfront,isenthalpy;
-	int frictionlaw,basalforcing_model;
+	int frictionlaw,basalforcing_model,materialstype;
 	int FrictionCoupling;
 	
@@ -134,4 +134,5 @@
 	iomodel->FindConstant(&ismovingfront,"md.transient.ismovingfront");
 	iomodel->FindConstant(&frictionlaw,"md.friction.law");
+	iomodel->FindConstant(&materialstype,"md.materials.type");
 
 	iomodel->FetchDataToInput(elements,"md.geometry.thickness",ThicknessEnum);
@@ -145,6 +146,4 @@
 		iomodel->FetchDataToInput(elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
 	}
-	iomodel->FetchDataToInput(elements,"md.materials.rheology_B",MaterialsRheologyBEnum);
-	iomodel->FetchDataToInput(elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
 	iomodel->FetchDataToInput(elements,"md.initialization.pressure",PressureEnum);
 	iomodel->FetchDataToInput(elements,"md.initialization.temperature",TemperatureEnum);
@@ -171,4 +170,21 @@
 			iomodel->FetchDataToInput(elements,"md.basalforcings.geothermalflux",BasalforcingsGeothermalfluxEnum);
 			break;
+	}
+
+	/*Rheology type*/
+	iomodel->FetchDataToInput(elements,"md.materials.rheology_B",MaterialsRheologyBEnum);
+	switch(materialstype){
+		case MatdamageiceEnum:
+			iomodel->FetchDataToInput(elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
+			break;
+		case MatestarEnum:
+			iomodel->FetchDataToInput(elements,"md.materials.rheology_Ec",MaterialsRheologyEcEnum);
+			iomodel->FetchDataToInput(elements,"md.materials.rheology_Es",MaterialsRheologyEsEnum);
+			break;
+		case MaticeEnum:
+			iomodel->FetchDataToInput(elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
+			break;
+		default:
+			_error_("not supported");
 	}
 	
@@ -589,5 +605,5 @@
 		}
 
-		/*Artifficial diffusivity*/
+		/*Artificial diffusivity*/
 		if(stabilization==1){
 			element->ElementSizes(&hx,&hy,&hz);
@@ -1452,4 +1468,5 @@
 	int         i,rheology_law;
 	IssmDouble  B_average,s_average,T_average=0.,P_average=0.;
+	IssmDouble  n=3.0;
 	int        *doflist   = NULL;
 	IssmDouble *xyz_list  = NULL;
@@ -1477,4 +1494,5 @@
 
 	/*Get all inputs and parameters*/
+	if(element->material->ObjectEnum()!=MatestarEnum) n=element->GetMaterialParameter(MaterialsRheologyNEnum);
 	element->GetInputValue(&converged,ConvergedEnum);
 	element->GetInputListOnNodes(&pressure[0],PressureEnum);
@@ -1506,5 +1524,5 @@
 				break;
 			case CuffeyTemperateEnum:
-				for(i=0;i<numnodes;i++) B[i]=CuffeyTemperate(temperature[i], waterfraction[i], element->GetMaterialParameter(MaterialsRheologyNEnum));
+				for(i=0;i<numnodes;i++) B[i]=CuffeyTemperate(temperature[i], waterfraction[i],n);
 				element->AddInput(MaterialsRheologyBEnum,&B[0],element->GetElementType());
 				break;
@@ -1515,9 +1533,9 @@
 			case ArrheniusEnum:
 				element->GetVerticesCoordinates(&xyz_list);
-				for(i=0;i<numnodes;i++) B[i]=Arrhenius(temperature[i],surface[i]-xyz_list[i*3+2],element->GetMaterialParameter(MaterialsRheologyNEnum));
+				for(i=0;i<numnodes;i++) B[i]=Arrhenius(temperature[i],surface[i]-xyz_list[i*3+2],n);
 				element->AddInput(MaterialsRheologyBEnum,&B[0],element->GetElementType());
 				break;
 			case LliboutryDuvalEnum:
-				for(i=0;i<numnodes;i++) B[i]=LliboutryDuval(values[i],pressure[i],element->GetMaterialParameter(MaterialsRheologyNEnum),element->GetMaterialParameter(MaterialsBetaEnum),element->GetMaterialParameter(ConstantsReferencetemperatureEnum),element->GetMaterialParameter(MaterialsHeatcapacityEnum),element->GetMaterialParameter(MaterialsLatentheatEnum)); 
+				for(i=0;i<numnodes;i++) B[i]=LliboutryDuval(values[i],pressure[i],n,element->GetMaterialParameter(MaterialsBetaEnum),element->GetMaterialParameter(ConstantsReferencetemperatureEnum),element->GetMaterialParameter(MaterialsHeatcapacityEnum),element->GetMaterialParameter(MaterialsLatentheatEnum)); 
 				element->AddInput(MaterialsRheologyBEnum,&B[0],element->GetElementType()); 
 				break; 
Index: /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp	(revision 21381)
+++ /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp	(revision 21382)
@@ -73,5 +73,5 @@
 void ThermalAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
 
-	int frictionlaw,basalforcing_model;
+	int frictionlaw,basalforcing_model,materialstype;
 	int FrictionCoupling;
 	/*Now, is the model 3d? otherwise, do nothing: */
@@ -93,4 +93,5 @@
 	iomodel->FindConstant(&ismovingfront,"md.transient.ismovingfront");
 	iomodel->FindConstant(&frictionlaw,"md.friction.law");
+	iomodel->FindConstant(&materialstype,"md.materials.type");
 
 	iomodel->FetchDataToInput(elements,"md.geometry.thickness",ThicknessEnum);
@@ -106,6 +107,4 @@
 	iomodel->FetchDataToInput(elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
 	iomodel->FetchDataToInput(elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
-	iomodel->FetchDataToInput(elements,"md.materials.rheology_B",MaterialsRheologyBEnum);
-	iomodel->FetchDataToInput(elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
 	iomodel->FetchDataToInput(elements,"md.initialization.pressure",PressureEnum);
 	iomodel->FetchDataToInput(elements,"md.initialization.temperature",TemperatureEnum);
@@ -116,4 +115,21 @@
 	InputUpdateFromConstantx(elements,0.,VyMeshEnum);
 	InputUpdateFromConstantx(elements,0.,VzMeshEnum);
+
+	/*Rheology type*/
+	iomodel->FetchDataToInput(elements,"md.materials.rheology_B",MaterialsRheologyBEnum);
+	switch(materialstype){
+		case MatdamageiceEnum:
+			iomodel->FetchDataToInput(elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
+			break;
+		case MatestarEnum:
+			iomodel->FetchDataToInput(elements,"md.materials.rheology_Ec",MaterialsRheologyEcEnum);
+			iomodel->FetchDataToInput(elements,"md.materials.rheology_Es",MaterialsRheologyEsEnum);
+			break;
+		case MaticeEnum:
+			iomodel->FetchDataToInput(elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
+			break;
+		default:
+			_error_("not supported");
+	}
 	if(ismovingfront){
 		iomodel->FetchDataToInput(elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum); // required for updating active nodes
@@ -722,4 +738,5 @@
 	int        *doflist   = NULL;
 	IssmDouble *xyz_list  = NULL;
+	IssmDouble  n=3.0;
 	bool        hack      = false;
 
@@ -756,4 +773,5 @@
 
 	/*Get all inputs and parameters*/
+	if(element->material->ObjectEnum()!=MatestarEnum) n=element->GetMaterialParameter(MaterialsRheologyNEnum);
 	element->GetInputValue(&converged,ConvergedEnum);
 	if(converged){
@@ -782,5 +800,5 @@
 			case ArrheniusEnum:{
 				element->GetVerticesCoordinates(&xyz_list);
-				for(i=0;i<numnodes;i++) B[i]=Arrhenius(values[i],surface[i]-xyz_list[i*3+2],element->GetMaterialParameter(MaterialsRheologyNEnum));
+				for(i=0;i<numnodes;i++) B[i]=Arrhenius(values[i],surface[i]-xyz_list[i*3+2],n);
 				element->AddInput(MaterialsRheologyBEnum,&B[0],element->GetElementType());
 				break;