Index: ../trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp	(revision 18503)
+++ ../trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp	(revision 18504)
@@ -9,6 +9,9 @@
 	return 2;
 }/*}}}*/
 void SeaiceAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
+	parameters->AddObject(iomodel->CopyConstantObject(SeaiceMinConcentrationEnum));
+	parameters->AddObject(iomodel->CopyConstantObject(SeaiceMinThicknessEnum));
+	parameters->AddObject(iomodel->CopyConstantObject(SeaiceMaxThicknessEnum));
 }/*}}}*/
 void SeaiceAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
 
@@ -33,12 +36,17 @@
 	iomodel->FetchDataToInput(elements,BasalforcingsOceanVyEnum);
 	iomodel->FetchDataToInput(elements,SurfaceforcingsWindVxEnum);
 	iomodel->FetchDataToInput(elements,SurfaceforcingsWindVyEnum);
-	iomodel->FetchDataToInput(elements,MaterialsDamageEnum);
+	iomodel->FetchDataToInput(elements,DamageEnum);
 	iomodel->FetchDataToInput(elements,VxStarEnum);
 	iomodel->FetchDataToInput(elements,VyStarEnum);
 	iomodel->FetchDataToInput(elements,VxEnum);
 	iomodel->FetchDataToInput(elements,VyEnum);
 	iomodel->FetchDataToInput(elements,SeaiceCoriolisFactorEnum);
+	iomodel->FetchDataToInput(elements,StressTensorPredictorxxEnum);
+	iomodel->FetchDataToInput(elements,StressTensorPredictoryyEnum);
+	iomodel->FetchDataToInput(elements,StressTensorPredictorxyEnum);
+	iomodel->FetchDataToInput(elements,MeshXEnum);
+	iomodel->FetchDataToInput(elements,MeshYEnum);
 }/*}}}*/
 void SeaiceAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
 
@@ -70,7 +78,8 @@
 
 	/* Check if there is ice in this element */
 	Input* concentration_input = element->GetInput(SeaiceConcentrationEnum); _assert_(concentration_input);
-	if(concentration_input->Max()==0.) return NULL;
+	IssmDouble c_min; element->FindParam(&c_min,SeaiceMinConcentrationEnum);
+	if(concentration_input->Max()<=c_min) return NULL;
 
 	/*Intermediaries */
 	IssmDouble  Jdet,D_scalar,dt,thickness,concentration;
@@ -95,13 +104,18 @@
 	element->FindParam(&ocean_lin_drag_coef,BasalforcingsOceanLinDragCoefEnum);
 	element->FindParam(&ocean_quad_drag_coef,BasalforcingsOceanQuadDragCoefEnum);
 	element->FindParam(&ocean_turning_angle,BasalforcingsOceanTurningAngleEnum);
-	IssmDouble rho_i           = element->GetMaterialParameter(MaterialsRhoIceEnum);
+	IssmDouble rho_i                  = element->GetMaterialParameter(MaterialsRhoIceEnum);
+	IssmDouble time_relaxation_stress = element->GetMaterialParameter(MaterialsTimeRelaxationStressEnum);
 	Input*     thickness_input = element->GetInput(SeaiceThicknessEnum);       _assert_(thickness_input);
 	Input*     vx_input        = element->GetInput(VxEnum);                    _assert_(vx_input);
 	Input*     vy_input        = element->GetInput(VyEnum);                    _assert_(vy_input);
 	Input*     oceanvx_input   = element->GetInput(BasalforcingsOceanVxEnum);  _assert_(oceanvx_input);
 	Input*     oceanvy_input   = element->GetInput(BasalforcingsOceanVyEnum);  _assert_(oceanvy_input);
 
+	/*Get minimum inertia to avoid 0 in the line, and time_ratio*/
+	IssmDouble minimum_inertia = rho_i*0.01;
+	IssmDouble time_ratio=(1-dt/time_relaxation_stress);
+
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=element->NewGauss(2);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
@@ -115,18 +129,18 @@
 
 		/*Create first part of the stiffness matrix*/
 		this->GetM(M,element,xyz_list,gauss);
-		D_scalar=rho_i*thickness*gauss->weight*Jdet;
+		D_scalar=(rho_i*thickness+minimum_inertia)/dt*gauss->weight*Jdet;
 		TripleMultiply(M,1,numdof,1,
 					&D_scalar,1,1,0,
 					M,1,numdof,0,
 					&Ke->values[0],1);
 
