Index: /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp	(revision 19137)
+++ /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp	(revision 19138)
@@ -8,7 +8,28 @@
 void ThermalAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
 
+	/*Intermediary*/
+	int finiteelement = P1Enum;
+
 	/*Only 3d mesh supported*/
-	int finiteelement = P1Enum;
-	if(iomodel->domaintype==Domain3DEnum){
+	if(iomodel->domaintype!=Domain3DEnum) _error_("not supported yet");
+
+	/*Specific case for PDD, we want the constaints to be updated by the PDD scheme itself*/
+	bool isdynamic = false;
+	if(iomodel->solution_enum==ThermalSolutionEnum){
+		/*No PDD scheme, keep default*/
+	}
+	else if (iomodel->solution_enum==TransientSolutionEnum){
+		int smb_model;
+		iomodel->Constant(&smb_model,SurfaceforcingsEnum);
+		if(smb_model==SMBpddEnum) isdynamic=true;
+	}
+	else{
+		_error_("Solution "<<EnumToStringx(iomodel->solution_enum)<<" not supported yet");
+	}
+
+	if(isdynamic){
+		IoModelToDynamicConstraintsx(constraints,iomodel,ThermalSpctemperatureEnum,ThermalAnalysisEnum,finiteelement);
+	}
+	else{
 		IoModelToConstraintsx(constraints,iomodel,ThermalSpctemperatureEnum,ThermalAnalysisEnum,finiteelement);
 	}
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 19137)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 19138)
@@ -161,5 +161,5 @@
 	/*Open input file on cpu 0 and create IoModel */
 	if(my_rank==0) IOMODEL = pfopen0(inputfilename ,"rb");
-	IoModel* iomodel = new IoModel(IOMODEL,trace,X);
+	IoModel* iomodel = new IoModel(IOMODEL,in_solution_type,trace,X);
 
 	/*Figure out what analyses are activated for this solution*/
Index: /issm/trunk-jpl/src/c/classes/IoModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/IoModel.cpp	(revision 19137)
+++ /issm/trunk-jpl/src/c/classes/IoModel.cpp	(revision 19138)
@@ -21,4 +21,5 @@
 IoModel::IoModel(){/*{{{*/
 	this->fid=NULL;
+	this->solution_enum=-1;
 	this->data=NULL;
 	this->independents=NULL;
@@ -50,5 +51,5 @@
 }
 /*}}}*/
