Index: /issm/trunk/src/c/Container/DataSet.cpp
===================================================================
--- /issm/trunk/src/c/Container/DataSet.cpp	(revision 9001)
+++ /issm/trunk/src/c/Container/DataSet.cpp	(revision 9002)
@@ -288,4 +288,11 @@
 				dataset->AddObject(spc);}
 				break;
+			case SpctEnum:{
+				Spct* spct=NULL;
+				spct=new Spct();
+				spct->Demarshall(&marshalled_dataset);
+				dataset->AddObject(spct);}
+				break;
+
 			case PengridEnum:{
 				Pengrid* pengrid=NULL;
Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 9001)
+++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 9002)
@@ -154,4 +154,5 @@
 	/*Spc: */
 	SpcEnum,
+	SpctEnum,
 	/*}}}*/
 	/*Geography {{{1*/
Index: /issm/trunk/src/c/Makefile.am
===================================================================
--- /issm/trunk/src/c/Makefile.am	(revision 9001)
+++ /issm/trunk/src/c/Makefile.am	(revision 9002)
@@ -230,6 +230,9 @@
 					./objects/Materials/Matpar.h\
 					./objects/Materials/Matpar.cpp\
+					./objects/Constraints/Constraint.h\
 					./objects/Constraints/Spc.cpp\
 					./objects/Constraints/Spc.h\
+					./objects/Constraints/Spct.cpp\
+					./objects/Constraints/Spct.h\
 					./objects/Loads/Penpair.cpp\
 					./objects/Loads/Penpair.h\
@@ -566,4 +569,6 @@
 					./modules/IoModelToInputsx/IoModelToInputsx.h\
 					./modules/IoModelToInputsx/IoModelToInputsx.cpp\
+					./modules/IoModelToConstraintsx/IoModelToConstraintsx.h\
+					./modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp\
 					./modules/InputUpdateFromConstantx/InputUpdateFromConstantx.h\
 					./modules/InputUpdateFromConstantx/InputUpdateFromConstantx.cpp\
@@ -596,4 +601,6 @@
 					./modules/SpcNodesx/SpcNodesx.h\
 					./modules/SpcNodesx/SpcNodesx.cpp\
+					./modules/UpdateConstraintsx/UpdateConstraintsx.h\
+					./modules/UpdateConstraintsx/UpdateConstraintsx.cpp\
 					./modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp\
 					./modules/InterpFromGridToMeshx/InterpFromGridToMeshx.h\
@@ -909,6 +916,9 @@
 					./objects/Materials/Matpar.h\
 					./objects/Materials/Matpar.cpp\
+					./objects/Constraints/Constraint.h\
 					./objects/Constraints/Spc.cpp\
 					./objects/Constraints/Spc.h\
+					./objects/Constraints/Spct.cpp\
+					./objects/Constraints/Spct.h\
 					./objects/Loads/Penpair.cpp\
 					./objects/Loads/Penpair.h\
@@ -1238,4 +1248,6 @@
 					./modules/IoModelToInputsx/IoModelToInputsx.h\
 					./modules/IoModelToInputsx/IoModelToInputsx.cpp\
+					./modules/IoModelToConstraintsx/IoModelToConstraintsx.h\
+					./modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp\
 					./modules/InputUpdateFromConstantx/InputUpdateFromConstantx.h\
 					./modules/InputUpdateFromConstantx/InputUpdateFromConstantx.cpp\
@@ -1268,4 +1280,6 @@
 					./modules/SpcNodesx/SpcNodesx.h\
 					./modules/SpcNodesx/SpcNodesx.cpp\
+					./modules/UpdateConstraintsx/UpdateConstraintsx.h\
+					./modules/UpdateConstraintsx/UpdateConstraintsx.cpp\
 					./modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp\
 					./modules/InterpFromGridToMeshx/InterpFromGridToMeshx.h\
Index: /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp
===================================================================
--- /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp	(revision 9001)
+++ /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp	(revision 9002)
@@ -130,4 +130,5 @@
 		case ExternalResultEnum : return "ExternalResult";
 		case SpcEnum : return "Spc";
+		case SpctEnum : return "Spct";
 		case GeographyEnum : return "Geography";
 		case IceSheetEnum : return "IceSheet";
Index: /issm/trunk/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp
===================================================================
--- /issm/trunk/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp	(revision 9002)
+++ /issm/trunk/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp	(revision 9002)
@@ -0,0 +1,111 @@
+/*!\file IoModelToConstraintsx
+ * \brief: used in ModelProcessor, to retrieve a nodal vector  and to create corresponding constraints using Spc objects
+ */
+
+#include "./IoModelToConstraintsx.h"
+#include "../../shared/shared.h"
+#include "../../io/io.h"
+#include "../../include/include.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+
+void IoModelToConstraintsx(Constraints* constraints,IoModel* iomodel,FILE* iomodel_handle,int vector_enum,int analysis_type){
+	
+
+	/*intermediary: */
+	int     i,j;
+	bool    transient     = false;
+
+	FILE   *fid           = NULL;
+	int     code          = 0;
+	int     vector_layout = 0;
+	int     counter;
+	int     nods;
+	double* times=NULL;
+	double* values=NULL;
+	bool    spcpresent=false;
+	int     count=0;
+
+	/*variables being fetched: */
+	double *doublevector  = NULL;
+	int     M,N;
+
+
+	/*First of, find the record for the enum, and get code  of data type: */
+	fid=SetFilePointerToData(&code, &vector_layout,iomodel_handle,vector_enum);
+
+	if(code!=7)_error_("%s%s"," expecting a double vector for constraints with enum ",EnumToStringx(vector_enum));
+	if(vector_layout!=1)_error_("%s%s"," expecting a nodal vector for constraints with enum ",EnumToStringx(vector_enum));
+
+	/*Fetch vector:*/
+	IoModelFetchData(&doublevector,&M,&N,iomodel_handle,vector_enum);
+
+	/*Transient or static?:*/
+	if(M==iomodel->numberofvertices){ //{{{1
+		/*static: just create Spc objects*/
+
+	
+		/*Initialize counter*/
+		count=0;
+	
+		/*Create spcs from x,y,z, as well as the spc values on those spcs: */
+		for (i=0;i<iomodel->numberofvertices;i++){
+			
+			/*keep only this partition's nodes:*/
+			if((iomodel->my_vertices[i])){
+
+				if (!isnan(doublevector[i])){
+
+					constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,doublevector[i],analysis_type));
+					count++;
+				}
+
+			} //if((my_nodes[i]==1))
+		}
+	} //}}}
+	else{ //{{{1
+		/*transient: create transient Spct objects. Same logic, except we need to retrieve 
+		 * various times and values to initialize an Spct object: */
+
+		/*Initialize counter*/
+		count=0;
+
+		/*figure out times: */
+		times=(double*)xmalloc(N*sizeof(double));
+		for(j=0;j<N;j++){
+			times[j]=doublevector[(M-1)*N+j];
+		}
+		/*unit conversion: */
+		UnitConversion(times,N,ExtToIuEnum,TimeEnum);
+
+
+
+		/*Create spcs from x,y,z, as well as the spc values on those spcs: */
+		for (i=0;i<iomodel->numberofvertices;i++){
+			
+			/*keep only this partition's nodes:*/
+			if((iomodel->my_vertices[i])){
+
+				/*figure out times and values: */
+				values=(double*)xmalloc(N*sizeof(double));
+				spcpresent=false;
+				for(j=0;j<N;j++){
+					values[j]=doublevector[i*N+j];
+					if(!isnan(values[j]))spcpresent=true; //NaN means no spc by default
+				}
+
+				if(spcpresent){
+					constraints->AddObject(new Spct(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,N,times,values,analysis_type));
+					count++;
+				}
+				xfree((void**)&values);
+			} //if((my_nodes[i]==1))
+		}
+	} //}}}
+
+	/*Free ressources:*/
+	xfree((void**)&doublevector);
+	xfree((void**)&times);
+	xfree((void**)&values);
+
+}
Index: /issm/trunk/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.h
===================================================================
--- /issm/trunk/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.h	(revision 9002)
+++ /issm/trunk/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.h	(revision 9002)
@@ -0,0 +1,13 @@
+/*!\file:  IoModelToConstraintsx.h
+ */ 
+
+#ifndef _IOMODEL_TO_CONSTRAINTS_H_
+#define _IOMODEL_TO_CONSTRAINTS_H_
+
+#include "../../Container/Container.h"
+#include "../../objects/objects.h"
+
+/* local prototypes: */
+void IoModelToConstraintsx(Constraints* constraints,IoModel* iomodel,FILE* iomodel_handle,int vector_enum,int analysis_type);
+
+#endif  /* _IOMODELTOELEMENTINPUTX_H */
Index: /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 9001)
+++ /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 9002)
@@ -73,4 +73,5 @@
 	parameters->AddObject(new BoolParam(IsprognosticEnum,iomodel->isprognostic));
 	parameters->AddObject(new BoolParam(IsthermalEnum,iomodel->isthermal));
