/*! \file IoModel.cpp
 * \brief  file containing the methods that will help in processing the input data coming 
 * into ISSM, from Matlab, or through a binary file opened for reading.
 */

#ifdef HAVE_CONFIG_H
	#include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>

#include "./classes.h"
#include "../shared/io/io.h"
#include "../shared/shared.h"

IoModel::IoModel(){/*{{{*/
	this->fid=NULL;
	this->data=NULL;
	this->independents=NULL;
	this->independent_objects=NULL;
	this->constants=NULL;

	this->my_elements=NULL;
	this->my_vertices=NULL;

	this->domaintype=-1;
	this->domaindim=-1;
	this->meshelementtype=-1;
	this->numberofvertices=-1;
	this->numberofelements=-1;
	this->numberoffaces=-1;
	this->numberofedges=-1;
	this->facescols=-1;
	this->elements=NULL;
	this->faces=NULL;
	this->edges=NULL;
	this->elementtofaceconnectivity      =NULL;
	this->elementtoedgeconnectivity      =NULL;
	this->singlenodetoelementconnectivity=NULL;
	this->numbernodetoelementconnectivity=NULL;

	this->nodecounter=0;
	this->loadcounter=0;
	this->constraintcounter=0;
}
/*}}}*/
IoModel::IoModel(FILE* iomodel_handle){/*{{{*/

	/*First, keep track of the file handle: */
	this->fid=iomodel_handle;

	/*Check that Enums are Synchronized*/
	this->CheckEnumSync();

	/*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 
	 iomodel->DeleteData:*/
	this->StartTrace();
	this->DeclareIndependents();

	/*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*/

	/*Initialize permanent data: */
	this->my_elements = NULL;
	this->my_vertices = NULL;

	FetchData(&this->domaintype,DomainTypeEnum);
	FetchData(&this->domaindim,DomainDimensionEnum);
	FetchData(&this->meshelementtype,MeshElementtypeEnum);
	FetchData(&this->numberofvertices,MeshNumberofverticesEnum);
	FetchData(&this->numberofelements,MeshNumberofelementsEnum);
	FetchData(&this->elements,NULL,NULL,MeshElementsEnum);
	this->facescols                       = -1;
	this->faces                           = NULL;
	this->edges                           = NULL;
	this->elementtofaceconnectivity       = NULL;
	this->elementtoedgeconnectivity       = NULL;
	this->singlenodetoelementconnectivity = NULL;
	this->numbernodetoelementconnectivity = NULL;

	this->nodecounter=0;
	this->loadcounter=0;
	this->constraintcounter=0;
}
/*}}}*/
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*/
	#ifdef _ISSM_DEBUG_
	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;

	xDelete<bool>(this->my_elements);
	xDelete<int>(this->my_vertices);

	xDelete<int>(this->elements);
	xDelete<int>(this->faces);
	xDelete<int>(this->edges);
	xDelete<int>(this->elementtofaceconnectivity);
	xDelete<int>(this->elementtoedgeconnectivity);
	xDelete<int>(this->singlenodetoelementconnectivity);
	xDelete<int>(this->numbernodetoelementconnectivity);
}
/*}}}*/