-IoModel::IoModel(FILE* iomodel_handle,bool trace,IssmPDouble* X){/*{{{*/
+IoModel::IoModel(FILE* iomodel_handle,int solution_enum_in,bool trace,IssmPDouble* X){/*{{{*/
 
 	bool autodiff=false;
@@ -60,4 +61,7 @@
 	/*Check that Enums are Synchronized*/
 	this->CheckEnumSync();
+
+	/*Keep track of solution*/
+	this->solution_enum = solution_enum_in;
 
 	/*Initialize data: */
@@ -240,5 +244,5 @@
 }
 /*}}}*/
-void IoModel::Constant(bool* poutput,int constant_enum){/*{{{*/
+void  IoModel::Constant(bool* poutput,int constant_enum){/*{{{*/
 
 	_assert_(constant_enum>=0);
@@ -248,5 +252,5 @@
 }
 /*}}}*/
-void IoModel::Constant(int* poutput,int constant_enum){/*{{{*/
+void  IoModel::Constant(int* poutput,int constant_enum){/*{{{*/
 
 	_assert_(constant_enum>=0);
@@ -256,5 +260,5 @@
 }
 /*}}}*/
-void IoModel::Constant(IssmDouble* poutput,int constant_enum){/*{{{*/
+void  IoModel::Constant(IssmDouble* poutput,int constant_enum){/*{{{*/
 
 	_assert_(constant_enum>=0);
@@ -264,5 +268,5 @@
 }
 /*}}}*/
-void IoModel::Constant(char** poutput,int constant_enum){/*{{{*/
+void  IoModel::Constant(char** poutput,int constant_enum){/*{{{*/
 
 	_assert_(constant_enum>=0);
@@ -291,5 +295,5 @@
 }
 /*}}}*/
-void IoModel::StartTrace(bool trace){/*{{{*/
+void  IoModel::StartTrace(bool trace){/*{{{*/
 
 	bool autodiff = false;
@@ -331,5 +335,5 @@
 }
 /*}}}*/
-void IoModel::DeclareIndependents(bool trace,IssmPDouble* X){/*{{{*/
+void  IoModel::DeclareIndependents(bool trace,IssmPDouble* X){/*{{{*/
 
 	int  i;
@@ -1516,5 +1520,5 @@
 }
 /*}}}*/
-void IoModel::FetchDataToInput(Elements* elements,int vector_enum,IssmDouble default_value){/*{{{*/
+void  IoModel::FetchDataToInput(Elements* elements,int vector_enum,IssmDouble default_value){/*{{{*/
 
 	/*intermediary: */
@@ -1541,5 +1545,5 @@
 }
 /*}}}*/
-void IoModel::FetchDataToInput(Elements* elements,int vector_enum){/*{{{*/
+void  IoModel::FetchDataToInput(Elements* elements,int vector_enum){/*{{{*/
 
 	/*intermediary: */
@@ -1613,5 +1617,5 @@
 }
 /*}}}*/
-void IoModel::LastIndex(int *pindex){/*{{{*/
+void  IoModel::LastIndex(int *pindex){/*{{{*/
 
 	int my_rank;
@@ -1649,6 +1653,5 @@
 }
 /*}}}*/
-/*FUNCTION IoModel::SetFilePointerToData{{{*/
-FILE* IoModel::SetFilePointerToData(int* pcode,int* pvector_type, int data_enum){
+FILE* IoModel::SetFilePointerToData(int* pcode,int* pvector_type, int data_enum){/*{{{*/
 
 	int my_rank;
Index: /issm/trunk-jpl/src/c/classes/IoModel.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/IoModel.h	(revision 19137)
+++ /issm/trunk-jpl/src/c/classes/IoModel.h	(revision 19138)
@@ -25,4 +25,7 @@
 		/*pointer to input file*/
 		FILE *fid;
+
+		/*Solution*/
+		int   solution_enum;
 
 		/*Partitioning*/
@@ -59,5 +62,5 @@
 		~IoModel();
 		IoModel();
-		IoModel(FILE* iomodel_handle,bool trace,IssmPDouble* X);
+		IoModel(FILE* iomodel_handle,int solution_enum_in,bool trace,IssmPDouble* X);
 
 		/*Input/Output*/
Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/transient_core.cpp	(revision 19137)
+++ /issm/trunk-jpl/src/c/cores/transient_core.cpp	(revision 19138)
@@ -23,14 +23,13 @@
 	bool       isstressbalance,ismasstransport,isFS,isthermal,isgroundingline,isgia,islevelset,isdamageevolution,ishydrology,iscalving;
 	bool       save_results,dakota_analysis;
-	bool       time_adapt=false;
+	bool       time_adapt;
 	int        output_frequency;
-	int        domaintype,groundingline_migration;
-	int        numoutputs         = 0;
-	Analysis  *analysis = NULL;
-	char**     requested_outputs = NULL;
-
+	int        domaintype,groundingline_migration,smb_model;
+	int        numoutputs;
+	Analysis  *analysis          = NULL;
+	char     **requested_outputs = NULL;
 
 	/*intermediary: */
-	int    step;
+	int        step;
 	IssmDouble time;
 
@@ -54,4 +53,5 @@
 	femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
 	femmodel->parameters->FindParam(&iscalving,TransientIscalvingEnum);
+	femmodel->parameters->FindParam(&smb_model,SurfaceforcingsEnum);
 	if(isgroundingline) femmodel->parameters->FindParam(&groundingline_migration,GroundinglineMigrationEnum);
 	femmodel->parameters->FindParam(&numoutputs,TransientNumRequestedOutputsEnum);
@@ -83,4 +83,5 @@
 		if(isthermal && domaintype==Domain3DEnum){
 			if(VerboseSolution()) _printf0_("   computing thermal regime\n");
+			if(smb_model==SMBpddEnum) ResetBoundaryConditions(femmodel,ThermalAnalysisEnum);
 			thermal_core(femmodel);
 		}