+	parameters->AddObject(new DoubleParam(TimeEnum,0.0));  //start at time 0 by default for all solutions.
 
 	/*Deal with more complex parameters*/
Index: /issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp	(revision 9001)
+++ /issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp	(revision 9002)
@@ -4,4 +4,5 @@
 
 #include "../../../Container/Container.h"
+#include "../../../modules/modules.h"
 #include "../../../io/io.h"
 #include "../../../toolkits/toolkits.h"
@@ -19,5 +20,4 @@
 	/*Output*/
 	Constraints* constraints = NULL;
-	Spc*    spc  = NULL;
 
 	/*Recover pointer: */
@@ -30,30 +30,9 @@
 	if (iomodel->dim==2) goto cleanup_and_return;
 
-	/*Fetch data: */
-	IoModelFetchData(&iomodel->spctemperature,NULL,NULL,iomodel_handle,SpctemperatureEnum);
-
-	/*Initialize counter*/
-	count=0;
-
-	/*Create spcs from x,y,z, as well as the spc values on those spcs: */
-	for (i=0;i<iomodel->numberofvertices;i++){
-		/*keep only this partition's nodes:*/
-		if((iomodel->my_vertices[i])){
-
-			if ((int)iomodel->spctemperature[2*i]){
-
-				constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,iomodel->spctemperature[2*i+1],ThermalAnalysisEnum));
-				count++;
-
-			}
-
-		} //if((my_nodes[i]==1))
-	}
-
-	/*Free data: */
-	xfree((void**)&iomodel->spctemperature);
+	/*Create constraints: */
+	IoModelToConstraintsx(constraints,iomodel,iomodel_handle,SpctemperatureEnum,ThermalAnalysisEnum);
 
 	cleanup_and_return:
-	
+
 	/*Assign output pointer: */
 	*pconstraints=constraints;