-		/*Create Elasitc part*/
+		/*Create Elastic part*/
 		this->GetB(B,element,xyz_list,gauss);
 		this->CreateCTensor(&C[0][0],element,gauss);
 		for(int i=0;i<3;i++){
 			for(int j=0;j<3;j++){
-				C[i][j] = gauss->weight*Jdet*dt*C[i][j];
+				C[i][j] = dt*time_ratio*gauss->weight*Jdet*C[i][j];
 			}
 		}
 		TripleMultiply(B,3,numdof,1,
@@ -140,7 +154,7 @@
 		vx_input->GetInputValue(&vx,gauss);
 		vy_input->GetInputValue(&vy,gauss);
 		vnorm = sqrt(pow(oceanvx-vx,2) + pow(oceanvy-vy,2));
-		D_scalar = dt*concentration*ocean_coef*rho_ocean*(ocean_lin_drag_coef+ocean_quad_drag_coef*vnorm)*cos(ocean_turning_angle)*gauss->weight*Jdet;
+		D_scalar = concentration*ocean_coef*rho_ocean*(ocean_lin_drag_coef+ocean_quad_drag_coef*vnorm)*cos(ocean_turning_angle)*gauss->weight*Jdet;
 		TripleMultiply(M,1,numdof,1,
 					&D_scalar,1,1,0,
 					M,1,numdof,0,
@@ -159,15 +173,17 @@
 
 	/* Check if there is ice in this element */
 	Input* concentration_input = element->GetInput(SeaiceConcentrationEnum); _assert_(concentration_input);
-	if(concentration_input->Max()==0.) return NULL;
+	IssmDouble c_min; element->FindParam(&c_min,SeaiceMinConcentrationEnum);
+	if(concentration_input->Max()<=c_min) return NULL;
 
 	/*Intermediaries */
-	IssmDouble  air_coef,ocean_coef,constant_part;
-	IssmDouble  rho_ice,rho_air,rho_ocean,gravity;
+	IssmDouble  air_coef,ocean_coef,constant_part,time_relaxation_stress;
+	IssmDouble  rho_ice,rho_air,rho_ocean,gravity,omega;
 	IssmDouble  vx,vy,vxstar,vystar,windvx,windvy,oceanvx,oceanvy,vnorm;
 	IssmDouble  air_lin_drag_coef,air_quad_drag_coef;
 	IssmDouble  ocean_lin_drag_coef,ocean_quad_drag_coef;
-	IssmDouble  concentration,thickness,coriolis_factor,dt,Jdet;
+	IssmDouble  concentration,thickness,coriolis_factor,dt,Jdet,D_scalar;
+	IssmDouble  sigma_xx,sigma_yy,sigma_xy,sigma_vec[3];
 	IssmDouble  ocean_turning_angle,dssh[2];
 	IssmDouble* xyz_list = NULL;
 
@@ -176,8 +192,10 @@
 	int numdof   = 2*numnodes;
 
 	/*Initialize Element vector and vectors*/
-	ElementVector* pe    = element->NewElementVector();
-	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+	ElementVector* pe      = element->NewElementVector();
+	IssmDouble*    B       = xNew<IssmDouble>(3*numdof);
+	IssmDouble*    basis   = xNew<IssmDouble>(numnodes);
+	IssmDouble     D[3][3] = {0.};
 
 	/*Retrieve all inputs and parameters*/
 	element->GetVerticesCoordinates(&xyz_list);
@@ -191,8 +209,10 @@
 	element->FindParam(&ocean_lin_drag_coef,BasalforcingsOceanLinDragCoefEnum);
 	element->FindParam(&ocean_quad_drag_coef,BasalforcingsOceanQuadDragCoefEnum);
 	element->FindParam(&ocean_turning_angle,BasalforcingsOceanTurningAngleEnum);
-	rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
-	gravity = element->GetMaterialParameter(ConstantsGEnum);
+	rho_ice                = element->GetMaterialParameter(MaterialsRhoIceEnum);
+	gravity                = element->GetMaterialParameter(ConstantsGEnum);
+	omega                  = element->GetMaterialParameter(ConstantsOmegaEnum);
+	time_relaxation_stress = element->GetMaterialParameter(MaterialsTimeRelaxationStressEnum);
 	Input* thickness_input     = element->GetInput(SeaiceThicknessEnum);       _assert_(thickness_input);
 	Input* coriolis_fact_input = element->GetInput(SeaiceCoriolisFactorEnum);  _assert_(coriolis_fact_input);
 	Input* vx_input            = element->GetInput(VxEnum);                    _assert_(vx_input);
@@ -204,7 +224,13 @@
 	Input* oceanvx_input       = element->GetInput(BasalforcingsOceanVxEnum);  _assert_(oceanvx_input);
 	Input* oceanvy_input       = element->GetInput(BasalforcingsOceanVyEnum);  _assert_(oceanvy_input);
 	Input* oceanssh_input      = element->GetInput(BasalforcingsOceanSshEnum); _assert_(oceanssh_input);
+	Input* sigma_xx_input      = element->GetInput(StressTensorxxEnum);        _assert_(sigma_xx_input);
+	Input* sigma_yy_input      = element->GetInput(StressTensoryyEnum);        _assert_(sigma_yy_input);
+	Input* sigma_xy_input      = element->GetInput(StressTensorxyEnum);        _assert_(sigma_xy_input);
 
+	/*Calculate time_ratio*/
+	IssmDouble time_ratio=(1-dt/time_relaxation_stress);
+
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=element->NewGauss(2);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
@@ -213,6 +239,7 @@
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
 		element->NodalFunctions(basis, gauss);
+		this->GetB(B,element,xyz_list,gauss);
 
 		/*Get all inputs on gauss point*/
 		thickness_input->GetInputValue(&thickness,gauss);
@@ -226,6 +253,9 @@
 		windvy_input->GetInputValue(&windvy,gauss);
 		oceanvx_input->GetInputValue(&oceanvx,gauss);
 		oceanvy_input->GetInputValue(&oceanvy,gauss);
+		sigma_xx_input->GetInputValue(&sigma_xx,gauss);
+		sigma_yy_input->GetInputValue(&sigma_yy,gauss);
+		sigma_xy_input->GetInputValue(&sigma_xy,gauss);
 		oceanssh_input->GetInputDerivativeValue(&dssh[0],xyz_list,gauss);
 
 		/*Previous speed (finite differencing)*/
@@ -238,34 +268,45 @@
 		vnorm    = sqrt(windvx*windvx + windvy*windvy);
 		constant_part = concentration*air_coef*rho_air*(air_lin_drag_coef+air_quad_drag_coef*vnorm);
 		for(int i=0;i<numnodes;i++){
-			pe->values[i*2+0] += dt*constant_part*windvx*Jdet*gauss->weight*basis[i];
-			pe->values[i*2+1] += dt*constant_part*windvy*Jdet*gauss->weight*basis[i];
+			pe->values[i*2+0] += constant_part*windvx*Jdet*gauss->weight*basis[i];
+			pe->values[i*2+1] += constant_part*windvy*Jdet*gauss->weight*basis[i];
 		}
 
 		/*Ocean forcing*/
 		vnorm = sqrt(pow(oceanvx-vx,2) + pow(oceanvy-vy,2));
 		constant_part = concentration*ocean_coef*rho_ocean*(ocean_lin_drag_coef+ocean_quad_drag_coef*vnorm);
 		for(int i=0;i<numnodes;i++){
-			pe->values[i*2+0] += dt*constant_part*(oceanvy-vy)*sin(ocean_turning_angle)*Jdet*gauss->weight*basis[i];
-			pe->values[i*2+1] += dt*constant_part*(vx-oceanvx)*sin(ocean_turning_angle)*Jdet*gauss->weight*basis[i];
+			pe->values[i*2+0] += constant_part*omega*(oceanvy-vy)*sin(ocean_turning_angle)*Jdet*gauss->weight*basis[i];
+			pe->values[i*2+1] += constant_part*omega*(vx-oceanvx)*sin(ocean_turning_angle)*Jdet*gauss->weight*basis[i];
 
-			pe->values[i*2+0] += dt*constant_part*oceanvx*cos(ocean_turning_angle)*Jdet*gauss->weight*basis[i];
-			pe->values[i*2+1] += dt*constant_part*oceanvy*cos(ocean_turning_angle)*Jdet*gauss->weight*basis[i];
+			pe->values[i*2+0] += constant_part*oceanvx*cos(ocean_turning_angle)*Jdet*gauss->weight*basis[i];
+			pe->values[i*2+1] += constant_part*oceanvy*cos(ocean_turning_angle)*Jdet*gauss->weight*basis[i];
+
+			pe->values[i*2+0] += (-rho_ice*thickness*gravity*omega*dssh[0])*Jdet*gauss->weight*basis[i];
+			pe->values[i*2+1] += (-rho_ice*thickness*gravity*omega*dssh[1])*Jdet*gauss->weight*basis[i];
 		}
 
 		/*Coriolis forces (use ustar)*/
 		for(int i=0;i<numnodes;i++){
-			pe->values[i*2+0] += dt*(-rho_ice*thickness*coriolis_factor*vystar)*Jdet*gauss->weight*basis[i];
-			pe->values[i*2+1] += dt*(+rho_ice*thickness*coriolis_factor*vxstar)*Jdet*gauss->weight*basis[i];
-
-			pe->values[i*2+0] += dt*(-rho_ice*thickness*gravity*dssh[0])*Jdet*gauss->weight*basis[i];
-			pe->values[i*2+1] += dt*(-rho_ice*thickness*gravity*dssh[1])*Jdet*gauss->weight*basis[i];
+			pe->values[i*2+0] += (-rho_ice*thickness*coriolis_factor*vystar)*Jdet*gauss->weight*basis[i];
+			pe->values[i*2+1] += (+rho_ice*thickness*coriolis_factor*vxstar)*Jdet*gauss->weight*basis[i];
 		}
+		
+		/*Add elastic part of previous time step*/
+		sigma_vec[0] = sigma_xx;
+		sigma_vec[1] = sigma_yy;
+		sigma_vec[2] = sigma_xy;
+		D[0][0] = D[1][1] = D[2][2] = time_ratio*thickness*Jdet*gauss->weight;
+		TripleMultiply(B,3,numdof,1,
+					&D[0][0],3,3,0,
+					&sigma_vec[0],3,1,0,
+					&pe->values[0],1);
 	}
 
 	/*Clean up and return*/
 	xDelete<IssmDouble>(xyz_list);
 	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(B);
 	delete gauss;
 	return pe;
 }/*}}}*/
