Index: issm/trunk-jpl/src/c/classes/IoModel.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/IoModel.cpp	(revision 20645)
+++ issm/trunk-jpl/src/c/classes/IoModel.cpp	(revision 20653)
@@ -3,4 +3,12 @@
  * into ISSM, from Matlab, or through a binary file opened for reading.
  */
+
+/*CODES:
+ * 1: boolean constant
+ * 2: integer constant
+ * 3: IssmDouble constant
+ * 5: boolean vector
+ * 6: int vector
+ * 7: IssmDouble vector*/
 
 #ifdef HAVE_CONFIG_H
@@ -19,11 +27,364 @@
 #include "../shared/shared.h"
 
+/*NEW NEW NEW NEW NEW NEW NEW NEW NEW NEW*/
+IoConstant::IoConstant(){/*{{{*/
+	this->isindependent = false;
+	this->data_enum     = -1;
+	this->name          = NULL;
+	this->constant      = NULL;
+}
+/*}}}*/
+IoConstant::~IoConstant(){/*{{{*/
+	xDelete<char>(this->name);
+	delete this->constant;
+}
+/*}}}*/
+IoConstant::IoConstant(bool value,int enum_in){/*{{{*/
+	this->isindependent = false;
+	this->data_enum     = enum_in;
+	this->name          = NULL;
+	this->constant      = new BoolParam(enum_in,value);
+}
+/*}}}*/
+IoConstant::IoConstant(int value,int enum_in){/*{{{*/
+	this->isindependent = false;
+	this->data_enum     = enum_in;
+	this->name          = NULL;
+	this->constant      = new IntParam(enum_in,value);
+}
+/*}}}*/
+IoConstant::IoConstant(IssmDouble value,int enum_in){/*{{{*/
+	this->isindependent = false;
+	this->data_enum     = enum_in;
+	this->name          = NULL;
+	this->constant      = new DoubleParam(enum_in,value);
+}
+/*}}}*/
+IoConstant::IoConstant(char* value,int enum_in){/*{{{*/
+	this->isindependent = false;
+	this->data_enum     = enum_in;
+	this->name          = NULL;
+	this->constant      = new StringParam(enum_in,value);
+}
+/*}}}*/
+
+IoData::IoData(){/*{{{*/
+	this->isindependent = false;
+	this->data_enum     = -1;
+	this->name          = NULL;
+	this->code          = -1;
+	this->layout        = -1;
+	this->M             = 0;
+	this->N             = 0;
+	this->data          = NULL;
+}
+/*}}}*/
+IoData::~IoData(){/*{{{*/
+	xDelete<char>(this->name);
+	xDelete<IssmDouble>(this->data);
+}
+/*}}}*/
+IoData::IoData(IssmDouble* matrix,int code_in,int layout_in,int M_in,int N_in,int enum_in){/*{{{*/
+	this->isindependent = false;
+	this->data_enum     = enum_in;
+	this->name          = NULL;
+	this->code          = code_in;
+	this->layout        = layout_in;
+	this->M             = M_in;
+	this->N             = N_in;
+	this->data          = matrix; /*do not copy*/
+	_assert_(code_in==5 ||  code_in==6 || code_in==7);
+}
+/*}}}*/
+
+void  IoModel::AddConstant(IoConstant* in_constant){/*{{{*/
+
+	_assert_(in_constant);
+
+	/*Go through dataset of constant and check whether it already exists */
+	vector<IoConstant*>::iterator iter;
+
+	for(iter=constants.begin();iter<constants.end();iter++){
+		if((*iter)->data_enum==in_constant->data_enum){
+			delete in_constant;
+			return;
+		}
+	}
+
+	this->constants.push_back(in_constant);
+}
+/*}}}*/
+void  IoModel::AddConstantIndependent(IoConstant* in_constant){/*{{{*/
+
+	_assert_(in_constant);
+
+	/*Set constant as independent*/
+	in_constant->isindependent = true;
+
+	/*Add to constnats*/
+	this->AddConstant(in_constant);
+}
+/*}}}*/
+void  IoModel::AddData(IoData* in_data){/*{{{*/
+
+	_assert_(in_data);
+
+	/*Go through dataset of data and check whether it already exists */
+	vector<IoData*>::iterator iter;
+
+	for(iter=data.begin();iter<data.end();iter++){
+		if((*iter)->data_enum==in_data->data_enum){
+			delete in_data;
+			return;
+		}
+	}
+
+	this->data.push_back(in_data);
+}
+/*}}}*/
+void  IoModel::AddDataIndependent(IoData* in_data){/*{{{*/
+
+	_assert_(in_data);
+
+	/*Set data as independent*/
+	in_data->isindependent = true;
+
+	/*Add to constnats*/
+	this->AddData(in_data);
+}
+/*}}}*/
+void  IoModel::FetchIndependentConstant(int* pXcount,IssmPDouble* X,int name){/*{{{*/
+
+	/*recover my_rank:*/
+	int my_rank=IssmComm::GetRank();
+
+	/*recover Xcount if X is not NULL:*/
+	int Xcount = 0;
+	if(X) Xcount=*pXcount;
+
+	#ifdef _HAVE_ADOLC_ //cannot come here unless you are running AD mode, from DeclaredIndependents:
+
+	/*output: */
+	IssmPDouble  pscalar;
+	IssmDouble   scalar; //same as pscalar, except it's an ADOLC independent variable
+	int          code;
+
+	/*Set file pointer to beginning of the data: */
+	fid=this->SetFilePointerToData(&code,NULL,name);
+	if(code!=3) _error_("expecting a IssmDouble for enum " << EnumToStringx(name));
+
+	/*We have to read a scalar from disk. First read the dimensions of the scalar, then the scalar: */
+	if(my_rank==0){
+		if(fread(&pscalar,sizeof(IssmPDouble),1,fid)!=1)_error_("could not read scalar ");
+
+		/*Now, before we even broadcast this to other nodes, declare the scalar  as an independent variable!. If we 
+		 *have been supplied an X vector, use it instead of what we just read: */
+		if(X){
+			scalar<<=X[Xcount];
+		}
+		else{
+			scalar<<=pscalar;
+		}
+	}
+
+	ISSM_MPI_Bcast(&scalar,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 
+	this->AddConstantIndependent(new IoConstant(scalar,name));
+
+	/*increment offset into X vector, now that we have read 1 value:*/
+	Xcount++; *pXcount=Xcount;
+	#endif
+}
+/*}}}*/
+void  IoModel::FetchIndependentData(int* pXcount,IssmPDouble* X,int name){/*{{{*/
+
+	/*recover my_rank:*/
+	int my_rank=IssmComm::GetRank();
+
+	/*recover Xcount if X is not NULL:*/
+	int Xcount = 0;
+	if(X) Xcount=*pXcount;
+
+	#ifdef _HAVE_ADOLC_ //cannot come here unless you are running AD mode, from DeclaredIndependents:
+
+	/*Intermediaries*/
+	int M,N;
+	IssmPDouble* buffer=NULL; //a buffer to read the data from disk
+	IssmDouble* matrix=NULL; //our independent variable
+	int code,layout;
+
+	/*Set file pointer to beginning of the data: */
+	fid=this->SetFilePointerToData(&code,&layout,name);
+	if((code!=5) && (code!=6) && (code!=7))_error_("expecting a IssmDouble, integer or boolean matrix for enum " << EnumToStringx(name));
+
+	/*We have to read a matrix from disk. First read the dimensions of the matrix, then the whole matrix: */
+	/*numberofelements: */
+	if(my_rank==0){  
+		if(fread(&M,sizeof(int),1,fid)!=1) _error_("could not read number of rows for matrix ");
+	}
+	ISSM_MPI_Bcast(&M,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 
+
+	if(my_rank==0){  
+		if(fread(&N,sizeof(int),1,fid)!=1) _error_("could not read number of columns for matrix ");
+	}
+	ISSM_MPI_Bcast(&N,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 
+
+	/*Now allocate matrix: */
+	if(M*N){
+		buffer=xNew<IssmPDouble>(M*N);
+		matrix=xNew<IssmDouble>(M*N);
+
+		/*Read matrix on node 0, then broadcast: */
+		if(my_rank==0){  
+			if(fread(buffer,M*N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read matrix ");
+
+			/*Now, before we even broadcast this to other nodes, declare the whole matrix as a independent variable!
+			  If we have been supplied an X vector, use it instead of what we just read: */
+			if(X){
+				for(int i=0;i<M*N;i++) matrix[i]<<=X[Xcount+i];  /*<<= ADOLC overloaded operator to declare independent*/
+			}
+			else{
+				for(int i=0;i<M*N;i++) matrix[i]<<=buffer[i];
+			}
+		}
+		ISSM_MPI_Bcast(matrix,M*N,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 
+
+		xDelete<IssmPDouble>(buffer);
+	}
+	else _error_("cannot declare the independent variable " << EnumToStringx(name) <<  "if it's empty!");
+
+	/*Add to data as independent*/
+	this->AddDataIndependent(new IoData(matrix,code,layout,M,N,name));
+
+	/*increment offset into X vector, now that we have read M*N values:*/
+	Xcount+=M*N; *pXcount=Xcount;
+	#endif
+}
+/*}}}*/
+void  IoModel::FindConstant(bool* pvalue,int constant_enum){/*{{{*/
+
+	/*Intermediary*/
+	vector<IoConstant*>::iterator iter;
+
+	for(iter=constants.begin();iter<constants.end();iter++){
+		IoConstant* ioconstant=*iter;
+
+		if(ioconstant->data_enum==constant_enum){
+			ioconstant->constant->GetParameterValue(pvalue);
+			return;
+		}
+	}
+
+	for(vector<IoConstant*>::iterator iter=constants.begin();iter<constants.end();iter++) (*iter)->constant->Echo();
+	_error_("Could not find constant \""<<EnumToStringx(constant_enum) <<"\"");
+}
+/*}}}*/
+void  IoModel::FindConstant(int* pvalue,int constant_enum){/*{{{*/
+
+	/*Intermediary*/
+	vector<IoConstant*>::iterator iter;
+
+	for(iter=constants.begin();iter<constants.end();iter++){
+		IoConstant* ioconstant=*iter;
+
+		if(ioconstant->data_enum==constant_enum){
+			ioconstant->constant->GetParameterValue(pvalue);
+			return;
+		}
+	}
+
+	_error_("Could not find constant \""<<EnumToStringx(constant_enum) <<"\"");
+}
+/*}}}*/
+void  IoModel::FindConstant(IssmDouble* pvalue,int constant_enum){/*{{{*/
+
+	/*Intermediary*/
+	vector<IoConstant*>::iterator iter;
+
+	for(iter=constants.begin();iter<constants.end();iter++){
+		IoConstant* ioconstant=*iter;
+
+		if(ioconstant->data_enum==constant_enum){
+			ioconstant->constant->GetParameterValue(pvalue);
+			return;
+		}
+	}
+
+	_error_("Could not find constant \""<<EnumToStringx(constant_enum) <<"\"");
+}
+/*}}}*/
+void  IoModel::FindConstant(char** pvalue,int constant_enum){/*{{{*/
+
+	/*Intermediary*/
+	vector<IoConstant*>::iterator iter;
+
+	for(iter=constants.begin();iter<constants.end();iter++){
+		IoConstant* ioconstant=*iter;
+
+		if(ioconstant->data_enum==constant_enum){
+			ioconstant->constant->GetParameterValue(pvalue);
+			return;
+		}
+	}
+
+	_error_("Could not find constant \""<<EnumToStringx(constant_enum) <<"\"");
+}
+/*}}}*/
+int   IoModel::NumIndependents(void){/*{{{*/
+
+	/*Initialize output*/
+	int num_independents = 0;
+
+	/*Process constants*/
+	for(vector<IoConstant*>::iterator iter=constants.begin();iter<constants.end();iter++){
+		if((*iter)->isindependent){
+			num_independents+= 1;
+		}
+	}
+
+	/*Process data*/
+	for(vector<IoData*>::iterator iter=data.begin();iter<data.end();iter++){
+		if((*iter)->isindependent){
+			num_independents+= (*iter)->M*(*iter)->N;
+		}
+	}
+
+	/*return*/
+	return num_independents;
+}
+/*}}}*/
+void  IoModel::FillIndependents(IssmDouble* xp){/*{{{*/
+
+	_assert_(xp);
+
+	/*Initialize local num ind*/
+	int local_num_ind = 0;
+
+	/*Process constants*/
+	for(vector<IoConstant*>::iterator iter=constants.begin();iter<constants.end();iter++){
+		if((*iter)->isindependent){
+			(*iter)->constant->GetParameterValue(&xp[local_num_ind]);
+			local_num_ind += 1;
+		}
+	}
+
+	/*Process data*/
+	for(vector<IoData*>::iterator iter=data.begin();iter<data.end();iter++){
+		if((*iter)->isindependent){
+			for(int i=0;i<(*iter)->M*(*iter)->N;i++){
+				xp[local_num_ind+i] = (*iter)->data[i];
+			}
+			local_num_ind += (*iter)->M*(*iter)->N;
+		}
+	}
+
+	_assert_(local_num_ind == this->NumIndependents());
+}
+/*}}}*/
+/*NEW NEW NEW NEW NEW NEW NEW NEW NEW NEW*/
+
 IoModel::IoModel(){/*{{{*/
+
 	this->fid=NULL;
 	this->solution_enum=-1;
-	this->data=NULL;
-	this->independents=NULL;
-	this->independent_objects=NULL;
-	this->constants=NULL;
 
 	this->my_elements=NULL;
@@ -65,8 +426,4 @@
 	this->solution_enum = solution_enum_in;
 
-	/*Initialize data: */
-	this->data=xNew<IssmDouble*>(MaximumNumberOfDefinitionsEnum);
-	for(int i=0;i<MaximumNumberOfDefinitionsEnum;i++) this->data[i]=NULL;
-
 	/*If we are running in AD mode, we need to start the trace and declare our independent variables now, 
 	 *and prevent them from being erased during successive calls to iomodel->FetchConstants, iomodel->FetchData and 
@@ -76,21 +433,19 @@
 
 	/*Initialize and read constants:*/
-	this->constants=new Parameters();
-	this->FetchConstants(); /*this routine goes through the input file, and fetches bools, ints, IssmDoubles and strings only, nothing memory intensive*/
-
-	/*Now some very specific piece of logic. We want to know whether we are going to carry out an autodiff run. There are 
-	 *several cases: either trace is true (which bypasses the isautodiff choice), or trace is false but isautodiff is on, in which 
-	 *case two cases are possible: are we running a control method or not? To make things simpler, we are going to pin down everyting 
-	 *on isautodiff. By the way, we could not do this before we had the constants dataset initialized! {{{*/
-	this->constants->FindParam(&autodiff,AutodiffIsautodiffEnum);
-	this->constants->FindParam(&iscontrol,InversionIscontrolEnum);
-
-	if(trace)autodiff=true;
+	this->FetchConstants(); /*this routine goes through the input file, and fetches bool, int, IssmDouble and string only, nothing memory intensive*/
+
+	/*Is this an autodiff run?*/
+	this->FindConstant(&autodiff,AutodiffIsautodiffEnum);
+	this->FindConstant(&iscontrol,InversionIscontrolEnum);
+	if(trace){
+		autodiff=true;
+	}
 	else{
-		if(autodiff && !iscontrol) autodiff=true;
-		else autodiff=false;
-	}
-	this->constants->SetParam(autodiff,AutodiffIsautodiffEnum);
-	/*}}}*/
+		if(autodiff && !iscontrol)
+		 autodiff=true;
+		else
+		 autodiff=false;
+	}
+	this->AddConstant(new IoConstant(autodiff,AutodiffIsautodiffEnum));
 
 	/*Initialize permanent data: */
@@ -119,30 +474,22 @@
 IoModel::~IoModel(){/*{{{*/
 
-	/*Delete independents*/
-	if(this->independents){
-		for(int i=0;i<MaximumNumberOfDefinitionsEnum;i++){
-			if(this->independents[i]){
-				IssmDouble* array=this->data[i];
-				xDelete<IssmDouble>(array);
-			}
-		}
-	}
-
-	/*checks in debugging mode*/
-	#if defined(_ISSM_DEBUG_) && !defined(_HAVE_ADOLC_)
-	if(this->data){
-		for(int i=0;i<MaximumNumberOfDefinitionsEnum;i++){
-			if(this->data[i]){
-				_printf0_("WARNING: previous pointer of " << EnumToStringx(i) << " has not been freed (DeleteData has not been called)\n");
-			}
-		}
-	}
-	#endif
-
-	if(this->constants) delete this->constants;
-
-	xDelete<IssmDouble*>(this->data);
-	xDelete<bool>(this->independents);
-	if(this->independent_objects)delete this->independent_objects;
+	/*Delete constants*/
+	vector<IoConstant*>::iterator iter1;
+	for(iter1=constants.begin();iter1<constants.end();iter1++){
+		delete *iter1;
+	}
+	this->constants.clear();
+
+	/*Delete data*/
+	vector<IoData*>::iterator iter2;
+	for(iter2=data.begin();iter2<data.end();iter2++){
+		#if defined(_ISSM_DEBUG_)
+		if(!(*iter2)->isindependent){
+			_printf0_("WARNING: IoData " << EnumToStringx((*iter2)->data_enum) << " has not been freed (DeleteData has not been called)\n");
+		}
+		#endif
+		delete *iter2;
+	}
+	this->data.clear();
 
 	xDelete<bool>(this->my_elements);
@@ -244,54 +591,32 @@
 }
 /*}}}*/
-void  IoModel::FindConstant(bool* poutput,int constant_enum){/*{{{*/
-
-	_assert_(constant_enum>=0);
-	_assert_(this->constants);
-
-	this->constants->FindParam(poutput,constant_enum);
-}
-/*}}}*/
-void  IoModel::FindConstant(int* poutput,int constant_enum){/*{{{*/
-
-	_assert_(constant_enum>=0);
-	_assert_(this->constants);
-
-	this->constants->FindParam(poutput,constant_enum);
-}
-/*}}}*/
-void  IoModel::FindConstant(IssmDouble* poutput,int constant_enum){/*{{{*/
-
-	_assert_(constant_enum>=0);
-	_assert_(this->constants);
-
-	this->constants->FindParam(poutput,constant_enum);
-}
-/*}}}*/
-void  IoModel::FindConstant(char** poutput,int constant_enum){/*{{{*/
-
-	_assert_(constant_enum>=0);
-	_assert_(this->constants);
-
-	this->constants->FindParam(poutput,constant_enum);
-}
-/*}}}*/
 Param* IoModel::CopyConstantObject(int constant_enum){/*{{{*/
 
-	_assert_(this->constants);
-
-	/*Find constant*/
-	//FIXME _error_("Needs to be FIXED");
-	Param* param=this->constants->FindParamObject(constant_enum);
-	if(!param) _error_("Constant " << EnumToStringx(constant_enum) << " not found in iomodel");
-
-	return xDynamicCast<Param*>(param->copy());
+	/*Intermediary*/
+	vector<IoConstant*>::iterator iter;
+
+	for(iter=constants.begin();iter<constants.end();iter++){
+		IoConstant* ioconstant=*iter;
+
+		if(ioconstant->data_enum==constant_enum){
+			return ioconstant->constant->copy();
+		}
+	}
+
+	_error_("Constant " << EnumToStringx(constant_enum) << " not found in iomodel");
+	return NULL;
 }
 /*}}}*/
 IssmDouble* IoModel::Data(int data_enum){/*{{{*/
 
-	_assert_(data_enum<MaximumNumberOfDefinitionsEnum);
-	_assert_(data_enum>=0);
-
-	return this->data[data_enum];
+	/*Intermediary*/
+	vector<IoData*>::iterator iter;
+
+	for(iter=data.begin();iter<data.end();iter++){
+		IoData* iodata=*iter;
+		if(iodata->data_enum==data_enum) return iodata->data;
+	}
+
+	return NULL;
 }
 /*}}}*/
@@ -338,7 +663,5 @@
 void  IoModel::DeclareIndependents(bool trace,IssmPDouble* X){/*{{{*/
 
-	int  i;
-	bool autodiff = false;
-	bool iscontrol = false;
+	bool autodiff,iscontrol;
 	int  num_independent_objects;
 	int  Xcount=0;
@@ -347,10 +670,5 @@
 	int *types = NULL;
 
-	int  dummy;
-
 	/*Initialize array detecting whether data[i] is an independent AD mode variable: */
-	this->independents=xNew<bool>(MaximumNumberOfDefinitionsEnum);
-	for(i=0;i<MaximumNumberOfDefinitionsEnum;i++) this->independents[i]=false;
-
 	this->FetchData(&autodiff,AutodiffIsautodiffEnum);
 	this->FetchData(&iscontrol,InversionIscontrolEnum);
@@ -359,23 +677,24 @@
 
 		#ifdef _HAVE_ADOLC_
-		/*build dataset made of independent objects:*/
-		this->independent_objects=new DataSet();
 		this->FetchData(&num_independent_objects,AutodiffNumIndependentObjectsEnum);
 		if(num_independent_objects){
-			this->FetchData(&names,&dummy,&dummy,AutodiffIndependentObjectNamesEnum);
-			this->FetchData(&types,&dummy,&dummy,AutodiffIndependentObjectTypesEnum);
+			this->FetchData(&names,NULL,NULL,AutodiffIndependentObjectNamesEnum);
+			this->FetchData(&types,NULL,NULL,AutodiffIndependentObjectTypesEnum);
 
 			/*create independent objects, and at the same time, fetch the corresponding independent variables, 
 			 *and declare them as such in ADOLC: */
-			for(i=0;i<num_independent_objects;i++){
-
-				IndependentObject* independent_object=NULL;
-				independent_object=new IndependentObject(names[i],types[i]);
-
-				/*add to independent_objects dataset:*/
-				this->independent_objects->AddObject(independent_object);
-
-				/*now go fetch the independent variable: */
-				independent_object->FetchIndependent(this,&Xcount,X); //supply the pointer to iomodel.
+			for(int i=0;i<num_independent_objects;i++){
+
+				if(types[i]==0){
+					/*Scalar*/
+					this->FetchIndependentConstant(&Xcount,X,names[i]);
+				}
+				else if(types[i]==1){
+					/* vector:*/
+					this->FetchIndependentData(&Xcount,X,names[i]);
+				}
+				else{
+					_error_("Independent cannot be of size " << types[i]);
+				}
 			}
 			xDelete<int>(names);
@@ -387,29 +706,46 @@
 		#endif
 	}
-	else this->independent_objects=NULL;
-
 }
 /*}}}*/
 void  IoModel::DeleteData(int num,...){/*{{{*/
 
+	/*Intermediaries*/
 	va_list ap;
 	int     dataenum;
-
-	/*Go through the entire list of enums and delete the corresponding data from the iomodel-data dataset: */
+	vector<IoData*>::iterator iter;
+
+	/*Go through the entire list of data and delete the corresponding data from the iomodel-data dataset: */
 	va_start(ap,num);
 	for(int i=0;i<num;i++){
 		dataenum=va_arg(ap,int);
-		_assert_(dataenum<MaximumNumberOfDefinitionsEnum);
-
-		/*do not erase independent variables for the AD mode computations!: */
-		if (!this->independents[dataenum]) xDelete<IssmDouble>(this->data[dataenum]);
+
+		for(iter=data.begin();iter<data.end();iter++){
+			IoData* iodata=*iter;
+			if(iodata->data_enum==dataenum && !iodata->isindependent){
+				delete *iter;
+				this->data.erase(iter);
+				break;
+			}
+		}
 	}
 	va_end(ap);
 } /*}}}*/
-void  IoModel::DeleteData(IssmDouble* vector, int dataenum){/*{{{*/
-
-	/*do not erase independent variables for the AD mode computations!: */
-	if(vector)if (!this->independents[dataenum]) xDelete<IssmDouble>(vector);
-
+void  IoModel::DeleteData(IssmDouble* vector_in, int dataenum){/*{{{*/
+
+	vector<IoData*>::iterator iter;
+
+	/*do not do anything if pointer is NULL*/
+	if(!vector_in) return;
+
+	/*do not delete if this is an independent variable*/
+	for(iter=data.begin();iter<data.end();iter++){
+		IoData* iodata=*iter;
+		if(iodata->data_enum==dataenum && iodata->isindependent){
+			return;
+		}
+	}
+
+	/*Go ahead and delete*/
+	xDelete<IssmDouble>(vector_in);
 } /*}}}*/
 void  IoModel::DeleteData(char*** pstringarray, int numstrings, int dataenum){/*{{{*/
@@ -429,6 +765,4 @@
 void  IoModel::FetchConstants(void){/*{{{*/
 
-	int my_rank;
-
 	/*record descriptions; */
 	int record_enum;
@@ -437,17 +771,16 @@
 
 	/*records: */
-	int  booleanint=0;
-	int  integer=0;
-	IssmPDouble pscalar=0;
-	IssmDouble scalar=0;
-	char* string=NULL;
-	int   string_size;
+	int          booleanint  = 0;
+	int          integer     = 0;
+	IssmPDouble  pscalar     = 0;
+	IssmDouble   scalar      = 0;
+	char        *string      = NULL;
+	int          string_size;
 
 	/*recover my_rank:*/
-	my_rank=IssmComm::GetRank();
+	int my_rank=IssmComm::GetRank();
 
 	/*Check that some fields have been allocated*/
 	_assert_(this->fid || my_rank);
-	_assert_(this->constants);
 
 	/*Go find in the binary file, the position of the data we want to fetch: */
@@ -487,5 +820,5 @@
 						/*create BoolParam: */
 						if(record_enum!=MaximumNumberOfDefinitionsEnum && record_enum!=MaximumNumberOfDefinitionsEnum+1)
-						 this->constants->AddObject(new BoolParam(record_enum,(bool)booleanint)); //cast to boolean
+						 this->AddConstant(new IoConstant((bool)booleanint,record_enum)); //cast to boolean
 
 						break;
@@ -496,23 +829,28 @@
 
 						/*create IntParam: */
-						this->constants->AddObject(new IntParam(record_enum,integer));
+						this->AddConstant(new IoConstant(integer,record_enum));
 
 						break;
 					case 3:
-						/*Read the scalar and broadcast it to other cpus. However, if this record has already 
-						 * been read in DeclareIndependents, we should not re-read it, but grab it from the iomodel->data 
-						 * slots: */
-						if(this->independents[record_enum]){
-							scalar=*(this->data[record_enum]);
-						}
-						else{
-							if(fread(&pscalar,sizeof(IssmPDouble),1,this->fid)!=1) _error_("could not read scalar ");
-							ISSM_MPI_Bcast(&pscalar,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); 
-							scalar=pscalar;
-						}
-
-						/*create DoubleParam: */
-						this->constants->AddObject(new DoubleParam(record_enum,scalar));
-
+						  {
+							/*IssmDouble, check whether it is already there (from "declare independents")*/
+							bool exists = false;
+							vector<IoConstant*>::iterator iter;
+							for(iter=constants.begin();iter<constants.end();iter++){
+								IoConstant* ioconstant=*iter;
+								if(ioconstant->data_enum==record_enum){
+									exists = true;
+									break;
+								}
+							}
+							if(!exists){
+								if(fread(&pscalar,sizeof(IssmPDouble),1,this->fid)!=1) _error_("could not read scalar ");
+								ISSM_MPI_Bcast(&pscalar,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); 
+								scalar=pscalar;
+
+								/*create DoubleParam: */
+								this->AddConstant(new IoConstant(scalar,record_enum));
+							}
+						  }
 						break;
 					case 4: 
@@ -535,36 +873,13 @@
 
 						/*Add string to parameters: */
-						this->constants->AddObject(new StringParam(record_enum,string));
+						this->AddConstant(new IoConstant(string,record_enum));
 
 						/*Free string*/
 						xDelete<char>(string);
-
 						break;
 					case 5: 
-							/*We are not interested in this record, too memory intensive. Skip it: */
-							/*skip: */
-							fseek(fid,-sizeof(int),SEEK_CUR); //backtrak 1 integer
-							fseek(fid,record_length,SEEK_CUR);
-							break;
 					case 6: 
-							/*We are not interested in this record, too memory intensive. Skip it: */
-							/*skip: */
-							fseek(fid,-sizeof(int),SEEK_CUR); //backtrak 1 integer
-							fseek(fid,record_length,SEEK_CUR);
-							break;
 					case 7: 
-							/*We are not interested in this record, too memory intensive. Skip it: */
-							/*skip: */
-							fseek(fid,-sizeof(int),SEEK_CUR); //backtrak 1 integer
-							fseek(fid,record_length,SEEK_CUR);
-							break;
-
 					case 8: 
-							/*We are not interested in this record, too memory intensive. Skip it: */
-							/*skip: */
-							fseek(fid,-sizeof(int),SEEK_CUR); //backtrak 1 integer
-							fseek(fid,record_length,SEEK_CUR);
-							break;
-
 					case 9: 
 							/*We are not interested in this record, too memory intensive. Skip it: */
@@ -573,8 +888,7 @@
 							fseek(fid,record_length,SEEK_CUR);
 							break;
-
 					default: 
 						_error_("unknown record type:" << record_code); 
-						break;;
+						break;
 				}
 			}
@@ -597,5 +911,5 @@
 					/*create BoolParam: */
 					if(record_enum!=MaximumNumberOfDefinitionsEnum && record_enum!=MaximumNumberOfDefinitionsEnum+1)
-					 this->constants->AddObject(new BoolParam(record_enum,(bool)booleanint)); //cast to a boolean
+					 this->AddConstant(new IoConstant((bool)booleanint,record_enum)); //cast to a boolean
 					break;
 
@@ -605,14 +919,27 @@
 
 					/*create IntParam: */
-					this->constants->AddObject(new IntParam(record_enum,integer));
+					this->AddConstant(new IoConstant(integer,record_enum));
 
 					break;
 				case 3:
 					/*scalar. get it from cpu 0 */
-					ISSM_MPI_Bcast(&pscalar,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
-					scalar=pscalar;
-					/*create DoubleParam: */
-					this->constants->AddObject(new DoubleParam(record_enum,scalar));
-
+					  {
+						/*IssmDouble, check whether it is already there (from "declare independents")*/
+						bool exists = false;
+						vector<IoConstant*>::iterator iter;
+						for(iter=constants.begin();iter<constants.end();iter++){
+							IoConstant* ioconstant=*iter;
+							if(ioconstant->data_enum==record_enum){
+								exists = true;
+								break;
+							}
+						}
+						if(!exists){
+							ISSM_MPI_Bcast(&pscalar,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
+							scalar=pscalar;
+							/*create DoubleParam: */
+							this->AddConstant(new IoConstant(scalar,record_enum));
+						}
+					  }
 					break;
 				case 4: 
@@ -630,5 +957,5 @@
 					}
 					/*Add string to parameters: */
-					this->constants->AddObject(new StringParam(record_enum,string));
+					this->AddConstant(new IoConstant(string,record_enum));
 
 					/*Free string*/
@@ -644,5 +971,5 @@
 				default: 
 					_error_("unknown record type:" << record_code); 
-					break;;
+					break;
 				}
 
@@ -653,6 +980,4 @@
 /*}}}*/
 void  IoModel::FetchData(bool* pboolean,int data_enum){/*{{{*/
-
-	int my_rank;
 
 	/*output: */
@@ -661,5 +986,5 @@
 
 	/*recover my_rank:*/
-	my_rank=IssmComm::GetRank();
+	int my_rank=IssmComm::GetRank();
 
 	/*Set file pointer to beginning of the data: */
@@ -682,6 +1007,4 @@
 void  IoModel::FetchData(int* pinteger,int data_enum){/*{{{*/
 
-	int my_rank;
-
 	/*output: */
 	int   integer;
@@ -689,5 +1012,5 @@
 
 	/*recover my_rank:*/
-	my_rank=IssmComm::GetRank();
+	int my_rank=IssmComm::GetRank();
 
 	/*Set file pointer to beginning of the data: */
@@ -709,6 +1032,4 @@
 void  IoModel::FetchData(IssmDouble* pscalar,int data_enum){/*{{{*/
 
-	int my_rank;
-
 	/*output: */
 	IssmPDouble   scalar;
@@ -716,5 +1037,5 @@
 
 	/*recover my_rank:*/
-	my_rank=IssmComm::GetRank();
+	int my_rank=IssmComm::GetRank();
 
 	/*Set file pointer to beginning of the data: */
@@ -735,6 +1056,4 @@
 /*}}}*/
 void  IoModel::FetchData(char** pstring,int data_enum){/*{{{*/
-
-	int my_rank;
 
 	/*output: */
@@ -744,5 +1063,5 @@
 
 	/*recover my_rank:*/
-	my_rank=IssmComm::GetRank();
+	int my_rank=IssmComm::GetRank();
 
 	/*Set file pointer to beginning of the data: */
@@ -781,6 +1100,4 @@
 /*}}}*/
 void  IoModel::FetchData(int** pmatrix,int* pM,int* pN,int data_enum){/*{{{*/
-
-	int my_rank;
 	int i,j;
 
@@ -792,5 +1109,5 @@
 
 	/*recover my_rank:*/
-	my_rank=IssmComm::GetRank();
+	int my_rank=IssmComm::GetRank();
 
 	/*Set file pointer to beginning of the data: */
@@ -850,13 +1167,23 @@
 void  IoModel::FetchData(IssmDouble** pmatrix,int* pM,int* pN,int data_enum){/*{{{*/
 
-	int my_rank;
-
+	/*First, look if has already been loaded (might be an independent variable)*/
+	vector<IoData*>::iterator iter;
+	for(iter=data.begin();iter<data.end();iter++){
+		IoData* iodata=*iter;
+		if(iodata->data_enum==data_enum){
+			*pmatrix=iodata->data;
+			if(pM) *pM=iodata->M;
+			if(pN) *pN=iodata->N;
+			return;
+		}
+	}
+	 
 	/*output: */
-	int M,N;
-	IssmPDouble* matrix=NULL;
-	int code=0;
+	int          M,N;
+	IssmPDouble *matrix = NULL;
+	int          code   = 0;
 
 	/*recover my_rank:*/
-	my_rank=IssmComm::GetRank();
+	int my_rank=IssmComm::GetRank();
 
 	/*Set file pointer to beginning of the data: */
@@ -868,5 +1195,5 @@
 	/*We have to read a matrix from disk. First read the dimensions of the matrix, then the whole matrix: */
 	/*numberofelements: */
-	if(my_rank==0){  
+	if(my_rank==0){
 		if(fread(&M,sizeof(int),1,fid)!=1) _error_("could not read number of rows for matrix ");
 	}
@@ -888,14 +1215,6 @@
 		ISSM_MPI_Bcast(matrix,M*N,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); 
 
-		_assert_(this->independents);
-		if(this->independents[data_enum]){
-			/*this data has already been checked out! So cancel all that we've done here, and return 
-			 * the data[data_enum] directly: */
-			*pmatrix=this->data[data_enum];
-		}
-		else{
-			*pmatrix=xNew<IssmDouble>(M*N);
-			for(int i=0;i<M*N;++i) (*pmatrix)[i]=matrix[i];
-		}
+		*pmatrix=xNew<IssmDouble>(M*N);
+		for(int i=0;i<M*N;++i) (*pmatrix)[i]=matrix[i];
 		xDelete<IssmPDouble>(matrix);
 	}
@@ -903,11 +1222,9 @@
 	  *pmatrix=NULL;
 	/*Assign output pointers: */
-	if (pM)*pM=M;
-	if (pN)*pN=N;
+	if(pM) *pM=M;
+	if(pN) *pN=N;
 }
 /*}}}*/
 void  IoModel::FetchData(char*** pstrings,int* pnumstrings,int data_enum){/*{{{*/
-
-	int my_rank;
 
 	int i;
@@ -923,5 +1240,5 @@
 
 	/*recover my_rank:*/
-	my_rank=IssmComm::GetRank();
+	int my_rank=IssmComm::GetRank();
 
 	/*Set file pointer to beginning of the data: */
@@ -975,7 +1292,4 @@
 
 	int i;
-
-	int my_rank;
-
 	/*output: */
 	IssmDouble** matrices=NULL;
@@ -990,5 +1304,5 @@
 
 	/*recover my_rank:*/
-	my_rank=IssmComm::GetRank();
+	int my_rank=IssmComm::GetRank();
 
 	/*Set file pointer to beginning of the data: */
@@ -1105,8 +1419,9 @@
 
 	va_list     ap;
-	int         dataenum;
+	int         dataenum,code,layout;
 	IssmDouble *matrix   = NULL;
 	int         M,N;
-	int         i;
+	bool        exists;
+	vector<IoData*>::iterator iter;
 
 	/*Go through the entire list of enums and fetch the corresponding data. Add it to the iomodel->data dataset. Everything
@@ -1114,31 +1429,31 @@
 
 	va_start(ap,num);
-	for(i=0; i<num; i++){
+	for(int i=0; i<num; i++){
 
 		dataenum=va_arg(ap, int);
-		_assert_(dataenum<MaximumNumberOfDefinitionsEnum); 
-		_assert_(dataenum>=0);
-
-		if (this->independents[dataenum]){
+		exists = false;
+
+		for(iter=data.begin();iter<data.end();iter++){
+			IoData* iodata=*iter;
+			if(iodata->data_enum==dataenum){
+				/*Already there, no need to fetch it*/
+				_assert_(iodata->isindependent);
+				exists = true;
+				break;
+			}
+		}
+
+		if(exists){
 			/*this data has already been checked out! Continue: */
 			continue;
 		}
-
-		/*Some checks in debugging mode*/
-		/*{{{*/
-		#ifdef _ISSM_DEBUG_
-		_assert_(dataenum<MaximumNumberOfDefinitionsEnum);
-		if(this->data[dataenum]){
-			_error_("Info: trying to fetch " << EnumToStringx(dataenum) << " but previous pointer has not been freed (DeleteData has not been called)");
-		}
-		#endif
-		/*}}}*/
-
-		/*Add to this->data: */
-		this->FetchData(&matrix,&M,&N,dataenum);
-		this->data[dataenum]=matrix;
+		else{
+			/*Add to this->data: */
+			this->SetFilePointerToData(&code,&layout,dataenum);
+			this->FetchData(&matrix,&M,&N,dataenum);
+			this->AddData(new IoData(matrix,code,layout,M,N,dataenum));
+		}
 	}
 	va_end(ap);
-
 }
 /*}}}*/
@@ -1373,14 +1688,13 @@
 				ISSM_MPI_Bcast(pmatrix,M*N,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); 
 
-				_assert_(this->independents);
-				if(this->independents[data_enum]){
-					/*this data has already been checked out! So cancel all that we've done here, and return 
-					 * the data[data_enum] directly: */
-					matrix=this->data[data_enum];
-				}
-				else{
+				//if(this->independents[data_enum]){ FIXME
+				//	/*this data has already been checked out! So cancel all that we've done here, and return 
+				//	 * the data[data_enum] directly: */
+				//	matrix=this->data[data_enum];
+				//}
+				//else{
 					matrix=xNew<IssmDouble>(M*N);
 					for (int i=0;i<M*N;++i) matrix[i]=pmatrix[i];
-				}
+				//}
 				xDelete<IssmPDouble>(pmatrix);
 			}
@@ -1479,14 +1793,13 @@
 				ISSM_MPI_Bcast(pmatrix,M*N,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); 
 
-				_assert_(this->independents);
-				if(this->independents[data_enum]){
-					/*this data has already been checked out! So cancel all that we've done here, and return 
-					 * the data[data_enum] directly: */
-					matrix=this->data[data_enum];
-					for (int i=0;i<M*N;++i) integer_matrix[i]=reCast<int>(matrix[i]);
-				}
-				else{
+				//if(this->independents[data_enum]){ FIXME
+				//	/*this data has already been checked out! So cancel all that we've done here, and return 
+				//	 * the data[data_enum] directly: */
+				//	matrix=this->data[data_enum];
+				//	for (int i=0;i<M*N;++i) integer_matrix[i]=reCast<int>(matrix[i]);
+				//}
+				//else{
 					for (int i=0;i<M*N;++i) integer_matrix[i]=pmatrix[i];
-				}
+				//}
 				xDelete<IssmPDouble>(pmatrix);
 			}
@@ -1525,4 +1838,18 @@
 void  IoModel::FetchDataToInput(Elements* elements,int vector_enum,IssmDouble default_value){/*{{{*/
 
+	/*First, look whether it is not already loaded in this->data*/
+	vector<IoData*>::iterator iter;
+	for(iter=data.begin();iter<data.end();iter++){
+		IoData* iodata=*iter;
+		if(iodata->data_enum==vector_enum){
+			_assert_(iodata->code==7);
+			for(int i=0;i<elements->Size();i++){
+				Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
+				element->InputCreate(iodata->data,this,iodata->M,iodata->N,iodata->layout,vector_enum,iodata->code);//we need i to index into elements.
+			}
+			return;
+		}
+	}
+
 	/*intermediary: */
 	int         code,vector_layout;
@@ -1544,9 +1871,22 @@
 	}
 
-	/*Free ressources. Pay attention to not freeing an AD mode independent variable though!:*/
-	if(!this->independents[vector_enum]) xDelete<IssmDouble>(doublearray);
+	/*Free ressources*/
+	xDelete<IssmDouble>(doublearray);
 }
 /*}}}*/
 void  IoModel::FetchDataToInput(Elements* elements,int vector_enum){/*{{{*/
+
+	/*First, look whether it is not already loaded in this->data*/
+	vector<IoData*>::iterator iter;
+	for(iter=data.begin();iter<data.end();iter++){
+		IoData* iodata=*iter;
+		if(iodata->data_enum==vector_enum){
+			for(int i=0;i<elements->Size();i++){
+				Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
+				element->InputCreate(iodata->data,this,iodata->M,iodata->N,iodata->layout,vector_enum,iodata->code);//we need i to index into elements.
+			}
+			return;
+		}
+	}
 
 	/*intermediary: */
@@ -1615,6 +1955,6 @@
 			break;
 	}
-	/*Free ressources. Pay attention to not freeing an AD mode independent variable though!:*/
-	if (!this->independents[vector_enum]) xDelete<IssmDouble>(doublearray);
+	/*Free ressources*/
+	xDelete<IssmDouble>(doublearray);
 	xDelete<char>(string);
 }