Index: /issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp	(revision 9001)
+++ /issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp	(revision 9002)
@@ -39,5 +39,5 @@
 		if((iomodel->my_vertices[i]==1)){
 
-			if (!iomodel->spctemperature[2*i]){ //No penalty applied on spc nodes!
+			if (isnan(iomodel->spctemperature[i])){ //No penalty applied on spc nodes!
 
 				loads->AddObject(new Pengrid(iomodel->loadcounter+i+1,i,iomodel,ThermalAnalysisEnum));
Index: /issm/trunk/src/c/modules/SpcNodesx/SpcNodesx.cpp
===================================================================
--- /issm/trunk/src/c/modules/SpcNodesx/SpcNodesx.cpp	(revision 9001)
+++ /issm/trunk/src/c/modules/SpcNodesx/SpcNodesx.cpp	(revision 9002)
@@ -10,37 +10,17 @@
 #include "../../EnumDefinitions/EnumDefinitions.h"
 
-void SpcNodesx(Nodes* nodes,Constraints* constraints,int analysis_type){
+void SpcNodesx(Nodes* nodes,Constraints* constraints,Parameters* parameters, int analysis_type){
 
 	int i;
-	Node* node=NULL;
-	int nodeid;
-	int dof;
-	double value;
+			
+	for(i=0;i<constraints->Size();i++){
+	
+		Constraint* constraint=(Constraint*)constraints->GetObjectByOffset(i);
 
-	for(i=0;i<constraints->Size();i++){
+		/*Check this constraint belongs to this analysis: */
+		if(constraint->InAnalysis(analysis_type)){
 
-		Object* object=(Object*)constraints->GetObjectByOffset(i);
-
-		/*Check constraints is a single point constraint (spc): */
-		if(object->Enum()==SpcEnum){
-
-			Spc* spc=(Spc*)object;
-
-			if(spc->InAnalysis(analysis_type)){
-
-				/*Ok, constraints object is a constraint. Get the nodeid from the node it applies to: */
-				nodeid=spc->GetNodeId();
-				dof=spc->GetDof();
-				value=spc->GetValue();
-
-				/*Now, chase through nodes and find the corect node: */
-				node=(Node*)nodes->GetObjectById(NULL,nodeid);
-
-				/*Apply constraint: */
-				if(node){ //in case the spc is dealing with a node on another cpu
-					node->ApplyConstraint(dof,value);
-				}
-			}
-
+			/*Ok, apply constraint onto corresponding node: */
+			constraint->ConstrainNode(nodes,parameters);
 		}
 	}
Index: /issm/trunk/src/c/modules/SpcNodesx/SpcNodesx.h
===================================================================
--- /issm/trunk/src/c/modules/SpcNodesx/SpcNodesx.h	(revision 9001)
+++ /issm/trunk/src/c/modules/SpcNodesx/SpcNodesx.h	(revision 9002)
@@ -11,5 +11,5 @@
 
 /* local prototypes: */
-void SpcNodesx(Nodes* nodesin,Constraints* constraints,int analysis_type);
+void SpcNodesx(Nodes* nodes, Constraints* constraints,Parameters* parameters,int analysis_type);
 
 #endif  /* _SPCNODESX_H */
Index: /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp
===================================================================
--- /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp	(revision 9001)
+++ /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp	(revision 9002)
@@ -128,4 +128,5 @@
 	else if (strcmp(name,"ExternalResult")==0) return ExternalResultEnum;
 	else if (strcmp(name,"Spc")==0) return SpcEnum;
+	else if (strcmp(name,"Spct")==0) return SpctEnum;
 	else if (strcmp(name,"Geography")==0) return GeographyEnum;
 	else if (strcmp(name,"IceSheet")==0) return IceSheetEnum;
Index: /issm/trunk/src/c/modules/UpdateConstraintsx/UpdateConstraintsx.cpp
===================================================================
--- /issm/trunk/src/c/modules/UpdateConstraintsx/UpdateConstraintsx.cpp	(revision 9002)
+++ /issm/trunk/src/c/modules/UpdateConstraintsx/UpdateConstraintsx.cpp	(revision 9002)
@@ -0,0 +1,32 @@
+/*!\file UpdateConstraintsx
+ * \brief: update constraints and nodes (especially degree of freedom distribution) at time t.
+ */
+
+#include "./UpdateConstraintsx.h"
+
+#include "../../shared/shared.h"
+#include "../../include/include.h"
+#include "../../modules/modules.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+
+void UpdateConstraintsx(Nodes* nodes,Constraints* constraints,Parameters* parameters){
+
+	double time;
+	int    analysis_type;
+
+	/*retrieve parameters: */
+	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+	parameters->FindParam(&time,TimeEnum);
+
+	/*start module: */
+	_printf_(VerboseModule(),"%s%g\n","   Updating constraints for time: ",time);
+	
+	/*First, update dof constraints in nodes, using Spc constraints: */
+	SpcNodesx(nodes,constraints,parameters,analysis_type); 
+
+	/*Now, update degrees of freedoms: */
+	NodesDofx(nodes,parameters,analysis_type);
+
+
+}
Index: /issm/trunk/src/c/modules/UpdateConstraintsx/UpdateConstraintsx.h
===================================================================
--- /issm/trunk/src/c/modules/UpdateConstraintsx/UpdateConstraintsx.h	(revision 9002)
+++ /issm/trunk/src/c/modules/UpdateConstraintsx/UpdateConstraintsx.h	(revision 9002)
@@ -0,0 +1,15 @@
+/*!\file:  UpdateConstraintsx.h
+ * \brief header file for node and constraints update
+ */ 
+
+#ifndef _UPDATECONSTRAINTSX_H
+#define _UPDATECONSTRAINTSX_H
+
+#include "../../objects/objects.h"
+#include "../../Container/Container.h"
+#include "../../toolkits/toolkits.h"
+
+/* local prototypes: */
+void UpdateConstraintsx(Nodes* nodes,Constraints* constraints,Parameters* parameters);
+
+#endif  /* _UPDATECONSTRAINTSX_H */
Index: /issm/trunk/src/c/modules/modules.h
===================================================================
--- /issm/trunk/src/c/modules/modules.h	(revision 9001)
+++ /issm/trunk/src/c/modules/modules.h	(revision 9002)
@@ -50,4 +50,5 @@
 #include "./InputArtificialNoisex/InputArtificialNoisex.h"
 #include "./IoModelToInputsx/IoModelToInputsx.h"
+#include "./IoModelToConstraintsx/IoModelToConstraintsx.h"
 #include "./KMLMeshWritex/KMLFileReadx.h"
 #include "./KMLMeshWritex/KMLMeshWritex.h"
@@ -95,4 +96,5 @@
 #include "./Solverx/Solverx.h"
 #include "./SpcNodesx/SpcNodesx.h"
+#include "./UpdateConstraintsx/UpdateConstraintsx.h"
 #include "./SurfaceAreax/SurfaceAreax.h"
 #include "./SystemMatricesx/SystemMatricesx.h"
Index: /issm/trunk/src/c/objects/Constraints/Constraint.h
===================================================================
--- /issm/trunk/src/c/objects/Constraints/Constraint.h	(revision 9002)
+++ /issm/trunk/src/c/objects/Constraints/Constraint.h	(revision 9002)
@@ -0,0 +1,30 @@
+/*!\file:  Constraint.h
+ * \brief abstract class for Constraint object
+ * This class is a place holder for the Spc and Spct constraints
+ * It is derived from Object, so DataSets can contain them.
+ */ 
+
+
+#ifndef _CONSTRAINT_H_
+#define _CONSTRAINT_H_
+
+/*Headers:*/
+/*{{{1*/
+#include "../Object.h"
+
+class Nodes;
+
+#include "../../toolkits/toolkits.h"
+/*}}}*/
+
+class Constraint: public Object{
+
+	public: 
+		
+		virtual      ~Constraint(){};
+		virtual void ConstrainNode(Nodes* nodes,Parameters* parameters)=0;
+		virtual bool InAnalysis(int analysis_type)=0;
+
+
+};
+#endif
Index: /issm/trunk/src/c/objects/Constraints/Spc.cpp
===================================================================
--- /issm/trunk/src/c/objects/Constraints/Spc.cpp	(revision 9001)
+++ /issm/trunk/src/c/objects/Constraints/Spc.cpp	(revision 9002)
@@ -147,4 +147,26 @@
 /*}}}1*/
 
+/*Constraint virtual functions definitions: */
+/*FUNCTION Spc::InAnalysis{{{1*/
+bool Spc::InAnalysis(int in_analysis_type){
+	if (in_analysis_type==this->analysis_type) return true;
+	else return false;
+}
+/*}}}*/
+/*FUNCTION Spc::ConstrainNode{{{1*/
+void Spc::ConstrainNode(Nodes* nodes,Parameters* parameters){
+
+	Node* node=NULL;
+
+	/*Chase through nodes and find the node to which this Spc applys: */
+	node=(Node*)nodes->GetObjectById(NULL,nodeid);
+
+	/*Apply constraint: */
+	if(node){ //in case the spc is dealing with a node on another cpu
+		node->ApplyConstraint(dof,value);
+	}
+}
+/*}}}*/
+
 /*Spc functions*/
 /*FUNCTION Spc::GetDof {{{1*/
@@ -164,8 +186,2 @@
 }
 /*}}}1*/
-/*FUNCTION Spc::InAnalysis{{{1*/
-bool Spc::InAnalysis(int in_analysis_type){
-	if (in_analysis_type==this->analysis_type) return true;
-	else return false;
-}
-/*}}}*/
Index: /issm/trunk/src/c/objects/Constraints/Spc.h
===================================================================
--- /issm/trunk/src/c/objects/Constraints/Spc.h	(revision 9001)
+++ /issm/trunk/src/c/objects/Constraints/Spc.h	(revision 9002)
@@ -12,5 +12,5 @@
 /*}}}*/
 
-class Spc: public Object{
+class Spc: public Constraint{
 
 	private: 
@@ -39,9 +39,12 @@
 		Object* copy();
 		/*}}}*/
+		/*Constraint virtual functions definitions: {{{1*/
+		void   ConstrainNode(Nodes* nodes,Parameters* parameters);
+		bool   InAnalysis(int analysis_type);
+		/*}}}*/
 		/*Spc management:{{{1 */
 		int    GetNodeId();
 		int    GetDof();
 		double GetValue();
-		bool    InAnalysis(int analysis_type);
 		/*}}}*/
 
Index: /issm/trunk/src/c/objects/Constraints/Spct.cpp
===================================================================
--- /issm/trunk/src/c/objects/Constraints/Spct.cpp	(revision 9002)
+++ /issm/trunk/src/c/objects/Constraints/Spct.cpp	(revision 9002)
@@ -0,0 +1,246 @@
+/*!\file Spct.c
+ * \brief: implementation of the Spct object
+ */
+
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "stdio.h"
+#include <string.h>
+#include "../../include/include.h"
+#include "../../shared/shared.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../Container/Container.h"
+#include "../objects.h"
+
+/*Spct constructors and destructor*/
+/*FUNCTION Spct::Spct(){{{1*/
+Spct::Spct(){
+	sid=-1;
+	nodeid=-1;
+	dof=-1;
+	values=NULL;
+	times=NULL;
+	nsteps=-1;
+	analysis_type=-1;
+	return;
+}
+/*}}}1*/
+/*FUNCTION Spct::Spct(int spc_sid,int spc_nodeid,...){{{1*/
+Spct::Spct(int spc_sid,int spc_nodeid, int spc_dof,int spc_nsteps, double* spc_times, double* spc_values,int spc_analysis_type){
+
+	sid=spc_sid;
+	nodeid=spc_nodeid;
+	dof=spc_dof;
+	nsteps=spc_nsteps;
+	if(spc_nsteps){
+		values=(double*)xmalloc(spc_nsteps*sizeof(double));
+		times=(double*)xmalloc(spc_nsteps*sizeof(double));
+		memcpy(values,spc_values,nsteps*sizeof(double));
+		memcpy(times,spc_times,nsteps*sizeof(double));
+	}
+	analysis_type=spc_analysis_type;
+	return;
+}
+/*}}}1*/
+/*FUNCTION Spct::~Spct{{{1*/
+Spct::~Spct(){
+	xfree((void**)&times);
+	xfree((void**)&values);
+	return;
+}
+/*}}}1*/
+		
+/*Object virtual functions definitions:*/
+/*FUNCTION Spct::Echo {{{1*/
+void Spct::Echo(void){
+
+	int i;
+	printf("Spct:\n");
+	printf("   sid: %i\n",sid);
+	printf("   nodeid: %i\n",nodeid);
+	printf("   dof: %i\n",dof);
+	printf("   nsteps: %i\n",nsteps);
+	printf("   analysis_type: %s\n",EnumToStringx(analysis_type));
+	printf("   steps|times|values\n");
+	for(i=0;i<nsteps;i++){
+		printf("%i-%g:%g\n",i,times[i],values[i]);
+	}
+	return;
+}
+/*}}}1*/
+/*FUNCTION Spct::DeepEcho {{{1*/
+void Spct::DeepEcho(void){
+	this->Echo();
+}		
+/*}}}1*/
+/*FUNCTION Spct::Id {{{1*/
+int    Spct::Id(void){ return sid; }
+/*}}}1*/
+/*FUNCTION Spct::MyRank {{{1*/
+int    Spct::MyRank(void){ 
+	extern int my_rank;
+	return my_rank; 
+}
+/*}}}1*/
+/*FUNCTION Spct::Marshall {{{1*/
+void  Spct::Marshall(char** pmarshalled_dataset){
+
+	char* marshalled_dataset=NULL;
+	int   enum_type=0;
+
+	/*recover marshalled_dataset: */
+	marshalled_dataset=*pmarshalled_dataset;
+
+	/*get enum type of Spct: */
+	enum_type=SpctEnum;
+	
+	/*marshall enum: */
+	memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
+	
+	/*marshall Spct data: */
+	memcpy(marshalled_dataset,&sid,sizeof(sid));marshalled_dataset+=sizeof(sid);
+	memcpy(marshalled_dataset,&nodeid,sizeof(nodeid));marshalled_dataset+=sizeof(nodeid);
+	memcpy(marshalled_dataset,&dof,sizeof(dof));marshalled_dataset+=sizeof(dof);
+	memcpy(marshalled_dataset,&nsteps,sizeof(nsteps));marshalled_dataset+=sizeof(nsteps);
+	memcpy(marshalled_dataset,&analysis_type,sizeof(analysis_type));marshalled_dataset+=sizeof(analysis_type);
+	if(nsteps){
+		memcpy(marshalled_dataset,values,nsteps*sizeof(double));marshalled_dataset+=nsteps*sizeof(double);
+		memcpy(marshalled_dataset,times,nsteps*sizeof(double));marshalled_dataset+=nsteps*sizeof(double);
+	}
+
+	*pmarshalled_dataset=marshalled_dataset;
+	return;
+}
+/*}}}1*/
+/*FUNCTION Spct::MarshallSize {{{1*/
+int   Spct::MarshallSize(){
+
+	return sizeof(sid)
+		+sizeof(nodeid)
+		+sizeof(dof)
+		+sizeof(nsteps)
+		+nsteps*2*sizeof(double)
+		+sizeof(analysis_type)
+		+sizeof(int); //sizeof(int) for enum type
+}
+/*}}}1*/
+/*FUNCTION Spct::Demarshall {{{1*/
+void  Spct::Demarshall(char** pmarshalled_dataset){
+
+	char* marshalled_dataset=NULL;
+
+	/*recover marshalled_dataset: */
+	marshalled_dataset=*pmarshalled_dataset;
+
+	/*this time, no need to get enum type, the pointer directly points to the beginning of the 
+	 *object data (thanks to DataSet::Demarshall):*/
+
+	memcpy(&sid,marshalled_dataset,sizeof(sid));marshalled_dataset+=sizeof(sid);
+	memcpy(&nodeid,marshalled_dataset,sizeof(nodeid));marshalled_dataset+=sizeof(nodeid);
+	memcpy(&dof,marshalled_dataset,sizeof(dof));marshalled_dataset+=sizeof(dof);
+	memcpy(&nsteps,marshalled_dataset,sizeof(nsteps));marshalled_dataset+=sizeof(nsteps);
+	memcpy(&analysis_type,marshalled_dataset,sizeof(analysis_type));marshalled_dataset+=sizeof(analysis_type);
+	if(nsteps){
+		values=(double*)xmalloc(nsteps*sizeof(double));
+		times=(double*)xmalloc(nsteps*sizeof(double));
+		memcpy(values,marshalled_dataset,nsteps*sizeof(double));marshalled_dataset+=nsteps*sizeof(double);
+		memcpy(times,marshalled_dataset,nsteps*sizeof(double));marshalled_dataset+=nsteps*sizeof(double);
+	}
+
+	/*return: */
+	*pmarshalled_dataset=marshalled_dataset;
+	return;
+}
+/*}}}1*/
+/*FUNCTION Spct::Enum {{{1*/
+int Spct::Enum(void){
+
+	return SpctEnum;
+
+}
+/*}}}1*/
+/*FUNCTION Spct::copy {{{1*/
+Object* Spct::copy() {
+	return new Spct(sid,nodeid,dof,nsteps,times,values,analysis_type);
+}
+/*}}}1*/
+
+/*Constraint virtual functions definitions:*/
+/*FUNCTION Spct::InAnalysis{{{1*/
+bool Spct::InAnalysis(int in_analysis_type){
+	
+	if (in_analysis_type==this->analysis_type) return true;
+	else return false;
+}
+/*}}}*/
+/*FUNCTION Spct::ConstrainNode{{{1*/
+void Spct::ConstrainNode(Nodes* nodes,Parameters* parameters){
+
+	Node* node=NULL;
+	double time=0;
+	int    i;
+	double alpha=-1;
+	double value;
+	bool   found=false;
+
+	/*Chase through nodes and find the node to which this Spct applys: */
+	node=(Node*)nodes->GetObjectById(NULL,nodeid);
+	
+	if(node){ //in case the spc is dealing with a node on another cpu
+
+		/*Retrieve time in parameters: */
+		parameters->FindParam(&time,TimeEnum);
+
+		/*Now, go fetch value for this time: */
+		if (time<=times[0]){
+			value=values[0];
+			found=true;
+		}
+		else if (time>=times[nsteps-1]){
+			value=values[nsteps-1];
+			found=true;
+		}
+		else{
+			for(i=0;i<nsteps-1;i++){
+				if (times[i]<=time && time<times[i+1]){
+					alpha=(time-times[i])/(times[i+1]-times[i]);
+					value=(1-alpha)*values[i]+alpha*values[i+1];
+					found=true;
+					break;
+				}
+			}
+		}
+
+		if(!found)_error_("could not find time segment for constraint");
+
+		/*Apply or relax constraint: */
+		if(isnan(value)){
+			node->RelaxConstraint(dof);
+		}
+		else node->ApplyConstraint(dof,value);
+	}
+}
+/*}}}*/
+
+/*Spct functions*/
+/*FUNCTION Spct::GetDof {{{1*/
+int Spct::GetDof(){
+	return dof;
+}
+/*}}}1*/
+/*FUNCTION Spct::GetNodeId {{{1*/
+int   Spct::GetNodeId(){
+	
+	return nodeid;
+}
+/*}}}1*/
+/*FUNCTION Spct::GetValue {{{1*/
+double Spct::GetValue(){
+	return values[0];
+}
+/*}}}1*/
+
Index: /issm/trunk/src/c/objects/Constraints/Spct.h
===================================================================
--- /issm/trunk/src/c/objects/Constraints/Spct.h	(revision 9002)
+++ /issm/trunk/src/c/objects/Constraints/Spct.h	(revision 9002)
@@ -0,0 +1,55 @@
+/*!\file Spct.h
+ * \brief: header file for spc object
+ */
+
+#ifndef _SPCT_H_
+#define _SPCT_H_
+
+/*Headers:*/
+/*{{{1*/
+#include "../Object.h"
+class DataSet;
+/*}}}*/
+
+class Spct: public Constraint{
+
+	private: 
+		int	sid; /*! id, to track it*/
+		int	nodeid; /*!node id*/
+		int dof; /*!component*/
+		double* values; /*different values in time*/
+		double* times; /*different time steps*/
+		int nsteps; /*number of time steps*/
+		int analysis_type;
+
+	public:
+
+		/*Spct constructors, destructors:{{{1*/
+		Spct();
+		Spct(int sid,int nodeid, int dof,int nsteps, double* times, double* values,int analysis_type);
+		~Spct();
+		/*}}}*/
+		/*Object virtual functions definitions:{{{1 */
+		void  Echo();
+		void  DeepEcho();
+		int   Id(); 
+		int   MyRank();
+		void  Marshall(char** pmarshalled_dataset);
+		int   MarshallSize();
+		void  Demarshall(char** pmarshalled_dataset);
+		int   Enum();
+		Object* copy();
+		/*}}}*/
+		/*Constraint virtual functions definitions: {{{1*/
+		void   ConstrainNode(Nodes* nodes,Parameters* parameters);
+		bool   InAnalysis(int analysis_type);
+		/*}}}*/
+		/*Spct management:{{{1 */
+		int    GetNodeId();
+		int    GetDof();
+		double GetValue();
+		/*}}}*/
+
+};
+
+#endif  /* _SPCT_H_ */
Index: /issm/trunk/src/c/objects/FemModel.cpp
===================================================================
--- /issm/trunk/src/c/objects/FemModel.cpp	(revision 9001)
+++ /issm/trunk/src/c/objects/FemModel.cpp	(revision 9002)
@@ -59,6 +59,6 @@
 
 		_printf_(VerboseMProcessor(),"      resolve node constraints\n");
-		SpcNodesx(nodes,constraints,analysis_type); 
-	
+		SpcNodesx(nodes,constraints,parameters,analysis_type); 
+
 		_printf_(VerboseMProcessor(),"      create nodal degrees of freedom\n");
 		NodesDofx(nodes,parameters,analysis_type);
Index: /issm/trunk/src/c/objects/Loads/Load.h
===================================================================
--- /issm/trunk/src/c/objects/Loads/Load.h	(revision 9001)
+++ /issm/trunk/src/c/objects/Loads/Load.h	(revision 9002)
@@ -23,4 +23,5 @@
 
 		virtual       ~Load(){};
+		
 		/*Virtual functions: {{{1*/
 		virtual void  Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0;
Index: /issm/trunk/src/c/objects/Node.cpp
===================================================================
--- /issm/trunk/src/c/objects/Node.cpp	(revision 9001)
+++ /issm/trunk/src/c/objects/Node.cpp	(revision 9002)
@@ -529,4 +529,12 @@
 }
 /*}}}*/
+/*FUNCTION Node::RelaxConstraint{{{1*/
+void  Node::RelaxConstraint(int dof){
+
+	/*Dof should be added to the f-set, and taken out of the s-set:*/
+	DofInFSet(dof-1);
+	this->indexing.svalues[dof-1]=NAN;
+}
+/*}}}*/
 /*FUNCTION Node::CreateVecSets {{{1*/
 void  Node::CreateVecSets(Vec pv_g,Vec pv_f,Vec pv_s){
@@ -593,4 +601,14 @@
 }
 /*}}}*/
+/*FUNCTION Node::DofInFSet {{{1*/
+void  Node::DofInFSet(int dof){
+
+	/*Put dof for this node into the f set (ie, this dof will NOT be constrained 
+	 * to a fixed value during computations. */
+
+	this->indexing.f_set[dof]=1; 
+	this->indexing.s_set[dof]=0;
+}
+/*}}}*/
 /*FUNCTION Node::FreezeDof{{{1*/
 void  Node::FreezeDof(int dof){
Index: /issm/trunk/src/c/objects/Node.h
===================================================================
--- /issm/trunk/src/c/objects/Node.h	(revision 9001)
+++ /issm/trunk/src/c/objects/Node.h	(revision 9002)
@@ -75,5 +75,7 @@
 		int   IsClone();
 		void  ApplyConstraint(int dof,double value);
+		void  RelaxConstraint(int dof);
 		void  DofInSSet(int dof);
+		void  DofInFSet(int dof);
 		int   GetDof(int dofindex,int setenum);
 		void  CreateVecSets(Vec pv_g,Vec pv_f,Vec pv_s);
Index: /issm/trunk/src/c/objects/objects.h
===================================================================
--- /issm/trunk/src/c/objects/objects.h	(revision 9001)
+++ /issm/trunk/src/c/objects/objects.h	(revision 9002)
@@ -20,5 +20,7 @@
 
 /*Constraints: */
+#include "./Constraints/Constraint.h"
 #include "./Constraints/Spc.h"
+#include "./Constraints/Spct.h"
 
 /*Gauss*/
Index: /issm/trunk/src/c/shared/Numerics/UnitConversion.cpp
===================================================================
--- /issm/trunk/src/c/shared/Numerics/UnitConversion.cpp	(revision 9001)
+++ /issm/trunk/src/c/shared/Numerics/UnitConversion.cpp	(revision 9002)
@@ -47,4 +47,5 @@
 	double scale;
 	switch(type_enum){
+		case TimeEnum:        scale=1.0/yts;break; //yr
 		case VxEnum:          scale=yts;break; //m/yr
 		case VxObsEnum:       scale=yts;break; //m/yr
Index: /issm/trunk/src/c/solutions/thermal_core.cpp
===================================================================
--- /issm/trunk/src/c/solutions/thermal_core.cpp	(revision 9001)
+++ /issm/trunk/src/c/solutions/thermal_core.cpp	(revision 9002)
@@ -43,4 +43,5 @@
 		if(nsteps)_printf_(VerboseSolution(),"time step: %i/%i\n",i+1,nsteps);
 		time=(i+1)*dt;
+		femmodel->parameters->SetParam(time,TimeEnum);
 
 		/*call thermal_core_step: */
Index: /issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp
===================================================================
--- /issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp	(revision 9001)
+++ /issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp	(revision 9002)
@@ -48,4 +48,7 @@
 	InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,reset_penalties,ResetPenaltiesEnum);
 	InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,false,ConvergedEnum);
+
+	UpdateConstraintsx(femmodel->nodes,femmodel->constraints,femmodel->parameters);
+
 	for(;;){
 
Index: /issm/trunk/src/m/enum/SpctEnum.m
===================================================================
--- /issm/trunk/src/m/enum/SpctEnum.m	(revision 9002)
+++ /issm/trunk/src/m/enum/SpctEnum.m	(revision 9002)
@@ -0,0 +1,11 @@
+function macro=SpctEnum()
+%SPCTENUM - Enum of Spct
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh
+%            Please read src/c/EnumDefinitions/README for more information
+%
+%   Usage:
+%      macro=SpctEnum()
+
+macro=StringToEnum('Spct');
Index: /issm/trunk/src/m/model/extrude.m
===================================================================
--- /issm/trunk/src/m/model/extrude.m	(revision 9001)
+++ /issm/trunk/src/m/model/extrude.m	(revision 9002)
@@ -136,49 +136,41 @@
 
 %drag_coefficient is limited to nodes that are on the bedrock.
-md.drag_coefficient=project3d(md,md.drag_coefficient,'node',1);
+md.drag_coefficient=project3d(md,'vector',md.drag_coefficient,'type','node','layer',1);
 
 %p and q (same deal, except for element that are on the bedrock: )
-md.drag_p=project3d(md,md.drag_p,'element');
-md.drag_q=project3d(md,md.drag_q,'element');
+md.drag_p=project3d(md,'vector',md.drag_p,'type','element');
+md.drag_q=project3d(md,'vector',md.drag_q,'type','element');
 
 %observations
-md.vx_obs=project3d(md,md.vx_obs,'node');
-md.vy_obs=project3d(md,md.vy_obs,'node');
-md.vel_obs=project3d(md,md.vel_obs,'node');
-md.vx_bal=project3d(md,md.vx_bal,'node');
-md.vy_bal=project3d(md,md.vy_bal,'node');
-md.vel_bal=project3d(md,md.vel_bal,'node');
-md.vel_obs_raw=project3d(md,md.vel_obs_raw,'node');
-md.surface_mass_balance=project3d(md,md.surface_mass_balance,'node');
-md.surface_accumulation_rate=project3d(md,md.surface_accumulation_rate,'node');
-md.surface_ablation_rate=project3d(md,md.surface_ablation_rate,'node');
-md.dhdt=project3d(md,md.dhdt,'node');
-md.firn_layer=project3d(md,md.firn_layer,'node',md.numlayers);
-if ~isempty(md.forcings),
-	forcingnames=fieldnames(md.forcings);
-	numforcings=length(fieldnames(md.forcings));
-	for i=1:numforcings
-		forcing=md.forcings.(forcingnames{i});
-		md.forcings.(forcingnames{i})=[project3d(md,forcing(1:end-1,:),'node');forcing(end,:)];
-	end
-end
+md.vx_obs=project3d(md,'vector',md.vx_obs,'type','node');
+md.vy_obs=project3d(md,'vector',md.vy_obs,'type','node');
+md.vel_obs=project3d(md,'vector',md.vel_obs,'type','node');
+md.vx_bal=project3d(md,'vector',md.vx_bal,'type','node');
+md.vy_bal=project3d(md,'vector',md.vy_bal,'type','node');
+md.vel_bal=project3d(md,'vector',md.vel_bal,'type','node');
+md.vel_obs_raw=project3d(md,'vector',md.vel_obs_raw,'type','node');
+md.surface_mass_balance=project3d(md,'vector',md.surface_mass_balance,'type','node');
+md.surface_accumulation_rate=project3d(md,'vector',md.surface_accumulation_rate,'type','node');
+md.surface_ablation_rate=project3d(md,'vector',md.surface_ablation_rate,'type','node');
+md.dhdt=project3d(md,'vector',md.dhdt,'type','node');
+md.firn_layer=project3d(md,'vector',md.firn_layer,'type','node','layer',md.numlayers);
 
 %results
-if ~isnan(md.vx),md.vx=project3d(md,md.vx,'node');end;
-if ~isnan(md.vy),md.vy=project3d(md,md.vy,'node');end;
-if ~isnan(md.vz),md.vz=project3d(md,md.vz,'node');end;
-if ~isnan(md.vel),md.vel=project3d(md,md.vel,'node');end;
-if ~isnan(md.temperature),md.temperature=project3d(md,md.temperature,'node');end;
-if ~isnan(md.waterfraction),md.waterfraction=project3d(md,md.waterfraction,'node');end;
-if ~isnan(md.surface_slopex),md.surface_slopex=project3d(md,md.surface_slopex,'node');end;
-if ~isnan(md.surface_slopey),md.surface_slopey=project3d(md,md.surface_slopey,'node');end;
-if ~isnan(md.bed_slopex),md.bed_slopex=project3d(md,md.bed_slopex,'node');end;
-if ~isnan(md.bed_slopey),md.bed_slopey=project3d(md,md.bed_slopey,'node');end;
+if ~isnan(md.vx),md.vx=project3d(md,'vector',md.vx,'type','node');end;
+if ~isnan(md.vy),md.vy=project3d(md,'vector',md.vy,'type','node');end;
+if ~isnan(md.vz),md.vz=project3d(md,'vector',md.vz,'type','node');end;
+if ~isnan(md.vel),md.vel=project3d(md,'vector',md.vel,'type','node');end;
+if ~isnan(md.temperature),md.temperature=project3d(md,'vector',md.temperature,'type','node');end;
+if ~isnan(md.waterfraction),md.waterfraction=project3d(md,'vector',md.waterfraction,'type','node');end;
+if ~isnan(md.surface_slopex),md.surface_slopex=project3d(md,'vector',md.surface_slopex,'type','node');end;
+if ~isnan(md.surface_slopey),md.surface_slopey=project3d(md,'vector',md.surface_slopey,'type','node');end;
+if ~isnan(md.bed_slopex),md.bed_slopex=project3d(md,'vector',md.bed_slopex,'type','node');end;
+if ~isnan(md.bed_slopey),md.bed_slopey=project3d(md,'vector',md.bed_slopey,'type','node');end;
 
 %bedinfo and surface info
-md.elementonbed=project3d(md,ones(md.numberofelements2d,1),'element',1);
-md.elementonsurface=project3d(md,ones(md.numberofelements2d,1),'element',md.numlayers-1);
-md.nodeonbed=project3d(md,ones(md.numberofnodes2d,1),'node',1);
-md.nodeonsurface=project3d(md,ones(md.numberofnodes2d,1),'node',md.numlayers);
+md.elementonbed=project3d(md,'vector',ones(md.numberofelements2d,1),'type','element','layer',1);
+md.elementonsurface=project3d(md,'vector',ones(md.numberofelements2d,1),'type','element','layer',md.numlayers-1);
+md.nodeonbed=project3d(md,'vector',ones(md.numberofnodes2d,1),'type','node','layer',1);
+md.nodeonsurface=project3d(md,'vector',ones(md.numberofnodes2d,1),'type','node','layer',md.numlayers);
 
 %elementstype
@@ -186,9 +178,9 @@
 	oldelements_type=md.elements_type2d;
 	md.elements_type=zeros(number_el3d,1);
-	md.elements_type=project3d(md,oldelements_type,'element');
-	md.nodeonhutter=project3d(md,md.nodeonhutter,'node');
-	md.nodeonmacayeal=project3d(md,md.nodeonmacayeal,'node');
-	md.nodeonpattyn=project3d(md,md.nodeonpattyn,'node');
-	md.nodeonstokes=project3d(md,md.nodeonstokes,'node');
+	md.elements_type=project3d(md,'vector',oldelements_type,'type','element');
+	md.nodeonhutter=project3d(md,'vector',md.nodeonhutter,'type','node');
+	md.nodeonmacayeal=project3d(md,'vector',md.nodeonmacayeal,'type','node');
+	md.nodeonpattyn=project3d(md,'vector',md.nodeonpattyn,'type','node');
+	md.nodeonstokes=project3d(md,'vector',md.nodeonstokes,'type','node');
 end
 
@@ -197,14 +189,14 @@
 	oldvertices_type=md.vertices_type2d;
 	md.vertices_type=zeros(number_nodes3d,1);
-	md.vertices_type=project3d(md,oldvertices_type,'node');
+	md.vertices_type=project3d(md,'vector',oldvertices_type,'type','node');
 end
 
 %boundary conditions
-md.spcvx=project3d(md,md.spcvx,'node');
-md.spcvy=project3d(md,md.spcvy,'node');
-md.spcvz=project3d(md,md.spcvz,'node');
-md.spctemperature=project3d(md,md.spctemperature,'node',md.numlayers);
-md.spcthickness=project3d(md,md.spcthickness,'node');
-md.diagnostic_ref=project3d(md,md.diagnostic_ref,'node');
+md.spcvx=project3d(md,'vector',md.spcvx,'type','node');
+md.spcvy=project3d(md,'vector',md.spcvy,'type','node');
+md.spcvz=project3d(md,'vector',md.spcvz,'type','node');
+md.spctemperature=project3d(md,'vector',md.spctemperature,'type','node','layer',md.numlayers,'padding',NaN);
+md.spcthickness=project3d(md,'vector',md.spcthickness,'type','node');
+md.diagnostic_ref=project3d(md,'vector',md.diagnostic_ref,'type','node');
 
 %in 3d, pressureload: [node1 node2 node3 node4 element]
@@ -226,19 +218,19 @@
 
 %materials
-md.rheology_B=project3d(md,md.rheology_B,'node');
-md.rheology_n=project3d(md,md.rheology_n,'element');
+md.rheology_B=project3d(md,'vector',md.rheology_B,'type','node');
+md.rheology_n=project3d(md,'vector',md.rheology_n,'type','element');
 
 %parameters
-md.surface=project3d(md,md.surface,'node');
-md.thickness=project3d(md,md.thickness,'node');
-md.bed=project3d(md,md.bed,'node');
-md.nodeonboundary=project3d(md,md.nodeonboundary,'node');
-md.elementoniceshelf=project3d(md,md.elementoniceshelf,'element');
-md.nodeoniceshelf=project3d(md,md.nodeoniceshelf,'node');
-md.elementonicesheet=project3d(md,md.elementonicesheet,'element');
-md.nodeonicesheet=project3d(md,md.nodeonicesheet,'node');
-md.elementonwater=project3d(md,md.elementonwater,'element');
-md.nodeonwater=project3d(md,md.nodeonwater,'node');
-if ~isnan(md.weights),md.weights=project3d(md,md.weights,'node');end;
+md.surface=project3d(md,'vector',md.surface,'type','node');
+md.thickness=project3d(md,'vector',md.thickness,'type','node');
+md.bed=project3d(md,'vector',md.bed,'type','node');
+md.nodeonboundary=project3d(md,'vector',md.nodeonboundary,'type','node');
+md.elementoniceshelf=project3d(md,'vector',md.elementoniceshelf,'type','element');
+md.nodeoniceshelf=project3d(md,'vector',md.nodeoniceshelf,'type','node');
+md.elementonicesheet=project3d(md,'vector',md.elementonicesheet,'type','element');
+md.nodeonicesheet=project3d(md,'vector',md.nodeonicesheet,'type','node');
+md.elementonwater=project3d(md,'vector',md.elementonwater,'type','element');
+md.nodeonwater=project3d(md,'vector',md.nodeonwater,'type','node');
+if ~isnan(md.weights),md.weights=project3d(md,'vector',md.weights,'type','node');end;
 
 %Put lithostatic pressure is there is an existing pressure
@@ -248,12 +240,12 @@
 
 %special for thermal modeling:
-md.basal_melting_rate=project3d(md,md.basal_melting_rate,'node',1); 
-md.observed_temperature=project3d(md,md.observed_temperature,'node'); 
+md.basal_melting_rate=project3d(md,'vector',md.basal_melting_rate,'type','node','layer',1); 
+md.observed_temperature=project3d(md,'vector',md.observed_temperature,'type','node'); 
 if ~isnan(md.geothermalflux)
-	md.geothermalflux=project3d(md,md.geothermalflux,'node',1); %bedrock only gets geothermal flux
+	md.geothermalflux=project3d(md,'vector',md.geothermalflux,'type','node','layer',1); %bedrock only gets geothermal flux
 end
 
 %mask:
-md.mask=project3d(md,md.mask,'node');
+md.mask=project3d(md,'vector',md.mask,'type','node');
 
 %increase connectivity if less than 25:
Index: /issm/trunk/src/m/model/ismodelselfconsistent.m
===================================================================
--- /issm/trunk/src/m/model/ismodelselfconsistent.m	(revision 9001)
+++ /issm/trunk/src/m/model/ismodelselfconsistent.m	(revision 9002)
@@ -485,5 +485,5 @@
 
 			%CHECK THAT WE ARE NOT FULLY CONSTRAINED
-			if isnan(find(~md.spctemperature(:,1))),
+			if ~any(~isnan(md.spctemperature))
 				message(['model not consistent: model ' md.name ' is totally constrained for temperature, no need to solve!']);
 			end
@@ -507,6 +507,12 @@
 
 				%INITIAL TEMPERATURE
-				fields={'temperature','spctemperature(:,2)','observed_temperature'};
+				fields={'temperature','observed_temperature'};
 				checkgreater(md,fields,0)
+
+				%CHECK SPCTEMPERATURE that are not NaN are >0.
+				if find(any(md.spctemperature(find(~isnan(md.spctemperature)))<=0)),
+					message(['model not consistent: model ' md.name ' is constrained with negative or nil temperatures!']);
+				end
+
 			end
 			%}}}
Index: /issm/trunk/src/m/model/project3d.m
===================================================================
--- /issm/trunk/src/m/model/project3d.m	(revision 9001)
+++ /issm/trunk/src/m/model/project3d.m	(revision 9002)
@@ -1,3 +1,3 @@
-function projected_vector=project3d(md3d,vector2d,type,varargin);
+function projected_vector=project3d(md3d,varargin);
 %PROJECT3D - vertically project a vector from 2d mesh
 %
@@ -5,23 +5,34 @@
 %   This vector can be a node vector of size (md3d.numberofnodes2d,N/A) or an 
 %   element vector of size (md3d.numberofelements2d,N/A). 
-%   type is 'element' or 'node'. layer an optional layer number where vector 
-%   should keep its values. If not specified, all layers adopt the value of the 
-%   2d vector.
+%   arguments: 
+%      'vector': 2d vector
+%      'type': 'element' or 'node'. 
+%   options: 
+%      'layer' a layer number where vector should keep its values. If not specified, all layers adopt the 
+%             value of the 2d vector.
+%      'padding': default to 0 (value adopted by other 3d layers not being projected0
 %
-%   Usage:
-%      extruded_vector=project3d(md3d,vector2d,type,layer);
+%   Egs:
+%      extruded_vector=project3d(md3d,'vector',vector2d,'type','node','layer',1,'padding',NaN);
+%      extruded_vector=project3d(md3d,'vector',vector2d,'type','element','padding',0);
+%      extruded_vector=project3d(md3d,'vector',vector2d,'type','node');
 
-if nargin==4,
-	layer=varargin{1};
-	if ((layer<1) || (layer>md3d.numlayers)),
-		error(['project3d error message: layer shoud be between 1 and ' num2str(md3d.numlayers)]);
-	end
-else
-	layer=0;
+
+if nargin==0,
+	Project3dUsage;
+	error();
 end
+
+options=pairoptions(varargin{:});
+
+%retrieve parameters from options.
+vector2d=getfieldvalue(options,'vector'); %mandatory
+type=getfieldvalue(options,'type'); %mandatory
+layer=getfieldvalue(options,'layer',0); %optional (do all layers otherwise)
+paddingvalue=getfieldvalue(options,'padding',0); %0 by default
 
 if strcmpi(type,'node'),
 
-	projected_vector=zeros(md3d.numberofnodes,size(vector2d,2));
+	projected_vector=paddingvalue*ones(md3d.numberofnodes,size(vector2d,2));
 	
 	if layer==0,
@@ -32,7 +43,7 @@
 		projected_vector(((layer-1)*md3d.numberofnodes2d+1):(layer*md3d.numberofnodes2d),:)=vector2d;
 	end
-else
+elseif strcmpi(type,'element'),
 
-	projected_vector=zeros(md3d.numberofelements,size(vector2d,2));
+	projected_vector=paddingvalue*ones(md3d.numberofelements,size(vector2d,2));
 
 	if layer==0,
@@ -44,3 +55,26 @@
 		projected_vector( ((layer-1)*md3d.numberofelements2d+1):(layer*md3d.numberofelements2d),:)=vector2d;
 	end
+else
+	error('project3d error message: unknown projection type');
 end
+
+function Project3dUsage;
+
+disp('function projected_vector=project3d(md3d,varargin);');
+disp('PROJECT3D - vertically project a vector from 2d mesh...');
+disp('');
+disp('   vertically project a vector from 2d mesh (split in noncoll and coll areas) into a 3d mesh.');
+disp('   This vector can be a node vector of size (md3d.numberofnodes2d,N/A) or an ');
+disp('   element vector of size (md3d.numberofelements2d,N/A). ');
+disp('   arguments: ');
+disp('      ''vector'': 2d vector');
+disp('      ''type'': ''element'' or ''node''. ');
+disp('   options: ');
+disp('      ''layer'' a layer number where vector should keep its values. If not specified, all layers adopt the ');
+disp('             value of the 2d vector.');
+disp('      ''padding'': default to 0 (value adopted by other 3d layers not being projected0');
+disp('');
+disp('   Egs:');
+disp('      extruded_vector=project3d(md3d,''vector'',vector2d,''type'',''node'',''layer'',1,''padding'',NaN);');
+disp('      extruded_vector=project3d(md3d,''vector'',vector2d,''type'',''element'',''padding'',0);');
+disp('      extruded_vector=project3d(md3d,''vector'',vector2d,''type'',''node'');');
Index: /issm/trunk/src/m/solutions/NewFemModel.m
===================================================================
--- /issm/trunk/src/m/solutions/NewFemModel.m	(revision 9001)
+++ /issm/trunk/src/m/solutions/NewFemModel.m	(revision 9002)
@@ -37,5 +37,5 @@
 
 		issmprintf(VerboseMProcessor(),'%s','      resolve node constraints');
-		[femmodel.nodes]=SpcNodes(femmodel.nodes,femmodel.constraints,analysis_type);
+		[femmodel.nodes]=SpcNodes(femmodel.nodes,femmodel.constraints,femmodel.parameters,analysis_type);
 
 		issmprintf(VerboseMProcessor(),'%s','      create nodal degrees of freedom');
Index: /issm/trunk/src/m/solvers/solver_thermal_nonlinear.m
===================================================================
--- /issm/trunk/src/m/solvers/solver_thermal_nonlinear.m	(revision 9001)
+++ /issm/trunk/src/m/solvers/solver_thermal_nonlinear.m	(revision 9002)
@@ -20,4 +20,7 @@
 	[femmodel.elements loads]=InputUpdateFromConstant(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,0,ConvergedEnum);
 	[femmodel.elements femmodel.loads]=InputUpdateFromConstant(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,1,ResetPenaltiesEnum);
+
+	%update constraints
+	[femmodel.nodes]=UpdateConstraints(femmodel.nodes,femmodel.constraints,femmodel.parameters);
 
 	while(~converged),
Index: /issm/trunk/src/m/utils/BC/SetIceSheetBC.m
===================================================================
--- /issm/trunk/src/m/utils/BC/SetIceSheetBC.m	(revision 9001)
+++ /issm/trunk/src/m/utils/BC/SetIceSheetBC.m	(revision 9002)
@@ -61,5 +61,6 @@
 
 if (length(md.observed_temperature)==md.numberofnodes),
-	md.spctemperature=[md.nodeonsurface md.observed_temperature]; %impose observed temperature on surface
+	md.spctemperature=NaN*ones(md.numberofnodes,1);
+	pos=find(md.nodeonsurface); md.spctemperature(pos)=md.observed_temperature(pos); %impose observed temperature on surface
 	if (length(md.geothermalflux)~=md.numberofnodes),
 		md.geothermalflux=50*10^-3*ones(md.numberofnodes,1); %50 mW/m^2
Index: /issm/trunk/src/m/utils/BC/SetIceShelfBC.m
===================================================================
--- /issm/trunk/src/m/utils/BC/SetIceShelfBC.m	(revision 9001)
+++ /issm/trunk/src/m/utils/BC/SetIceShelfBC.m	(revision 9002)
@@ -93,5 +93,6 @@
 
 if (length(md.observed_temperature)==md.numberofnodes),
-	md.spctemperature=[md.nodeonsurface md.observed_temperature]; %impose observed temperature on surface
+	md.spctemperature=NaN*ones(md.numberofnodes,1);
+	pos=find(md.nodeonsurface); md.spctemperature(pos)=md.observed_temperature(pos); %impose observed temperature on surface
 	if (length(md.geothermalflux)~=md.numberofnodes),
 		md.geothermalflux=zeros(md.numberofnodes,1);
Index: /issm/trunk/src/m/utils/BC/SetMarineIceSheetBC.m
===================================================================
--- /issm/trunk/src/m/utils/BC/SetMarineIceSheetBC.m	(revision 9001)
+++ /issm/trunk/src/m/utils/BC/SetMarineIceSheetBC.m	(revision 9002)
@@ -104,5 +104,6 @@
 
 if (length(md.observed_temperature)==md.numberofnodes),
-	md.spctemperature=[md.nodeonsurface md.observed_temperature]; %impose observed temperature on surface
+	md.spctemperature=NaN*ones(md.numberofnodes,1);
+	pos=find(md.nodeonsurface); md.spctemperature(pos)=md.observed_temperature(pos); %impose observed temperature on surface
 	if (length(md.geothermalflux)~=md.numberofnodes),
 		md.geothermalflux=zeros(md.numberofnodes,1);
Index: /issm/trunk/src/mex/Makefile.am
===================================================================
--- /issm/trunk/src/mex/Makefile.am	(revision 9001)
+++ /issm/trunk/src/mex/Makefile.am	(revision 9002)
@@ -69,4 +69,5 @@
 				SparseToVector\
 				SpcNodes\
+				UpdateConstraints\
 				SurfaceArea\
 				SystemMatrices\
@@ -301,4 +302,7 @@
 			  SpcNodes/SpcNodes.h
 
+UpdateConstraints_SOURCES = UpdateConstraints/UpdateConstraints.cpp\
+			  UpdateConstraints/UpdateConstraints.h
+
 SystemMatrices_SOURCES = SystemMatrices/SystemMatrices.cpp\
 			  SystemMatrices/SystemMatrices.h
Index: /issm/trunk/src/mex/SpcNodes/SpcNodes.cpp
===================================================================
--- /issm/trunk/src/mex/SpcNodes/SpcNodes.cpp	(revision 9001)
+++ /issm/trunk/src/mex/SpcNodes/SpcNodes.cpp	(revision 9002)
@@ -10,4 +10,5 @@
 	/*input datasets: */
 	Nodes* nodes=NULL;
+	Parameters* parameters=NULL;
 	Constraints* constraints=NULL;
 	int      analysis_type;
@@ -22,8 +23,9 @@
 	FetchMatlabData((DataSet**)&nodes,NODESIN);
 	FetchMatlabData((DataSet**)&constraints,CONSTRAINTS);
+	FetchMatlabData(&parameters,PARAMETERS);
 	FetchMatlabData(&analysis_type,ANALYSISTYPE);
 	
 	/*!Generate internal degree of freedom numbers: */
-	SpcNodesx(nodes,constraints,analysis_type); 
+	SpcNodesx(nodes,constraints,parameters,analysis_type); 
 
 	/*write output datasets: */
@@ -33,4 +35,5 @@
 	delete nodes;
 	delete constraints;
+	delete parameters;
 
 	/*end module: */
@@ -41,5 +44,5 @@
 {
 	_printf_(true,"\n");
-	_printf_(true,"   usage: [m.node]=%s(m.nodes,m.constraints);\n",__FUNCT__);
+	_printf_(true,"   usage: [m.node]=%s(m.nodes,m.constraints,m.parameters);\n",__FUNCT__);
 	_printf_(true,"\n");
 }
Index: /issm/trunk/src/mex/SpcNodes/SpcNodes.h
===================================================================
--- /issm/trunk/src/mex/SpcNodes/SpcNodes.h	(revision 9001)
+++ /issm/trunk/src/mex/SpcNodes/SpcNodes.h	(revision 9002)
@@ -21,5 +21,6 @@
 #define NODESIN (mxArray*)prhs[0]
 #define CONSTRAINTS (mxArray*)prhs[1]
-#define ANALYSISTYPE (mxArray*)prhs[2]
+#define PARAMETERS (mxArray*)prhs[2]
+#define ANALYSISTYPE (mxArray*)prhs[3]
 
 /* serial output macros: */
@@ -30,5 +31,5 @@
 #define NLHS  1
 #undef NRHS
-#define NRHS  3
+#define NRHS  4
 
 
Index: /issm/trunk/src/mex/UpdateConstraints/UpdateConstraints.cpp
===================================================================
--- /issm/trunk/src/mex/UpdateConstraints/UpdateConstraints.cpp	(revision 9002)
+++ /issm/trunk/src/mex/UpdateConstraints/UpdateConstraints.cpp	(revision 9002)
@@ -0,0 +1,46 @@
+/*\file UpdateConstraints.c
+ *\brief: build degrees of freedom for every node.
+ */
+
+#include "./UpdateConstraints.h"
+
+
+void mexFunction( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]){
+
+	/*input datasets: */
+	Nodes* nodes=NULL;
+	Parameters* parameters=NULL;
+	Constraints* constraints=NULL;
+
+	/*Boot module: */
+	MODULEBOOT();
+
+	/*checks on arguments on the matlab side: */
+	CheckNumMatlabArguments(nlhs,NLHS,nrhs,NRHS,__FUNCT__,&UpdateConstraintsUsage);
+
+	/*Input datasets: */
+	FetchMatlabData((DataSet**)&nodes,NODESIN);
+	FetchMatlabData((DataSet**)&constraints,CONSTRAINTS);
+	FetchMatlabData(&parameters,PARAMETERS);
+	
+	/*!Generate internal degree of freedom numbers: */
+	UpdateConstraintsx(nodes,constraints,parameters);
+
+	/*write output datasets: */
+	WriteMatlabData(NODES,nodes);
+
+	/*Free ressources: */
+	delete nodes;
+	delete constraints;
+	delete parameters;
+
+	/*end module: */
+	MODULEEND();
+}
+
+void UpdateConstraintsUsage(void)
+{
+	_printf_(true,"\n");
+	_printf_(true,"   usage: [m.node]=%s(m.nodes,m.constraints,m.parameters);\n",__FUNCT__);
+	_printf_(true,"\n");
+}
Index: /issm/trunk/src/mex/UpdateConstraints/UpdateConstraints.h
===================================================================
--- /issm/trunk/src/mex/UpdateConstraints/UpdateConstraints.h	(revision 9002)
+++ /issm/trunk/src/mex/UpdateConstraints/UpdateConstraints.h	(revision 9002)
@@ -0,0 +1,36 @@
+
+/*
+	UpdateConstraints.h
+*/
+
+
+#ifndef _UPDATECONSTRAINTS_H
+#define _UPDATECONSTRAINTS_H
+
+/* local prototypes: */
+void UpdateConstraintsUsage(void);
+
+#include "../../c/modules/modules.h"
+#include "../../c/Container/Container.h"
+#include "../../c/shared/shared.h"
+
+#undef __FUNCT__ 
+#define __FUNCT__  "UpdateConstraints"
+
+/* serial input macros: */
+#define NODESIN (mxArray*)prhs[0]
+#define CONSTRAINTS (mxArray*)prhs[1]
+#define PARAMETERS (mxArray*)prhs[2]
+
+/* serial output macros: */
+#define NODES (mxArray**)&plhs[0]
+
+/* serial arg counts: */
+#undef NLHS
+#define NLHS  1
+#undef NRHS
+#define NRHS  3
+
+
+#endif  /* _UPDATECONSTRAINTS_H */
+