@@ -322,9 +363,13 @@
 void SeaiceAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
 
 	/*Intermediaries*/
+	IssmDouble          c_min;
 	Vector<IssmDouble>* mask        = NULL;
 	IssmDouble*         serial_mask = NULL;
 
+	/*Get minimum concentration allowed*/
+	femmodel->parameters->FindParam(&c_min,SeaiceMinConcentrationEnum);
+
 	/*Step 1: update mask of active nodes*/
 	mask=new Vector<IssmDouble>(femmodel->nodes->NumberOfNodes(SeaiceAnalysisEnum));
 	for (int i=0;i<femmodel->elements->Size();i++){
@@ -333,7 +378,7 @@
 		/*Get current concentration of element and decide whether it is an active element*/
 		int    numnodes            = element->GetNumberOfNodes();
 		Input* concentration_input = element->GetInput(SeaiceConcentrationEnum); _assert_(concentration_input);
-		if(concentration_input->Max()>0.){
+		if(concentration_input->Max()>c_min){
 			for(int in=0;in<numnodes;in++) mask->SetValue(element->nodes[in]->Sid(),1.,INS_VAL);
 		}
 	}
@@ -378,7 +423,7 @@
 	IssmDouble ridging_exponent = element->GetMaterialParameter(MaterialsRidgingExponentEnum);
 
 	/*Get damage input at this location*/
-	Input* damage_input        = element->GetInput(MaterialsDamageEnum);     _assert_(damage_input);
+	Input* damage_input        = element->GetInput(DamageEnum);              _assert_(damage_input);
 	Input* thickness_input     = element->GetInput(SeaiceThicknessEnum);     _assert_(thickness_input);
 	Input* concentration_input = element->GetInput(SeaiceConcentrationEnum); _assert_(concentration_input);
 	damage_input->GetInputValue(&damage,gauss);
@@ -517,15 +562,12 @@
 		tan_phi          = element->GetMaterialParameter(MaterialsInternalFrictionCoefEnum); //mu = tan(phi)
 
 		/*Get current stress state*/
-		Input*   sigma_xx_input = element->GetInput(StressTensorxxEnum);
-		Input*   sigma_yy_input = element->GetInput(StressTensoryyEnum);
-		Input*   sigma_xy_input = element->GetInput(StressTensorxyEnum);
-		if(sigma_xx_input) sigma_xx_input->GetInputAverage(&sigma_xx);
-		else               sigma_xx = 0.;
-		if(sigma_yy_input) sigma_yy_input->GetInputAverage(&sigma_yy);
-		else               sigma_yy = 0.;
-		if(sigma_xy_input) sigma_xy_input->GetInputAverage(&sigma_xy);
-		else               sigma_xy = 0.;
+		Input*   sigma_xx_input = element->GetInput(StressTensorPredictorxxEnum); _assert_(sigma_xx_input);
+		Input*   sigma_yy_input = element->GetInput(StressTensorPredictoryyEnum); _assert_(sigma_yy_input);
+		Input*   sigma_xy_input = element->GetInput(StressTensorPredictorxyEnum); _assert_(sigma_xy_input);
+		sigma_xx_input->GetInputAverage(&sigma_xx);
+		sigma_yy_input->GetInputAverage(&sigma_yy);
+		sigma_xy_input->GetInputAverage(&sigma_xy);
 
 		/* Compute the invariants of the elastic deformation and instantaneous deformation rate */
 		sigma_s=sqrt(pow((sigma_xx-sigma_yy)/2.,2)+pow(sigma_xy,2));
@@ -536,11 +578,10 @@
 		compression_max=compression_coef*cohesion;
 
 		/* estimate the internal constraints using the current elastic deformation */
-		Input* damage_input    = element->GetInput(MaterialsDamageEnum); _assert_(damage_input);
-		Input* damagenew_input = element->GetInput(MaterialsDamageEnum); _assert_(damagenew_input);//FIXME
+		Input* damage_input    = element->GetInput(DamageEnum); _assert_(damage_input);
 		damage_input->GetInputAverage(&damage);
-		damagenew_input->GetInputAverage(&damage_new);
 		damage_test = damage;
+		damage_new  = damage;
 		if(sigma_n>traction || sigma_n<-compression_max){
 			if(sigma_n>traction){
 				sigma_target=traction;
@@ -563,7 +604,7 @@
 
 		/* The damage variable is changed */
 		damage_new=damage_test;
-		element->AddInput(MaterialsDamageEnum,&damage_test,P0Enum);//FIXME
+		element->AddInput(DamageEnum,&damage_new,P0Enum);
 
 		/* Recompute the internal stress*/
 		if(damage<1.){
Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18503)
+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18504)
@@ -49,6 +49,7 @@
 	BedEnum,
 	BaseEnum,
 	ConstantsGEnum,
+	ConstantsOmegaEnum,
 	ConstantsReferencetemperatureEnum,
 	ConstantsYtsEnum,
 	DependentObjectEnum,
@@ -197,6 +198,7 @@
 	DamageEquivStressEnum,
 	DamageEvolutionNumRequestedOutputsEnum,
 	DamageEvolutionRequestedOutputsEnum,
+	DamageEnum,
 	NewDamageEnum,
 	MaterialsRhoIceEnum,
 	MaterialsRhoSeawaterEnum,
@@ -748,6 +750,9 @@
 	SeaiceoceanEnum,
 	SeaiceThicknessEnum,
 	SeaiceConcentrationEnum,
+	SeaiceMinConcentrationEnum,
+	SeaiceMinThicknessEnum,
+	SeaiceMaxThicknessEnum,
 	SeaiceSpcvxEnum,
 	SeaiceSpcvyEnum,
 	SeaiceCoriolisFactorEnum,
@@ -768,7 +773,8 @@
 	MatseaiceEnum,
 	MaterialsPoissonEnum,
 	MaterialsYoungModulusEnum,
-	MaterialsDamageEnum,
+	MaterialsTimeRelaxationStressEnum,
+	MaterialsTimeRelaxationDamageEnum,
 	MaterialsRidgingExponentEnum,
 	MaterialsCohesionEnum,
 	MaterialsInternalFrictionCoefEnum,
@@ -776,6 +782,9 @@
 	MaterialsTractionCoefEnum,
 	VxStarEnum,
 	VyStarEnum,
+	StressTensorPredictorxxEnum,
+	StressTensorPredictoryyEnum,
+	StressTensorPredictorxyEnum,
 	/*}}}*/
 	MaximumNumberOfDefinitionsEnum
 };
Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18503)
+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18504)
@@ -57,6 +57,7 @@
 		case BedEnum : return "Bed";
 		case BaseEnum : return "Base";
 		case ConstantsGEnum : return "ConstantsG";