void  IoModel::CheckEnumSync(void){/*{{{*/

	int record_enum,record_length;

	/*recover my_rank:*/
	int my_rank=IssmComm::GetRank();

	/*Check that some fields have been allocated*/
	_assert_(this->fid || my_rank);

	/*Go find in the binary file, the position of the data we want to fetch: */
	if(my_rank==0){ //cpu 0

		/*First set FILE* position to the beginning of the file: */
		fseek(this->fid,0,SEEK_SET);

		/*Get first Enum*/
		if(fread(&record_enum,sizeof(int),1,this->fid)==0){
			_error_("Marshalled file is empty");
		}
		else{
			if(record_enum!=MaximumNumberOfDefinitionsEnum){
				_printf0_("\n");
				_printf0_("=========================================================================\n");
				_printf0_(" Enums in marshalled file are not compatible with compiled code          \n");
				_printf0_("                                                                         \n");
				_printf0_("   * If you are running ISSM on your local machine:                      \n");
				_printf0_("     make sure that all the code is compiled and installed (make install)\n");
				_printf0_("   * If you are running ISSM on a remote cluster:                        \n");
				_printf0_("     make sure that you are using the same version of ISSM on your local \n");
				_printf0_("     machine and remote cluster (you might need to run svn update)       \n");
				_printf0_("   * If you are a developer and just added a new Enum:                   \n");
				_printf0_("     you might need to run ./Synchronize.sh in src/c/shared/Enum         \n");
				_printf0_("     and recompile                                                       \n");
				_printf0_("=========================================================================\n");
				_printf0_("\n");
				_error_("Enums not consistent (See error message above)");
			}
		}

		/*Get last enum*/
		bool found     = false;
		int  last_enum = 0;

		/*First set FILE* position to the beginning of the file: */
		fseek(fid,0,SEEK_SET);
		for(;;){
			/*Have we reached the end of file ?*/
			if(fread(&record_enum,sizeof(int),1,fid)==0){
				break;
			}

			/*Have we found the last Enum ?*/
			if(record_enum==MaximumNumberOfDefinitionsEnum+1){
				found = true;
				break;
			}

			/*Check that record_enum is an enum*/
			if(record_enum>=0 && record_enum<=MaximumNumberOfDefinitionsEnum){
				if(record_enum>0) last_enum = record_enum;
			}

			/*Go to next Enum*/
			if(fread(&record_length,sizeof(int),1,fid)!=1) _error_("Could not read record_length");
			fseek(fid,record_length,SEEK_CUR);
		}
		if(!found){
			_printf0_("\n");
			_printf0_("=========================================================================\n");
			_printf0_(" Marshalled file is corrupted                                            \n");
			_printf0_("                                                                         \n");
			_printf0_("   * If you are running an old model, send it to the ISSM developers     \n");
			_printf0_("     so that a check is added before marshall                            \n");
			_printf0_("   * Last Enum found: " << EnumToStringx(last_enum)<<"Enum ("<<last_enum<<")\n");
			_printf0_("     the corresponding model field has probably been marshalled          \n");
			_printf0_("     incorrectly                                                         \n");
			_printf0_("                                                                         \n");
			_printf0_("=========================================================================\n");
			_printf0_("\n");
			_error_("Binary file corrupted (See error message above)");
		}
	}
}
/*}}}*/
void IoModel::Constant(bool* poutput,int constant_enum){/*{{{*/

	_assert_(constant_enum>=0);
	_assert_(this->constants);

	this->constants->FindParam(poutput,constant_enum);
}
/*}}}*/
void IoModel::Constant(int* poutput,int constant_enum){/*{{{*/

	_assert_(constant_enum>=0);
	_assert_(this->constants);

	this->constants->FindParam(poutput,constant_enum);
}
/*}}}*/
void IoModel::Constant(IssmDouble* poutput,int constant_enum){/*{{{*/

	_assert_(constant_enum>=0);
	_assert_(this->constants);

	this->constants->FindParam(poutput,constant_enum);
}
/*}}}*/
void IoModel::Constant(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*/
	Param* param=dynamic_cast<Param*>(this->constants->FindParamObject(constant_enum));
	if(!param) _error_("Constant " << EnumToStringx(constant_enum) << " not found in iomodel");

	return dynamic_cast<Param*>(param->copy());
}
/*}}}*/
IssmDouble* IoModel::Data(int data_enum){/*{{{*/

	_assert_(data_enum<MaximumNumberOfDefinitionsEnum);
	_assert_(data_enum>=0);

	return this->data[data_enum];
}
/*}}}*/
void IoModel::StartTrace(void){/*{{{*/

	bool autodiff = false;
	bool keep=false;
	IssmDouble gcTriggerRatio;
	IssmDouble gcTriggerMaxSize;
	IssmDouble obufsize;
	IssmDouble lbufsize;
	IssmDouble cbufsize;
	IssmDouble tbufsize;

	int my_rank=IssmComm::GetRank();

	this->FetchData(&autodiff,AutodiffIsautodiffEnum);
	if(autodiff){

		#ifdef _HAVE_ADOLC_
		/*Retrieve parameters: */
		this->FetchData(&keep,AutodiffKeepEnum);
		int keepTaylors=keep?1:0;
		this->FetchData(&gcTriggerRatio,AutodiffGcTriggerRatioEnum);
		this->FetchData(&gcTriggerMaxSize,AutodiffGcTriggerMaxSizeEnum);
		this->FetchData(&obufsize,AutodiffObufsizeEnum);
		this->FetchData(&lbufsize,AutodiffLbufsizeEnum);
		this->FetchData(&cbufsize,AutodiffCbufsizeEnum);
		this->FetchData(&tbufsize,AutodiffTbufsizeEnum);

		/*Set garbage collection parameters: */
		setStoreManagerControl(reCast<IssmPDouble>(gcTriggerRatio),reCast<size_t>(gcTriggerMaxSize));

		/*Start trace: */
		int skipFileDeletion=1;
		trace_on(my_rank,keepTaylors,reCast<size_t>(obufsize),reCast<size_t>(lbufsize),reCast<size_t>(cbufsize),reCast<size_t>(tbufsize),skipFileDeletion);
		#endif
	}

}
/*}}}*/
void IoModel::DeclareIndependents(void){/*{{{*/

	int  i;
	bool autodiff = false;
	int  num_independent_objects;

	int *names = NULL;
	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);
	if(autodiff){

		#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);

			/*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); //supply the pointer to iomodel.
			}
			xDelete<int>(names);
			xDelete<int>(types);
		}
		#else
		/*if we asked for AD computations, we have a problem!: */
		_error_("Cannot carry out AD mode computations without support of ADOLC compiled in!");
		#endif
	}
	else this->independent_objects=NULL;

}
/*}}}*/
void  IoModel::DeleteData(int num,...){/*{{{*/

	va_list ap;
	int     dataenum;

	/*Go through the entire list of enums 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]);
	}
	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(char*** pstringarray, int numstrings, int dataenum){/*{{{*/

	int i;
	char** stringarray=*pstringarray;
	
	if (numstrings){
		for (i=0;i<numstrings;i++){
			char* string=stringarray[i]; 
			xDelete<char>(string);
		}
		xDelete<char*>(stringarray);
	}
	*pstringarray=NULL;
} /*}}}*/
void  IoModel::FetchConstants(void){/*{{{*/

	int my_rank;

	/*record descriptions; */
	int record_enum;
	int record_length;
	int record_code; //1 to 7 number

	/*records: */
	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();

	/*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: */
	if(my_rank==0){ //cpu 0{{{

		/*First set FILE* position to the beginning of the file: */
		fseek(this->fid,0,SEEK_SET);

		/*Now march through file looking for the correct data identifiers (bool,int,IssmDouble or string): */
		for(;;){
			if(fread(&record_enum,sizeof(int),1,this->fid)==0){

				/*Ok, we have reached the end of the file. break: */
				record_code=0; //0 means bailout
				ISSM_MPI_Bcast(&record_code,1,ISSM_MPI_INT,0,IssmComm::GetComm());  /*tell others cpus we are bailing: */
				break;
			}
			else{

				/* Read the record length and the data type code: */
				if(fread(&record_length,sizeof(int),1,this->fid)!=1) _error_("Cound not read record_length");
				if(fread(&record_code  ,sizeof(int),1,this->fid)!=1) _error_("Cound not read record_code");

				/*Tell other cpus what we are doing: */
				ISSM_MPI_Bcast(&record_code,1,ISSM_MPI_INT,0,IssmComm::GetComm());  /*tell other cpus what we are going to do: */

				/*Tell other cpus the name of the data, then branch according to the data type: */
				ISSM_MPI_Bcast(&record_enum,1,ISSM_MPI_INT,0,IssmComm::GetComm());  
				ISSM_MPI_Bcast(&record_length,1,ISSM_MPI_INT,0,IssmComm::GetComm());  

				switch(record_code){
					case 1: 
						/*Read the boolean and broadcast it to other cpus:*/
						if(fread(&booleanint,sizeof(int),1,this->fid)!=1) _error_("could not read boolean ");
						ISSM_MPI_Bcast(&booleanint,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

						/*create BoolParam: */
						this->constants->AddObject(new BoolParam(record_enum,(bool)booleanint)); //cast to boolean

						break;
					case 2:
						/*Read the integer and broadcast it to other cpus:*/
						if(fread(&integer,sizeof(int),1,this->fid)!=1) _error_("could not read integer ");
						ISSM_MPI_Bcast(&integer,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

						/*create IntParam: */
						this->constants->AddObject(new IntParam(record_enum,integer));

						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));

						break;
					case 4: 
						/*We have to read a string from disk. First read the dimensions of the string, then the string: */
						if(fread(&string_size,sizeof(int),1,this->fid)!=1) _error_("could not read length of string ");
						ISSM_MPI_Bcast(&string_size,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

						if(string_size){
							string=xNew<char>(string_size+1);
							string[string_size]='\0';

							/*Read string, then broadcast: */
							if(fread(string,string_size*sizeof(char),1,this->fid)!=1)_error_(" could not read string ");
							ISSM_MPI_Bcast(string,string_size,ISSM_MPI_CHAR,0,IssmComm::GetComm()); 
						}
						else{
							string=xNew<char>(1);
							string[0]='\0';
						}

						/*Add string to parameters: */
						this->constants->AddObject(new StringParam(record_enum,string));

						/*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: */
							/*skip: */
							fseek(fid,-sizeof(int),SEEK_CUR); //backtrak 1 integer
							fseek(fid,record_length,SEEK_CUR);
							break;

					default: 
						_error_("unknown record type:" << record_code); 
						break;;
				}
			}
		}
	} //}}}
	else{ //cpu ~0 {{{
		for(;;){ //wait on cpu 0
			ISSM_MPI_Bcast(&record_code,1,ISSM_MPI_INT,0,IssmComm::GetComm());  /*get from cpu 0 what we are going to do: */
			if(record_code==0){
				break; //we are done, break from the loop
			}
			else{
				ISSM_MPI_Bcast(&record_enum,1,ISSM_MPI_INT,0,IssmComm::GetComm());   //get from cpu 0 name of the data
				ISSM_MPI_Bcast(&record_length,1,ISSM_MPI_INT,0,IssmComm::GetComm());  
				switch(record_code){
				case 1: 
					/*boolean. get it from cpu 0 */
					ISSM_MPI_Bcast(&booleanint,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

					/*create BoolParam: */
					this->constants->AddObject(new BoolParam(record_enum,(bool)booleanint)); //cast to a boolean
					break;

				case 2:
					/*integer. get it from cpu 0 */
					ISSM_MPI_Bcast(&integer,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

					/*create IntParam: */
					this->constants->AddObject(new IntParam(record_enum,integer));

					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));

					break;
				case 4: 
					ISSM_MPI_Bcast(&string_size,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 
					if(string_size){
						string=xNew<char>((string_size+1));
						string[string_size]='\0';

						/*Read string from cpu 0: */
						ISSM_MPI_Bcast(string,string_size,ISSM_MPI_CHAR,0,IssmComm::GetComm()); 
					}
					else{
						string=xNew<char>(1);
						string[0]='\0';
					}
					/*Add string to parameters: */
					this->constants->AddObject(new StringParam(record_enum,string));

					/*Free string*/
					xDelete<char>(string);

					break;
				case 5: break; //do nothing. not interested in this type of data, which is memory intensive.
				case 6: break; //do nothing. not interested in this type of data, which is memory intensive.
				case 7: break; //do nothing. not interested in this type of data, which is memory intensive.
				case 8: break; //do nothing. not interested in this type of data, which is memory intensive.
				case 9: break; //do nothing. not interested in this type of data, which is memory intensive.

				default: 
					_error_("unknown record type:" << record_code); 
					break;;
				}

			}
		}
	} //}}}
}
/*}}}*/
void  IoModel::FetchData(bool* pboolean,int data_enum){/*{{{*/

	int my_rank;

	/*output: */
	int   booleanint;
	int   code;

	/*recover my_rank:*/
	my_rank=IssmComm::GetRank();

	/*Set file pointer to beginning of the data: */
	fid=this->SetFilePointerToData(&code,NULL,data_enum);

	if(code!=1)_error_("expecting a boolean for enum " << EnumToStringx(data_enum));

	/*We have to read a boolean from disk. */
	if(my_rank==0){  
		if(fread(&booleanint,sizeof(int),1,fid)!=1) _error_("could not read boolean ");
	}
	ISSM_MPI_Bcast(&booleanint,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

	/*cast to bool: */
	/*Assign output pointers: */
	*pboolean=(bool)booleanint;

}
/*}}}*/
void  IoModel::FetchData(int* pinteger,int data_enum){/*{{{*/

	int my_rank;

	/*output: */
	int   integer;
	int   code;

	/*recover my_rank:*/
	my_rank=IssmComm::GetRank();

	/*Set file pointer to beginning of the data: */
	fid=this->SetFilePointerToData(&code,NULL,data_enum);

	if(code!=2)_error_("expecting an integer for enum " << EnumToStringx(data_enum));

	/*We have to read a integer from disk. First read the dimensions of the integer, then the integer: */
	if(my_rank==0){  
		if(fread(&integer,sizeof(int),1,fid)!=1) _error_("could not read integer ");
	}

	ISSM_MPI_Bcast(&integer,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

	/*Assign output pointers: */
	*pinteger=integer;
}
/*}}}*/
void  IoModel::FetchData(IssmDouble* pscalar,int data_enum){/*{{{*/

	int my_rank;

	/*output: */
	IssmPDouble   scalar;
	int      code;

	/*recover my_rank:*/
	my_rank=IssmComm::GetRank();

	/*Set file pointer to beginning of the data: */
	fid=this->SetFilePointerToData(&code,NULL,data_enum);

	if(code!=3)_error_("expecting a IssmDouble for enum " << EnumToStringx(data_enum));

	/*We have to read a scalar from disk. First read the dimensions of the scalar, then the scalar: */
	if(my_rank==0){
		if(fread(&scalar,sizeof(IssmPDouble),1,fid)!=1)_error_("could not read scalar ");
	}
	ISSM_MPI_Bcast(&scalar,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); 

	/*Assign output pointers: */
	*pscalar=scalar;

}
/*}}}*/
void  IoModel::FetchData(char** pstring,int data_enum){/*{{{*/

	int my_rank;

	/*output: */
	char* string=NULL;
	int   string_size;
	int code=0;

	/*recover my_rank:*/
	my_rank=IssmComm::GetRank();

	/*Set file pointer to beginning of the data: */
	fid=this->SetFilePointerToData(&code,NULL,data_enum);

	if(code!=4)_error_("expecting a string for enum " << EnumToStringx(data_enum));

	/*Now fetch: */

	/*We have to read a string from disk. First read the dimensions of the string, then the string: */
	if(my_rank==0){  
		if(fread(&string_size,sizeof(int),1,fid)!=1) _error_("could not read length of string ");
	}

	ISSM_MPI_Bcast(&string_size,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

	/*Now allocate string: */
	if(string_size){
		string=xNew<char>((string_size+1));
		string[string_size]='\0';

		/*Read string on node 0, then broadcast: */
		if(my_rank==0){  
			if(fread(string,string_size*sizeof(char),1,fid)!=1)_error_(" could not read string ");
		}
		ISSM_MPI_Bcast(string,string_size,ISSM_MPI_CHAR,0,IssmComm::GetComm()); 
	}
	else{
		string=xNew<char>(1);
		string[0]='\0';
	}

	/*Assign output pointers: */
	*pstring=string;
}
/*}}}*/
void  IoModel::FetchData(int** pmatrix,int* pM,int* pN,int data_enum){/*{{{*/

	int my_rank;
	int i,j;

	/*output: */
	int M,N;
	IssmPDouble* matrix=NULL;
	int*    integer_matrix=NULL;
	int code=0;

	/*recover my_rank:*/
	my_rank=IssmComm::GetRank();

	/*Set file pointer to beginning of the data: */
	fid=this->SetFilePointerToData(&code,NULL,data_enum);

	if(code!=5 && code!=6 && code!=7)_error_("expecting a IssmDouble, integer or boolean matrix for enum " << EnumToStringx(data_enum)<<" (Code is "<<code<<")");

	/*Now fetch: */

	/*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){
		matrix=xNew<IssmPDouble>(M*N);

		/*Read matrix on node 0, then broadcast: */
		if(my_rank==0){  
			if(fread(matrix,M*N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read matrix ");
		}

		ISSM_MPI_Bcast(matrix,M*N,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); 
	}

	/*Now cast to integer: */
	if(M*N){
		integer_matrix=xNew<int>(M*N);
		for (i=0;i<M;i++){
			for (j=0;j<N;j++){
				integer_matrix[i*N+j]=(int)matrix[i*N+j];
			}
		}
	}
	else{
		integer_matrix=NULL;
	}
	/*Free ressources:*/
	xDelete<IssmPDouble>(matrix);

	/*Assign output pointers: */
	*pmatrix=integer_matrix;
	if (pM)*pM=M;
	if (pN)*pN=N;

}
/*}}}*/
void  IoModel::FetchData(IssmDouble** pmatrix,int* pM,int* pN,int data_enum){/*{{{*/

	int my_rank;

	/*output: */
	int M,N;
	IssmPDouble* matrix=NULL;
	int code=0;

	/*recover my_rank:*/
	my_rank=IssmComm::GetRank();

	/*Set file pointer to beginning of the data: */
	fid=this->SetFilePointerToData(&code,NULL,data_enum);
	if(code!=5 && code!=6 && code!=7)_error_("expecting a IssmDouble, integer or boolean matrix for enum " << EnumToStringx(data_enum)<<" (Code is "<<code<<")");

	/*Now fetch: */

	/*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){
		matrix=xNew<IssmPDouble>(M*N);

		/*Read matrix on node 0, then broadcast: */
		if(my_rank==0){  
			if(fread(matrix,M*N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read matrix ");
		}
		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];
		}
		xDelete<IssmPDouble>(matrix);
	}
	else
	  *pmatrix=NULL;
	/*Assign output pointers: */
	if (pM)*pM=M;
	if (pN)*pN=N;
}
/*}}}*/
void  IoModel::FetchData(char*** pstrings,int* pnumstrings,int data_enum){/*{{{*/

	int my_rank;

	int i;

	/*output: */
	int   numstrings=0;
	char** strings=NULL;

	/*intermediary: */
	char* string=NULL;
	int   string_size;
	int   code;

	/*recover my_rank:*/
	my_rank=IssmComm::GetRank();

	/*Set file pointer to beginning of the data: */
	fid=this->SetFilePointerToData(&code,NULL,data_enum);

	if(code!=9)_error_("expecting a string array for enum " << EnumToStringx(data_enum));

	/*We have to read a bunch of strings from disk. First read the number of strings, and allocate: */
	if(my_rank==0){  
		if(fread(&numstrings,sizeof(int),1,fid)!=1) _error_("could not read length of string array");
	}
	ISSM_MPI_Bcast(&numstrings,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

	/*Now allocate string array: */
	if(numstrings){
		strings=xNew<char*>(numstrings);
		for(i=0;i<numstrings;i++)strings[i]=NULL;

		/*Go through strings, and read: */
		for(i=0;i<numstrings;i++){

			if(my_rank==0){  
				if(fread(&string_size,sizeof(int),1,fid)!=1) _error_("could not read length of string ");
			}
			ISSM_MPI_Bcast(&string_size,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 
			if(string_size){
				string=xNew<char>((string_size+1));
				string[string_size]='\0';

				/*Read string on node 0, then broadcast: */
				if(my_rank==0){  
					if(fread(string,string_size*sizeof(char),1,fid)!=1)_error_(" could not read string ");
				}
				ISSM_MPI_Bcast(string,string_size,ISSM_MPI_CHAR,0,IssmComm::GetComm()); 
			}
			else{
				string=xNew<char>(1);
				string[0]='\0';
			}

			strings[i]=string;
		}
	}

	/*Assign output pointers: */
	*pstrings=strings;
	*pnumstrings=numstrings;
}
/*}}}*/
void  IoModel::FetchData(IssmDouble*** pmatrices,int** pmdims,int** pndims, int* pnumrecords,int data_enum){/*{{{*/

	int i;

	int my_rank;

	/*output: */
	IssmDouble** matrices=NULL;
	int*     mdims=NULL;
	int*     ndims=NULL;
	int      numrecords=0;

	/*intermediary: */
	int     M, N;
	IssmPDouble *matrix = NULL;
	int     code;

	/*recover my_rank:*/
	my_rank=IssmComm::GetRank();

	/*Set file pointer to beginning of the data: */
	fid=this->SetFilePointerToData(&code,NULL,data_enum);
	if(code!=8)_error_("expecting a IssmDouble mat array for enum " << EnumToStringx(data_enum));

	/*Now fetch: */
	if(my_rank==0){  
		if(fread(&numrecords,sizeof(int),1,fid)!=1) _error_("could not read number of records in matrix array ");
	}
	ISSM_MPI_Bcast(&numrecords,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

	if(numrecords){

		/*Allocate matrices :*/
		matrices=xNew<IssmDouble*>(numrecords);
		mdims=xNew<int>(numrecords);
		ndims=xNew<int>(numrecords);

		for(i=0;i<numrecords;i++){
			matrices[i]=NULL;
			mdims[i]=0;
			ndims[i]=0;
		}

		/*Loop through records and fetch matrix: */
		for(i=0;i<numrecords;i++){

			if(my_rank==0){  
				if(fread(&M,sizeof(int),1,fid)!=1) _error_("could not read number of rows in " << i << "th matrix of matrix array");
			}
			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 in " << i << "th matrix of matrix array");
			}
			ISSM_MPI_Bcast(&N,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

			/*Now allocate matrix: */
			if(M*N){
				matrix=xNew<IssmPDouble>(M*N);

				/*Read matrix on node 0, then broadcast: */
				if(my_rank==0){  
					if(fread(matrix,M*N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read matrix ");
				}

				ISSM_MPI_Bcast(matrix,M*N,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); 
				matrices[i]=xNew<IssmDouble>(M*N);
				for (int j=0;j<M*N;++j) {matrices[i][j]=matrix[j];}
				xDelete<IssmPDouble>(matrix);
			}
			else
			  matrices[i]=NULL;
			/*Assign: */
			mdims[i]=M;
			ndims[i]=N;
		}
	}

	/*Assign output pointers: */
	*pmatrices=matrices;
	*pmdims=mdims;
	*pndims=ndims;
	*pnumrecords=numrecords;
}
/*}}}*/
void  IoModel::FetchData(Option** poption,int index){/*{{{*/

	/*output: */
	int     code;
	char   *name        = NULL;

	/*First get option name*/
	this->FetchData(&name,index);

	/*Get option value*/
	fid=this->SetFilePointerToData(&code,NULL,index+1);
	switch(code){
		case 3: {//IssmDouble
					GenericOption<IssmDouble>* option;
					IssmDouble value;
					FetchData(&value,index+1);
					option = new GenericOption<IssmDouble>();
					option->value = value;
					option->name  = name;
					option->numel = 1;
					option->ndims = 1;
					option->size  = NULL;
					/*Assign output pointers: */
					*poption=option;
					break;
				}
		case 4: {//char
					GenericOption<char*>* option;
					char* value = NULL;
					FetchData(&value,index+1);
					option = new GenericOption<char*>();
					option->value = value;
					option->name  = name;
					option->numel = 1;
					option->ndims = 1;
					option->size  = NULL;
					*poption=option;
					break;
				}
		default:
			  _error_("Option of format " << code << " not supported yet");
	}

}
/*}}}*/
void  IoModel::FetchData(int num,...){/*{{{*/

	va_list     ap;
	int         dataenum;
	IssmDouble *matrix   = NULL;
	int         M,N;
	int         i;

	/*Go through the entire list of enums and fetch the corresponding data. Add it to the iomodel->data dataset. Everything
	 *we fetch is a IssmDouble* : */

	va_start(ap,num);
	for(i=0; i<num; i++){

		dataenum=va_arg(ap, int);
		_assert_(dataenum<MaximumNumberOfDefinitionsEnum); 
		_assert_(dataenum>=0);

		if (this->independents[dataenum]){
			/*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;
	}
	va_end(ap);

}
/*}}}*/
void  IoModel::FetchMultipleData(char*** pstrings,int* pnumstrings,int data_enum){/*{{{*/

	int  num_instances;

	/*output: */
	int    numstrings = 0;
	char **strings    = NULL;

	/*intermediary: */
	char   *string         = NULL;
	int     string_size;
	int    *codes          = NULL;
	int    *code           = NULL;
	fpos_t *file_positions = NULL;

	/*recover my_rank:*/
	int my_rank=IssmComm::GetRank();

	/*Get file pointers to beginning of the data (multiple instances of it): */
	file_positions=this->SetFilePointersToData(&codes,NULL,&num_instances,data_enum);

	if(num_instances){
		strings=xNew<char*>(num_instances);

		for(int i=0;i<num_instances;i++){

			if(my_rank==0){
				/*check we are indeed finding a string, not something else: */
				if(codes[i]!=4)_error_("expecting a string for enum " << EnumToStringx(data_enum));
		
				/*We have to read a string from disk. First read the dimensions of the string, then the string: */
				fsetpos(fid,file_positions+i);
				if(fread(&string_size,sizeof(int),1,fid)!=1) _error_("could not read length of string ");
			}

			ISSM_MPI_Bcast(&string_size,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

			/*Now allocate string: */
			if(string_size){
				string=xNew<char>((string_size+1));
				string[string_size]='\0';

				/*Read string on node 0, then broadcast: */
				if(my_rank==0){
					if(fread(string,string_size*sizeof(char),1,fid)!=1)_error_(" could not read string ");
				}
				ISSM_MPI_Bcast(string,string_size,ISSM_MPI_CHAR,0,IssmComm::GetComm()); 
			}
			else{
				string=xNew<char>(1);
				string[0]='\0';
			}
			strings[i]=string;
		}
	}
	/*Free ressources:*/
	xDelete<int>(codes);
	xDelete<fpos_t>(file_positions);
	
	/*Assign output pointers: */
	*pstrings=strings;
	*pnumstrings=num_instances;
}
/*}}}*/
void  IoModel::FetchMultipleData(int** pvector, int* pnum_instances,int data_enum){/*{{{*/

	int     num_instances;
	fpos_t* file_positions=NULL;

	/*output: */
	int* vector=NULL;

	/*intermediary: */
	int          integer;
	int         *codes   = NULL;
	int          code;

	/*recover my_rank:*/
	int my_rank=IssmComm::GetRank();
	
	/*Get file pointers to beginning of the data (multiple instances of it): */
	file_positions=this->SetFilePointersToData(&codes,NULL,&num_instances,data_enum);

	if(num_instances){

		/*Allocate vector :*/
		vector=xNew<int>(num_instances);

		for(int i=0;i<num_instances;i++){

			if(my_rank==0){
				code=codes[i];

				if(code!=2)_error_("expecting an integer for enum " << EnumToStringx(data_enum));
				
				/*We have to read a integer from disk. First read the dimensions of the integer, then the integer: */
				fsetpos(fid,file_positions+i);
				if(my_rank==0){  
					if(fread(&integer,sizeof(int),1,fid)!=1) _error_("could not read integer ");
				}
			}
			ISSM_MPI_Bcast(&integer,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

			/*Assign: */
			vector[i]=integer;
		}
	}
			
	/*Free ressources:*/
	xDelete<fpos_t>(file_positions);
	xDelete<int>(codes);

	/*Assign output pointers: */
	*pvector=vector;
	*pnum_instances=num_instances;
}
/*}}}*/
void  IoModel::FetchMultipleData(IssmDouble*** pmatrices,int** pmdims,int** pndims, int* pnumrecords,int data_enum){/*{{{*/

	int     num_instances;
	fpos_t* file_positions=NULL;

	/*output: */
	IssmDouble **matrices = NULL;
	int         *mdims    = NULL;
	int         *ndims    = NULL;

	/*intermediary: */
	int          M, N;
	IssmPDouble *pmatrix = NULL;
	IssmDouble  *matrix  = NULL;
	int         *codes   = NULL;
	int          code;

	/*recover my_rank:*/
	int my_rank=IssmComm::GetRank();
	
	/*Get file pointers to beginning of the data (multiple instances of it): */
	file_positions=this->SetFilePointersToData(&codes,NULL,&num_instances,data_enum);

	if(num_instances){

		/*Allocate matrices :*/
		matrices=xNew<IssmDouble*>(num_instances);
		mdims=xNew<int>(num_instances);
		ndims=xNew<int>(num_instances);

		for(int i=0;i<num_instances;i++){

			if(my_rank==0){
				code=codes[i];

				if((code!=5) && (code!=6) && (code!=7))_error_("expecting a IssmDouble, integer or boolean matrix for enum " << EnumToStringx(data_enum));

				/*We have to read a matrix from disk. First read the dimensions of the matrix, then the whole matrix: */
				/*numberofelements: */
				fsetpos(fid,file_positions+i);
				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){
				pmatrix=xNew<IssmPDouble>(M*N);

				/*Read matrix on node 0, then broadcast: */
				if(my_rank==0){  
					if(fread(pmatrix,M*N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read matrix ");
				}
				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{
					matrix=xNew<IssmDouble>(M*N);
					for (int i=0;i<M*N;++i) matrix[i]=pmatrix[i];
				}
				xDelete<IssmPDouble>(pmatrix);
			}
			else
				matrix=NULL;
			
			
			/*Assign: */
			mdims[i]=M;
			matrices[i]=matrix;
			ndims[i]=N;
		}
	}
			
	/*Free ressources:*/
	xDelete<fpos_t>(file_positions);
	xDelete<int>(codes);

	/*Assign output pointers: */
	*pmatrices=matrices;
	if(pmdims){
		*pmdims=mdims;
	}
	else{
		xDelete<int>(mdims);
	}
	if(pndims){
		*pndims=ndims;
	}
	else{
		xDelete<int>(ndims);
	}
	*pnumrecords=num_instances;
}
/*}}}*/
void  IoModel::FetchMultipleData(int*** pmatrices,int** pmdims,int** pndims, int* pnumrecords,int data_enum){/*{{{*/

	int     num_instances;
	fpos_t* file_positions=NULL;

	/*output: */
	int        **matrices = NULL;
	int         *mdims    = NULL;
	int         *ndims    = NULL;

	/*intermediary: */
	int          M, N;
	IssmPDouble *pmatrix = NULL;
	IssmDouble  *matrix  = NULL;
	int         *integer_matrix=NULL;
	int         *codes   = NULL;
	int          code;

	/*recover my_rank:*/
	int my_rank=IssmComm::GetRank();
	
	/*Get file pointers to beginning of the data (multiple instances of it): */
	file_positions=this->SetFilePointersToData(&codes,NULL,&num_instances,data_enum);

	if(num_instances){

		/*Allocate matrices :*/
		matrices=xNew<int*>(num_instances);
		mdims=xNew<int>(num_instances);
		ndims=xNew<int>(num_instances);

		for(int i=0;i<num_instances;i++){

			if(my_rank==0){
				code=codes[i];

				if((code!=5) && (code!=6) && (code!=7))_error_("expecting a IssmDouble, integer or boolean matrix for enum " << EnumToStringx(data_enum));

				/*We have to read a matrix from disk. First read the dimensions of the matrix, then the whole matrix: */
				/*numberofelements: */
				fsetpos(fid,file_positions+i);
				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){
				pmatrix=xNew<IssmPDouble>(M*N);
				integer_matrix=xNew<int>(M*N);

				/*Read matrix on node 0, then broadcast: */
				if(my_rank==0){  
					if(fread(pmatrix,M*N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read matrix ");
				}
				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{
					for (int i=0;i<M*N;++i) integer_matrix[i]=pmatrix[i];
				}
				xDelete<IssmPDouble>(pmatrix);
			}
			else
				integer_matrix=NULL;
			
			
			/*Assign: */
			mdims[i]=M;
			matrices[i]=integer_matrix;
			ndims[i]=N;
		}
	}
			
	/*Free ressources:*/
	xDelete<fpos_t>(file_positions);
	xDelete<int>(codes);

	/*Assign output pointers: */
	*pmatrices=matrices;
	if(pmdims){
		*pmdims=mdims;
	}
	else{
		xDelete<int>(mdims);
	}
	if(pndims){
		*pndims=ndims;
	}
	else{
		xDelete<int>(ndims);
	}
	*pnumrecords=num_instances;
}
/*}}}*/
void IoModel::FetchDataToInput(Elements* elements,int vector_enum,IssmDouble default_value){/*{{{*/

	/*intermediary: */
	int         code,vector_layout;
	IssmDouble *doublearray = NULL;
	int         M,N;

	/*First of, find the record for the enum, and get code  of data type: */
	this->SetFilePointerToData(&code, &vector_layout,vector_enum);

	/*Defaulting only supported for double arrays*/
	if(code!=7) _error_(EnumToStringx(vector_enum)<<" is not a double array");

	this->FetchData(&doublearray,&M,&N,vector_enum);

	for(int i=0;i<elements->Size();i++){
		Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
		if(!doublearray) element->AddInput(vector_enum,&default_value,P0Enum); 
		else             element->InputCreate(doublearray,this,M,N,vector_layout,vector_enum,code);//we need i to index into elements.
	}

	/*Free ressources. Pay attention to not freeing an AD mode independent variable though!:*/
	if(!this->independents[vector_enum]) xDelete<IssmDouble>(doublearray);
}
/*}}}*/
void IoModel::FetchDataToInput(Elements* elements,int vector_enum){/*{{{*/

	/*intermediary: */
	int     i;
	int     code,vector_layout;

	/*variables being fetched: */
	bool        boolean;
	int         integer;
	IssmDouble  scalar;
	char       *string           = NULL;
	IssmDouble *doublearray = NULL;
	int         M,N;

	/*First of, find the record for the enum, and get code  of data type: */
	this->SetFilePointerToData(&code, &vector_layout,vector_enum);

	switch(code){
		case 1: //boolean constant
			this->FetchData(&boolean,vector_enum);
			for(i=0;i<elements->Size();i++){
				Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
				element->InputUpdateFromConstant(boolean,vector_enum);
			}
			break;
		case 2: //integer constant
			this->FetchData(&integer,vector_enum);
			for(i=0;i<elements->Size();i++){
				Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
				element->InputUpdateFromConstant(integer,vector_enum);
			}
			break;
		case 3: //IssmDouble constant
			this->FetchData(&scalar,vector_enum);
			for(i=0;i<elements->Size();i++){
				Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
				element->InputUpdateFromConstant(scalar,vector_enum);
			}
			break; 
		case 5: //boolean vector
			this->FetchData(&doublearray,&M,&N,vector_enum); //we still have a doublearray, because it might include times in transient mode
			if(!doublearray) _error_(EnumToStringx(vector_enum)<<" not found in binary file");
			for(i=0;i<elements->Size();i++){
				Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
				element->InputCreate(doublearray,this,M,N,vector_layout,vector_enum,code);//we need i to index into elements.
			}
			break;
		case 6: //int vector
			this->FetchData(&doublearray,&M,&N,vector_enum); //we still have a doublearray, because it might include times in transient mode
			if(!doublearray) _error_(EnumToStringx(vector_enum)<<" not found in binary file");
			for(i=0;i<elements->Size();i++){
				Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
				element->InputCreate(doublearray,this,M,N,vector_layout,vector_enum,code);//we need i to index into elements.
			}
			break;
		case 7: //IssmDouble vector
			this->FetchData(&doublearray,&M,&N,vector_enum);
			if(!doublearray) _error_(EnumToStringx(vector_enum)<<" not found in binary file");
			for(i=0;i<elements->Size();i++){
				Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
				element->InputCreate(doublearray,this,M,N,vector_layout,vector_enum,code);//we need i to index into elements.
			}
			break;
		default:
			_error_("data code " << code << " not supported yet!");
			break;
	}
	/*Free ressources. Pay attention to not freeing an AD mode independent variable though!:*/
	if (!this->independents[vector_enum]) xDelete<IssmDouble>(doublearray);
	xDelete<char>(string);
}
/*}}}*/
void IoModel::LastIndex(int *pindex){/*{{{*/

	int my_rank;
	int lastindex,index;
	int record_length;

	/*recover my_rank:*/
	my_rank=IssmComm::GetRank();

	/*Go find in the binary file, the position of the data we want to fetch: */
	if(my_rank==0){

		/*First set FILE* position to the beginning of the file: */
		fseek(fid,0,SEEK_SET);

		/*Now march through file looking for the correct data identifier: */
		for(;;){
			/*Read enum for this size of first string name: */
			if(fread(&index,sizeof(int),1,fid)==0){
				/*Ok, we have reached the end of the file. break: */
				break;
			}

			/*read the record length, and use it to skip this record: */
			if(fread(&record_length,sizeof(int),1,fid)!=1) _error_("Could not read record_length");
			fseek(fid,record_length,SEEK_CUR);
			lastindex=index;
		}
	}
	/*Broadcast code and vector type: */
	ISSM_MPI_Bcast(&lastindex,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

	/*Assign output pointers:*/
	*pindex=lastindex;
}
/*}}}*/
/*FUNCTION IoModel::SetFilePointerToData{{{*/
FILE* IoModel::SetFilePointerToData(int* pcode,int* pvector_type, int data_enum){

	int my_rank;

	int found         = 0;
	int record_enum;
	int record_length;
	int record_code;       //1 to 7 number
	int vector_type   = 0; //nodal or elementary

	/*recover my_rank:*/
	my_rank=IssmComm::GetRank();

	/*Go find in the binary file, the position of the data we want to fetch: */
	if(my_rank==0){

		/*First set FILE* position to the beginning of the file: */
		fseek(fid,0,SEEK_SET);

		/*Now march through file looking for the correct data identifier: */
		for(;;){
			/*Read enum for this size of first string name: */
			if(fread(&record_enum,sizeof(int),1,fid)==0){
				/*Ok, we have reached the end of the file. break: */
				found=0;
				break;
			}

			/*Is this the record sought for? : */
			if (data_enum==record_enum){
				/*Ok, we have found the correct string. Pass the record length, and read data type code: */
				fseek(fid,sizeof(int),SEEK_CUR);
				if(fread(&record_code,sizeof(int),1,fid)!=1) _error_("Could not read record_code");

				/*if record_code points to a vector, get its type (nodal or elementary): */
				if(5<=record_code && record_code<=7){
					if(fread(&vector_type,sizeof(int),1,fid)!=1) _error_("Could not read vector_type");
				}
				found=1;
				break;
			}
			else{
				/*This is not the correct string, read the record length, and use it to skip this record: */
				if(fread(&record_length,sizeof(int),1,fid)!=1) _error_("Could not read record_length");
				/*skip: */
				fseek(fid,record_length,SEEK_CUR);
			}
		}
	}
	ISSM_MPI_Bcast(&found,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 
	if(!found)_error_("could not find data with name " << EnumToStringx(data_enum) << " in binary file");

	/*Broadcast code and vector type: */
	ISSM_MPI_Bcast(&record_code,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 
	ISSM_MPI_Bcast(&vector_type,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

	/*Assign output pointers:*/
	*pcode=record_code;
	if(pvector_type)*pvector_type=vector_type;

	return fid;
}
/*}}}*/
fpos_t* IoModel::SetFilePointersToData(int** pcodes,int** pvector_types, int* pnum_instances, int data_enum){/*{{{*/

	int     found          = 0;
	int     record_enum;
	int     record_length;
	int     record_code;           //1 to 7 number
	int     vector_type;           //1 to 7 number
	int    *vector_types   = NULL;
	int    *codes          = NULL;
	int     num_instances  = 0;
	int     counter;
	fpos_t *file_positions = NULL;

	/*recover my_rank:*/
	int my_rank=IssmComm::GetRank();

	/*Go find in the binary file, the data we want to fetch and count the number of 
	 * instances it appears: */
	if(my_rank==0){

		/*First set FILE* position to the beginning of the file: */
		fseek(fid,0,SEEK_SET);

		/*Now march through file looking for the correct data identifier: */
		for(;;){
			/*Read enum for this size of first string name: */
			if(fread(&record_enum,sizeof(int),1,fid)==0){
				/*Ok, we have reached the end of the file. break: */
				break;
			}

			/*Is this the record sought for? : */
			if (data_enum==record_enum) num_instances++;

			/*Read the record length, and use it to skip the record: */
			if(fread(&record_length,sizeof(int),1,fid)!=1) _error_("Could not read record_length");
			fseek(fid,record_length,SEEK_CUR);
		}

		/*Ok, initialize the number of file handles we are going to return: */
		if(num_instances){
			file_positions = xNew<fpos_t>(num_instances);
			codes          = xNew<int>(num_instances);
			vector_types   = xNew<int>(num_instances);
		}
	
		/*Reset FILE* position to the beginning of the file, and start again, this time saving the data information 
		 * as we find it: */
		counter=0;
		fseek(fid,0,SEEK_SET);

		for(;;){
			/*Read enum for this size of first string name: */
			if(fread(&record_enum,sizeof(int),1,fid)==0){
				/*Ok, we have reached the end of the file. break: */
				break;
			}

			/*Is this the record sought for? : */
			if (data_enum==record_enum){
				/*Ok, we have found the correct string. Pass the record length, and read data type code: */
				fseek(fid,sizeof(int),SEEK_CUR);
				if(fread(&record_code,sizeof(int),1,fid)!=1) _error_("Could not read record_code");

				/*if record_code points to a vector, get its type (nodal or elementary): */
				if(5<=record_code && record_code<=7){
					if(fread(&vector_type,sizeof(int),1,fid)!=1) _error_("Could not read vector_type");
				}
				codes[counter]        = record_code;
				vector_types[counter] = vector_type;
				fgetpos(fid,file_positions+counter);
				
				/*backup and skip over the record, as we have more work to do: */
				if(5<=record_code && record_code<=7) fseek(fid,-sizeof(int),SEEK_CUR);
				fseek(fid,-sizeof(int),SEEK_CUR);
				fseek(fid,-sizeof(int),SEEK_CUR);
				
				/*increment counter: */
				counter++;
			}

			/*Read the record length, and use it to skip this record, as it has already been processed: */
			if(fread(&record_length,sizeof(int),1,fid)!=1) _error_("Could not read record_length");
			/*skip: */
			fseek(fid,record_length,SEEK_CUR);
		}
	}

	/*Broadcast data: */
	ISSM_MPI_Bcast(&num_instances,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

	/*Assign output pointers:*/
	*pcodes         = codes;
	*pnum_instances = num_instances;
	if(pvector_types){
		*pvector_types=vector_types;
	}
	else{
		xDelete<int>(vector_types);
	}
	return file_positions;
}
/*}}}*/