+		case ConstantsOmegaEnum : return "ConstantsOmega";
 		case ConstantsReferencetemperatureEnum : return "ConstantsReferencetemperature";
 		case ConstantsYtsEnum : return "ConstantsYts";
 		case DependentObjectEnum : return "DependentObject";
@@ -205,6 +206,7 @@
 		case DamageEquivStressEnum : return "DamageEquivStress";
 		case DamageEvolutionNumRequestedOutputsEnum : return "DamageEvolutionNumRequestedOutputs";
 		case DamageEvolutionRequestedOutputsEnum : return "DamageEvolutionRequestedOutputs";
+		case DamageEnum : return "Damage";
 		case NewDamageEnum : return "NewDamage";
 		case MaterialsRhoIceEnum : return "MaterialsRhoIce";
 		case MaterialsRhoSeawaterEnum : return "MaterialsRhoSeawater";
@@ -707,6 +709,9 @@
 		case SeaiceoceanEnum : return "Seaiceocean";
 		case SeaiceThicknessEnum : return "SeaiceThickness";
 		case SeaiceConcentrationEnum : return "SeaiceConcentration";
+		case SeaiceMinConcentrationEnum : return "SeaiceMinConcentration";
+		case SeaiceMinThicknessEnum : return "SeaiceMinThickness";
+		case SeaiceMaxThicknessEnum : return "SeaiceMaxThickness";
 		case SeaiceSpcvxEnum : return "SeaiceSpcvx";
 		case SeaiceSpcvyEnum : return "SeaiceSpcvy";
 		case SeaiceCoriolisFactorEnum : return "SeaiceCoriolisFactor";
@@ -727,7 +732,8 @@
 		case MatseaiceEnum : return "Matseaice";
 		case MaterialsPoissonEnum : return "MaterialsPoisson";
 		case MaterialsYoungModulusEnum : return "MaterialsYoungModulus";
-		case MaterialsDamageEnum : return "MaterialsDamage";
+		case MaterialsTimeRelaxationStressEnum : return "MaterialsTimeRelaxationStress";
+		case MaterialsTimeRelaxationDamageEnum : return "MaterialsTimeRelaxationDamage";
 		case MaterialsRidgingExponentEnum : return "MaterialsRidgingExponent";
 		case MaterialsCohesionEnum : return "MaterialsCohesion";
 		case MaterialsInternalFrictionCoefEnum : return "MaterialsInternalFrictionCoef";
@@ -735,6 +741,9 @@
 		case MaterialsTractionCoefEnum : return "MaterialsTractionCoef";
 		case VxStarEnum : return "VxStar";
 		case VyStarEnum : return "VyStar";
+		case StressTensorPredictorxxEnum : return "StressTensorPredictorxx";
+		case StressTensorPredictoryyEnum : return "StressTensorPredictoryy";
+		case StressTensorPredictorxyEnum : return "StressTensorPredictorxy";
 		case MaximumNumberOfDefinitionsEnum : return "MaximumNumberOfDefinitions";
 		default : return "unknown";
 
Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18503)
+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18504)
@@ -57,6 +57,7 @@
 	      else if (strcmp(name,"Bed")==0) return BedEnum;
 	      else if (strcmp(name,"Base")==0) return BaseEnum;
 	      else if (strcmp(name,"ConstantsG")==0) return ConstantsGEnum;
+	      else if (strcmp(name,"ConstantsOmega")==0) return ConstantsOmegaEnum;
 	      else if (strcmp(name,"ConstantsReferencetemperature")==0) return ConstantsReferencetemperatureEnum;
 	      else if (strcmp(name,"ConstantsYts")==0) return ConstantsYtsEnum;
 	      else if (strcmp(name,"DependentObject")==0) return DependentObjectEnum;
@@ -135,11 +136,11 @@
 	      else if (strcmp(name,"HydrologydcIsefficientlayer")==0) return HydrologydcIsefficientlayerEnum;
 	      else if (strcmp(name,"HydrologydcSedimentlimitFlag")==0) return HydrologydcSedimentlimitFlagEnum;
 	      else if (strcmp(name,"HydrologydcSedimentlimit")==0) return HydrologydcSedimentlimitEnum;
-	      else if (strcmp(name,"HydrologydcTransferFlag")==0) return HydrologydcTransferFlagEnum;
          else stage=2;
    }
    if(stage==2){
-	      if (strcmp(name,"HydrologydcLeakageFactor")==0) return HydrologydcLeakageFactorEnum;
+	      if (strcmp(name,"HydrologydcTransferFlag")==0) return HydrologydcTransferFlagEnum;
+	      else if (strcmp(name,"HydrologydcLeakageFactor")==0) return HydrologydcLeakageFactorEnum;
 	      else if (strcmp(name,"HydrologydcPenaltyFactor")==0) return HydrologydcPenaltyFactorEnum;
 	      else if (strcmp(name,"HydrologydcPenaltyLock")==0) return HydrologydcPenaltyLockEnum;
 	      else if (strcmp(name,"HydrologydcBasalMoulinInput")==0) return HydrologydcBasalMoulinInputEnum;
@@ -208,6 +209,7 @@
 	      else if (strcmp(name,"DamageEquivStress")==0) return DamageEquivStressEnum;
 	      else if (strcmp(name,"DamageEvolutionNumRequestedOutputs")==0) return DamageEvolutionNumRequestedOutputsEnum;
 	      else if (strcmp(name,"DamageEvolutionRequestedOutputs")==0) return DamageEvolutionRequestedOutputsEnum;
+	      else if (strcmp(name,"Damage")==0) return DamageEnum;
 	      else if (strcmp(name,"NewDamage")==0) return NewDamageEnum;
 	      else if (strcmp(name,"MaterialsRhoIce")==0) return MaterialsRhoIceEnum;
 	      else if (strcmp(name,"MaterialsRhoSeawater")==0) return MaterialsRhoSeawaterEnum;
@@ -257,12 +259,12 @@
 	      else if (strcmp(name,"MassFluxSegments")==0) return MassFluxSegmentsEnum;
 	      else if (strcmp(name,"MassFluxSegmentsPresent")==0) return MassFluxSegmentsPresentEnum;
 	      else if (strcmp(name,"QmuMassFluxSegmentsPresent")==0) return QmuMassFluxSegmentsPresentEnum;
-	      else if (strcmp(name,"QmuNumberofpartitions")==0) return QmuNumberofpartitionsEnum;
-	      else if (strcmp(name,"QmuNumberofresponses")==0) return QmuNumberofresponsesEnum;
          else stage=3;
    }
    if(stage==3){
-	      if (strcmp(name,"QmuPartition")==0) return QmuPartitionEnum;
+	      if (strcmp(name,"QmuNumberofpartitions")==0) return QmuNumberofpartitionsEnum;
+	      else if (strcmp(name,"QmuNumberofresponses")==0) return QmuNumberofresponsesEnum;
+	      else if (strcmp(name,"QmuPartition")==0) return QmuPartitionEnum;
 	      else if (strcmp(name,"QmuResponsedescriptors")==0) return QmuResponsedescriptorsEnum;
 	      else if (strcmp(name,"QmuVariabledescriptors")==0) return QmuVariabledescriptorsEnum;
 	      else if (strcmp(name,"RiftsNumrifts")==0) return RiftsNumriftsEnum;
@@ -380,12 +382,12 @@
 	      else if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum;
 	      else if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum;
 	      else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum;
-	      else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum;
-	      else if (strcmp(name,"SurfaceNormalVelocity")==0) return SurfaceNormalVelocityEnum;
          else stage=4;
    }
    if(stage==4){
-	      if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum;
+	      if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum;
+	      else if (strcmp(name,"SurfaceNormalVelocity")==0) return SurfaceNormalVelocityEnum;
+	      else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum;
 	      else if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum;
 	      else if (strcmp(name,"DepthAverageAnalysis")==0) return DepthAverageAnalysisEnum;
 	      else if (strcmp(name,"SteadystateSolution")==0) return SteadystateSolutionEnum;
@@ -503,12 +505,12 @@
 	      else if (strcmp(name,"Friction")==0) return FrictionEnum;
 	      else if (strcmp(name,"Internal")==0) return InternalEnum;
 	      else if (strcmp(name,"MassFlux")==0) return MassFluxEnum;
-	      else if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum;
-	      else if (strcmp(name,"Misfit")==0) return MisfitEnum;
          else stage=5;
    }
    if(stage==5){
-	      if (strcmp(name,"Pressure")==0) return PressureEnum;
+	      if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum;
+	      else if (strcmp(name,"Misfit")==0) return MisfitEnum;
+	      else if (strcmp(name,"Pressure")==0) return PressureEnum;
 	      else if (strcmp(name,"PressurePicard")==0) return PressurePicardEnum;
 	      else if (strcmp(name,"AndroidFrictionCoefficient")==0) return AndroidFrictionCoefficientEnum;
 	      else if (strcmp(name,"ResetPenalties")==0) return ResetPenaltiesEnum;
@@ -626,12 +628,12 @@
 	      else if (strcmp(name,"MisfitTimeinterpolation")==0) return MisfitTimeinterpolationEnum;
 	      else if (strcmp(name,"MisfitWeights")==0) return MisfitWeightsEnum;
 	      else if (strcmp(name,"MisfitWeightsEnum")==0) return MisfitWeightsEnumEnum;
-	      else if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum;
-	      else if (strcmp(name,"WeightsSurfaceObservation")==0) return WeightsSurfaceObservationEnum;
          else stage=6;
    }
    if(stage==6){
-	      if (strcmp(name,"MinVel")==0) return MinVelEnum;
+	      if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum;
+	      else if (strcmp(name,"WeightsSurfaceObservation")==0) return WeightsSurfaceObservationEnum;
+	      else if (strcmp(name,"MinVel")==0) return MinVelEnum;
 	      else if (strcmp(name,"MaxVel")==0) return MaxVelEnum;
 	      else if (strcmp(name,"MinVx")==0) return MinVxEnum;
 	      else if (strcmp(name,"MaxVx")==0) return MaxVxEnum;
@@ -722,6 +724,9 @@
 	      else if (strcmp(name,"Seaiceocean")==0) return SeaiceoceanEnum;
 	      else if (strcmp(name,"SeaiceThickness")==0) return SeaiceThicknessEnum;
 	      else if (strcmp(name,"SeaiceConcentration")==0) return SeaiceConcentrationEnum;
+	      else if (strcmp(name,"SeaiceMinConcentration")==0) return SeaiceMinConcentrationEnum;
+	      else if (strcmp(name,"SeaiceMinThickness")==0) return SeaiceMinThicknessEnum;
+	      else if (strcmp(name,"SeaiceMaxThickness")==0) return SeaiceMaxThicknessEnum;
 	      else if (strcmp(name,"SeaiceSpcvx")==0) return SeaiceSpcvxEnum;
 	      else if (strcmp(name,"SeaiceSpcvy")==0) return SeaiceSpcvyEnum;
 	      else if (strcmp(name,"SeaiceCoriolisFactor")==0) return SeaiceCoriolisFactorEnum;
@@ -742,16 +747,23 @@
 	      else if (strcmp(name,"Matseaice")==0) return MatseaiceEnum;
 	      else if (strcmp(name,"MaterialsPoisson")==0) return MaterialsPoissonEnum;
 	      else if (strcmp(name,"MaterialsYoungModulus")==0) return MaterialsYoungModulusEnum;
-	      else if (strcmp(name,"MaterialsDamage")==0) return MaterialsDamageEnum;
+	      else if (strcmp(name,"MaterialsTimeRelaxationStress")==0) return MaterialsTimeRelaxationStressEnum;
+	      else if (strcmp(name,"MaterialsTimeRelaxationDamage")==0) return MaterialsTimeRelaxationDamageEnum;
 	      else if (strcmp(name,"MaterialsRidgingExponent")==0) return MaterialsRidgingExponentEnum;
 	      else if (strcmp(name,"MaterialsCohesion")==0) return MaterialsCohesionEnum;
-	      else if (strcmp(name,"MaterialsInternalFrictionCoef")==0) return MaterialsInternalFrictionCoefEnum;
+         else stage=7;
+   }
+   if(stage==7){
+	      if (strcmp(name,"MaterialsInternalFrictionCoef")==0) return MaterialsInternalFrictionCoefEnum;
 	      else if (strcmp(name,"MaterialsCompressionCoef")==0) return MaterialsCompressionCoefEnum;
 	      else if (strcmp(name,"MaterialsTractionCoef")==0) return MaterialsTractionCoefEnum;
 	      else if (strcmp(name,"VxStar")==0) return VxStarEnum;
 	      else if (strcmp(name,"VyStar")==0) return VyStarEnum;
+	      else if (strcmp(name,"StressTensorPredictorxx")==0) return StressTensorPredictorxxEnum;
+	      else if (strcmp(name,"StressTensorPredictoryy")==0) return StressTensorPredictoryyEnum;
+	      else if (strcmp(name,"StressTensorPredictorxy")==0) return StressTensorPredictorxyEnum;
 	      else if (strcmp(name,"MaximumNumberOfDefinitions")==0) return MaximumNumberOfDefinitionsEnum;
-         else stage=7;
+         else stage=8;
    }
 	/*If we reach this point, the string provided has not been found*/
 	if(notfounderror)
Index: ../trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp	(revision 18503)
+++ ../trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp	(revision 18504)
@@ -79,7 +79,7 @@
 			}
 			break;
 		case MatseaiceEnum:
-			iomodel->FetchDataToInput(elements,MaterialsDamageEnum);
+			iomodel->FetchDataToInput(elements,DamageEnum);
 			for (i=0;i<iomodel->numberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matseaice(i+1,i,iomodel));
 			break;
 		default:
Index: ../trunk-jpl/src/c/cores/seaice_core.cpp
===================================================================
--- ../trunk-jpl/src/c/cores/seaice_core.cpp	(revision 18503)
+++ ../trunk-jpl/src/c/cores/seaice_core.cpp	(revision 18504)
@@ -19,11 +19,13 @@
 
 	/*Launch solution sequence for the only analysis we are interested in*/
 	femmodel->SetCurrentConfiguration(SeaiceAnalysisEnum);
+	SeaiceAnalysis* analysis = new SeaiceAnalysis();
+	analysis->UpdateDamageAndStress(femmodel);
 	solutionsequence_linear(femmodel);
 
 	if(save_results){
 		if(VerboseSolution()) _printf0_("   saving results\n");
-		int outputs[6] = {VxEnum,VyEnum,VelEnum,SeaiceConcentrationEnum,SeaiceThicknessEnum,MaterialsDamageEnum};
+		int outputs[6] = {VxEnum,VyEnum,VelEnum,SeaiceConcentrationEnum,SeaiceThicknessEnum,DamageEnum};
 		femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],6);
 	}
 
Index: ../trunk-jpl/src/c/classes/Materials/Matpar.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Materials/Matpar.cpp	(revision 18503)
+++ ../trunk-jpl/src/c/classes/Materials/Matpar.cpp	(revision 18504)
@@ -104,6 +104,9 @@
 			iomodel->Constant(&this->compression_coef,MaterialsCompressionCoefEnum);
 			iomodel->Constant(&this->traction_coef,MaterialsTractionCoefEnum);
 			iomodel->Constant(&this->g,ConstantsGEnum);
+			iomodel->Constant(&this->omega,ConstantsOmegaEnum);
+			iomodel->Constant(&this->time_relaxation_stress,MaterialsTimeRelaxationStressEnum);
+			iomodel->Constant(&this->time_relaxation_damage,MaterialsTimeRelaxationDamageEnum);
 			break;
 		default:
 			_error_("Material not supported yet");
@@ -318,6 +321,9 @@
 		case MaterialsInternalFrictionCoefEnum:		return this->internal_friction_coef;
 		case MaterialsCompressionCoefEnum:				return this->compression_coef;
 		case MaterialsTractionCoefEnum:					return this->traction_coef;
+		case ConstantsOmegaEnum:							return this->omega;
+		case MaterialsTimeRelaxationStressEnum:		return this->time_relaxation_stress;
+		case MaterialsTimeRelaxationDamageEnum:		return this->time_relaxation_damage;
 		default: _error_("Enum "<<EnumToStringx(enum_in)<<" not supported yet");
 	}
 
Index: ../trunk-jpl/src/c/classes/Materials/Matpar.h
===================================================================
--- ../trunk-jpl/src/c/classes/Materials/Matpar.h	(revision 18503)
+++ ../trunk-jpl/src/c/classes/Materials/Matpar.h	(revision 18504)
@@ -29,6 +29,7 @@
 		IssmDouble  mixed_layer_capacity;
 		IssmDouble  thermal_exchange_velocity;
 		IssmDouble  g;
+		IssmDouble  omega;
 		IssmDouble  desfac;
 		IssmDouble  s0p;
 
@@ -58,6 +59,8 @@
 		IssmDouble internal_friction_coef;
 		IssmDouble compression_coef;
 		IssmDouble traction_coef;
+		IssmDouble time_relaxation_stress;
+		IssmDouble time_relaxation_damage;
 
 	public:
 		Matpar();
