/*! \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 "./objects.h"
#include "../io/io.h"
#include "./Container/Parameters.h"
#include "../shared/shared.h"
#include "../io/io.h"
#include "../include/include.h"

/*FUNCTION IoModel::IoModel(){{{*/
IoModel::IoModel(){
	this->fid=NULL;
	this->data=NULL;
	this->constants=NULL;
	
	this->my_elements=NULL;
	this->my_nodes=NULL;
	this->my_vertices=NULL;
	this->singlenodetoelementconnectivity=NULL;
	this->numbernodetoelementconnectivity=NULL;
	
	this->nodecounter=0;
	this->loadcounter=0;
	this->constraintcounter=0;
}
/*}}}*/
/*FUNCTION IoModel::IoModel(FILE*  iomodel_handle){{{*/
IoModel::IoModel(FILE* iomodel_handle){
	
	/*First, keep track of the file handle: */
	this->fid=iomodel_handle;

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

	/*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 data: */
	this->data=xNew<IssmDouble*>(MaximumNumberOfEnums);
	for(int i=0;i<MaximumNumberOfEnums;i++) this->data[i]=NULL;
	
	/*Initialize permanent data: */
	this->my_elements=NULL;
	this->my_nodes=NULL;
	this->my_vertices=NULL;
	this->singlenodetoelementconnectivity=NULL;
	this->numbernodetoelementconnectivity=NULL;
	
	this->nodecounter=0;
	this->loadcounter=0;
	this->constraintcounter=0;
}
/*}}}*/
/*FUNCTION IoModel::~IoModel(){{{*/
IoModel::~IoModel(){

	if(this->constants) delete this->constants;

	/*Some checks in debugging mode*/
	#ifdef _ISSM_DEBUG_
	if(this->data){
		for(int i=0;i<MaximumNumberOfEnums;i++){
			if(this->data[i]){
				_pprintLine_("Info: previous pointer of " << EnumToStringx(i) << " has not been freed (DeleteData has not been called)");
			}
		}
	}
	#endif

	xDelete<IssmDouble*>(this->data);
	xDelete<bool>(this->my_elements);
	xDelete<bool>(this->my_nodes);
	xDelete<int>(this->my_vertices);
	xDelete<int>(this->singlenodetoelementconnectivity);
	xDelete<int>(this->numbernodetoelementconnectivity);
}
/*}}}*/

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

	extern int my_rank;
	int record_enum = 0;

	/*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){
			_error2_("Marshalled file is empty");
		}
		else{
			if(record_enum!=MaximumNumberOfEnums){
				_printLine_("");
				_printLine_("=========================================================================");
				_printLine_(" Enums in marshalled file are not compatible with compiled code          ");
				_printLine_("                                                                         ");
				_printLine_("   * If you are running ISSM on a remote cluster:                        ");
				_printLine_("     make sure that you are using the same version of ISSM on your local ");
				_printLine_("     machine and remote cluster (you might need to run svn update)       ");
				_printLine_("   * If you are running ISSM on your local machine:                      ");
				_printLine_("     make sure that all the code is compiled (modules and executables)   ");
				_printLine_("   * If you are a developer and just added a new Enum:                   ");
				_printLine_("     you might need to run ./Synchronize.sh in src/c/EnumDefinitions     ");
				_printLine_("     and recompile                                                       ");
				_printLine_("=========================================================================");
				_printLine_("");
				_error2_("Enums not consistent (See error message above)");
			}
		}
	}
}
/*}}}*/
/*FUNCTION IoModel::Constant(bool* poutput,int constant_enum){{{*/
void IoModel::Constant(bool* poutput,int constant_enum){

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

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

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

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

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

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

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

	this->constants->FindParam(poutput,constant_enum);
}
/*}}}*/
/*FUNCTION IoModel::CopyConstantObject{{{*/
Param* IoModel::CopyConstantObject(int constant_enum){

	_assert_(this->constants);

	/*Find constant*/
	Param* param=(Param*)this->constants->FindParamObject(constant_enum);
	if(!param) _error2_("Constant " << EnumToStringx(constant_enum) << " not found in iomodel");

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

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

	return this->data[data_enum];
}
/*}}}*/
/*FUNCTION IoModel::DeleteData{{{*/
void  IoModel::DeleteData(int num,...){

	va_list ap;
	int     dataenum;
	int     i;
	DoubleMatParam* parameter=NULL;

	/*Go through the entire list of enums and delete the corresponding data from the iomodel-data dataset: */

	va_start(ap,num);
	for(i = 0; i <num; i++){
		dataenum=va_arg(ap, int);
		_assert_(dataenum<MaximumNumberOfEnums);
		xDelete<IssmDouble>(this->data[dataenum]);
	}
	va_end(ap);
} /*}}}*/
/*FUNCTION IoModel::FetchConstants{{{*/
void  IoModel::FetchConstants(void){

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

	/*records: */
	int  booleanint=0;
	int  integer=0;
	IssmPDouble scalar=0;
	char* string=NULL;
	int   string_size;

	/*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
				#ifdef _HAVE_MPI_
				MPI_Bcast(&record_code,1,MPI_INT,0,MPI_COMM_WORLD);  /*tell others cpus we are bailing: */
				#endif
				break;
			}
			else{
			
				/* Read the record length and the data type code: */
				fread(&record_length,sizeof(int),1,this->fid);
				fread(&record_code,sizeof(int),1,this->fid);
					
				#ifdef _HAVE_MPI_
				/*Tell other cpus what we are doing: */
				MPI_Bcast(&record_code,1,MPI_INT,0,MPI_COMM_WORLD);  /*tell other cpus what we are going to do: */

				/*Tell other cpus the name of the data, then branch according to the data type: */
				MPI_Bcast(&record_enum,1,MPI_INT,0,MPI_COMM_WORLD);  
				MPI_Bcast(&record_length,1,MPI_INT,0,MPI_COMM_WORLD);  
				#endif
				
				switch(record_code){
					case 1: 
						/*Read the boolean and broadcast it to other cpus:*/
						if(fread(&booleanint,sizeof(int),1,this->fid)!=1) _error2_("could not read boolean ");
						#ifdef _HAVE_MPI_
						MPI_Bcast(&booleanint,1,MPI_INT,0,MPI_COMM_WORLD); 
						#endif

						/*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) _error2_("could not read integer ");
						#ifdef _HAVE_MPI_
						MPI_Bcast(&integer,1,MPI_INT,0,MPI_COMM_WORLD); 
						#endif

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

						break;
					case 3:
						/*Read the scalar and broadcast it to other cpus:*/
						if(fread(&scalar,sizeof(IssmPDouble),1,this->fid)!=1) _error2_("could not read scalar ");
						#ifdef _HAVE_MPI_
						MPI_Bcast(&scalar,1,MPI_DOUBLE,0,MPI_COMM_WORLD); 
						#endif

						/*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) _error2_("could not read length of string ");
						#ifdef _HAVE_MPI_
						MPI_Bcast(&string_size,1,MPI_INT,0,MPI_COMM_WORLD); 
						#endif

						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)_error2_(" could not read string ");
							#ifdef _HAVE_MPI_
							MPI_Bcast(string,string_size,MPI_CHAR,0,MPI_COMM_WORLD); 
							#endif
						}
						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: 
						_error2_("unknown record type:" << record_code); 
						break;;
				}
			}
		}
	} //}}}
	#ifdef _HAVE_MPI_
	else{ //cpu ~0 {{{
		for(;;){ //wait on cpu 0
			MPI_Bcast(&record_code,1,MPI_INT,0,MPI_COMM_WORLD);  /*get from cpu 0 what we are going to do: */
			if(record_code==0){
				break; //we are done, break from the loop
			}
			else{
				MPI_Bcast(&record_enum,1,MPI_INT,0,MPI_COMM_WORLD);   //get from cpu 0 name of the data
				MPI_Bcast(&record_length,1,MPI_INT,0,MPI_COMM_WORLD);  
				switch(record_code){
				case 1: 
					/*boolean. get it from cpu 0 */
					MPI_Bcast(&booleanint,1,MPI_INT,0,MPI_COMM_WORLD); 
						
					/*create BoolParam: */
					this->constants->AddObject(new BoolParam(record_enum,(bool)booleanint)); //cast to a boolean
					break;

				case 2:
					/*integer. get it from cpu 0 */
					MPI_Bcast(&integer,1,MPI_INT,0,MPI_COMM_WORLD); 
						
					/*create IntParam: */
					this->constants->AddObject(new IntParam(record_enum,integer));

					break;
				case 3:
					/*scalar. get it from cpu 0 */
					MPI_Bcast(&scalar,1,MPI_DOUBLE,0,MPI_COMM_WORLD); 
						
					/*create DoubleParam: */
					this->constants->AddObject(new DoubleParam(record_enum,scalar));

					break;
				case 4: 
					MPI_Bcast(&string_size,1,MPI_INT,0,MPI_COMM_WORLD); 
					if(string_size){
						string=xNew<char>((string_size+1));
						string[string_size]='\0';

						/*Read string from cpu 0: */
						MPI_Bcast(string,string_size,MPI_CHAR,0,MPI_COMM_WORLD); 
					}
					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: 
					_error2_("unknown record type:" << record_code); 
					break;;
				}


			}
		}
	} //}}}
	#endif
}
/*}}}*/
/*FUNCTION IoModel::FetchData(bool*     pbool,int data_enum){{{*/
void  IoModel::FetchData(bool* pboolean,int data_enum){

	extern int my_rank;
	extern int num_procs;
	

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

	if(code!=1)_error2_("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) _error2_("could not read boolean ");
	}
	#ifdef _HAVE_MPI_
	MPI_Bcast(&booleanint,1,MPI_INT,0,MPI_COMM_WORLD); 
	#endif

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

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

	extern int my_rank;
	extern int num_procs;

	/*output: */
	int   integer;
	int   code;
	
	/*Set file pointer to beginning of the data: */
	fid=this->SetFilePointerToData(&code,NULL,data_enum);
	
	if(code!=2)_error2_("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) _error2_("could not read integer ");
	}

	#ifdef _HAVE_MPI_
	MPI_Bcast(&integer,1,MPI_INT,0,MPI_COMM_WORLD); 
	#endif

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


	extern int my_rank;
	extern int num_procs;
	

	/*output: */
	IssmPDouble   scalar;
	int      code;
	
	/*Set file pointer to beginning of the data: */
	fid=this->SetFilePointerToData(&code,NULL,data_enum);
	
	if(code!=3)_error2_("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)_error2_("could not read scalar ");
	}
	#ifdef _HAVE_MPI_
	MPI_Bcast(&scalar,1,MPI_DOUBLE,0,MPI_COMM_WORLD); 
	#endif

	/*Assign output pointers: */
	*pscalar=scalar;
		 
}
/*}}}*/
/*FUNCTION IoModel::FetchData(char**    pstring,int data_enum){{{*/
void  IoModel::FetchData(char** pstring,int data_enum){

	extern int my_rank;
	extern int num_procs;
	

	/*output: */
	char* string=NULL;
	int   string_size;
	int code=0;
	
	/*Set file pointer to beginning of the data: */
	fid=this->SetFilePointerToData(&code,NULL,data_enum);
	
	if(code!=4)_error2_("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) _error2_("could not read length of string ");
	}

	#ifdef _HAVE_MPI_
	MPI_Bcast(&string_size,1,MPI_INT,0,MPI_COMM_WORLD); 
	#endif

	/*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)_error2_(" could not read string ");
		}
		#ifdef _HAVE_MPI_
		MPI_Bcast(string,string_size,MPI_CHAR,0,MPI_COMM_WORLD); 
		#endif
	}
	else{
		string=xNew<char>(1);
		string[0]='\0';
	}


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

	extern int my_rank;
	extern int num_procs;
	int i,j;

	/*output: */
	int M,N;
	IssmPDouble* matrix=NULL;
	int*    integer_matrix=NULL;
	int code=0;
	int vector_type=0;
	
	
	/*Set file pointer to beginning of the data: */
	fid=this->SetFilePointerToData(&code,&vector_type,data_enum);

	if((code!=5) && (code!=6) && (code!=7))_error2_("expecting a IssmDouble, integer or boolean matrix for enum " << EnumToStringx(data_enum));
	
	/*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) _error2_("could not read number of rows for matrix ");
	}

	#ifdef _HAVE_MPI_
	MPI_Bcast(&M,1,MPI_INT,0,MPI_COMM_WORLD); 
	#endif

	if(my_rank==0){  
		if(fread(&N,sizeof(int),1,fid)!=1) _error2_("could not read number of columns for matrix ");
	}
	#ifdef _HAVE_MPI_
	MPI_Bcast(&N,1,MPI_INT,0,MPI_COMM_WORLD);
	#endif

	/*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) _error2_("could not read matrix ");
		}
		
		#ifdef _HAVE_MPI_
		MPI_Bcast(matrix,M*N,MPI_DOUBLE,0,MPI_COMM_WORLD); 
		#endif
	}

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

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

	extern int my_rank;
	extern int num_procs;

	/*output: */
	int M,N;
	IssmPDouble* matrix=NULL;
	int code=0;
	int vector_type=0;
	
	/*Set file pointer to beginning of the data: */
	fid=this->SetFilePointerToData(&code,&vector_type,data_enum);
	if((code!=5) && (code!=6) && (code!=7))_error2_("expecting a IssmDouble, integer or boolean matrix for enum " << EnumToStringx(data_enum));
	
	/*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) _error2_("could not read number of rows for matrix ");
	}
	#ifdef _HAVE_MPI_
	MPI_Bcast(&M,1,MPI_INT,0,MPI_COMM_WORLD); 
	#endif

	if(my_rank==0){  
		if(fread(&N,sizeof(int),1,fid)!=1) _error2_("could not read number of columns for matrix ");
	}
	#ifdef _HAVE_MPI_
	MPI_Bcast(&N,1,MPI_INT,0,MPI_COMM_WORLD); 
	#endif

	/*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) _error2_("could not read matrix ");
		}
		#ifdef _HAVE_MPI_
		MPI_Bcast(matrix,M*N,MPI_DOUBLE,0,MPI_COMM_WORLD); 
		#endif
	        *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;
}
/*}}}*/
/*FUNCTION IoModel::FetchData(char***   pstrings,int* pnumstrings,int data_enum){{{*/
void  IoModel::FetchData(char*** pstrings,int* pnumstrings,int data_enum){

	extern int my_rank;
	extern int num_procs;
	
	int i;

	/*output: */
	int   numstrings=0;
	char** strings=NULL;
	
	/*intermediary: */
	char* string=NULL;
	int   string_size;
	int   code;
	
	/*Set file pointer to beginning of the data: */
	fid=this->SetFilePointerToData(&code,NULL,data_enum);
	
	if(code!=9)_error2_("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) _error2_("could not read length of string array");
	}
	#ifdef _HAVE_MPI_
	MPI_Bcast(&numstrings,1,MPI_INT,0,MPI_COMM_WORLD); 
	#endif

	/*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) _error2_("could not read length of string ");
			}
			#ifdef _HAVE_MPI_
			MPI_Bcast(&string_size,1,MPI_INT,0,MPI_COMM_WORLD); 
			#endif
			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)_error2_(" could not read string ");
				}
				#ifdef _HAVE_MPI_
				MPI_Bcast(string,string_size,MPI_CHAR,0,MPI_COMM_WORLD); 
				#endif
			}
			else{
				string=xNew<char>(1);
				string[0]='\0';
			}

			strings[i]=string;
		}
	}

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

	int i;

	extern int my_rank;
	extern int num_procs;

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

	/*intermediary: */
	int     M, N;
	IssmPDouble *matrix = NULL;
	int     code;
	
	/*Set file pointer to beginning of the data: */
	fid=this->SetFilePointerToData(&code,NULL,data_enum);
	if(code!=8)_error2_("expecting a IssmDouble mat array for enum " << EnumToStringx(data_enum));
	
	/*Now fetch: */
	if(my_rank==0){  
		if(fread(&numrecords,sizeof(int),1,fid)!=1) _error2_("could not read number of records in matrix array ");
	}
	#ifdef _HAVE_MPI_
	MPI_Bcast(&numrecords,1,MPI_INT,0,MPI_COMM_WORLD); 
	#endif

	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) _error2_("could not read number of rows in " << i << "th matrix of matrix array");
			}
			#ifdef _HAVE_MPI_
			MPI_Bcast(&M,1,MPI_INT,0,MPI_COMM_WORLD); 
			#endif

			if(my_rank==0){  
				if(fread(&N,sizeof(int),1,fid)!=1) _error2_("could not read number of columns in " << i << "th matrix of matrix array");
			}
			#ifdef _HAVE_MPI_
			MPI_Bcast(&N,1,MPI_INT,0,MPI_COMM_WORLD); 
			#endif

			/*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) _error2_("could not read matrix ");
				}

				#ifdef _HAVE_MPI_
				MPI_Bcast(matrix,M*N,MPI_DOUBLE,0,MPI_COMM_WORLD); 
				#endif
				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;
}
/*}}}*/
/*FUNCTION IoModel::FetchData(Option**  poption,int data_enum){{{*/
void  IoModel::FetchData(Option** poption,int index){

	extern int my_rank;
	extern int num_procs;

	/*output: */
	int     code;
	Option *option      = NULL;
	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
			  IssmDouble *value = NULL;
			  value=xNew<IssmDouble>(1);
			  FetchData(value,index+1);
			  option = new OptionDouble();
			  ((OptionDouble*)option)->values = value;
			  option->name  = name;
			  option->numel = 1;
			  option->ndims = 1;
			  option->size  = NULL;
			  break;
			  }
		case 4: {//char
			  char* value = NULL;
			  FetchData(&value,index+1);
			  option = new OptionChar();
			  ((OptionChar*)option)->values = value;
			  option->name  = name;
			  option->numel = 1;
			  option->ndims = 1;
			  option->size  = NULL;
			  break;
			  }
		default:
			  _error2_("Option of format " << code << " not supported yet");
	}

	/*Assign output pointers: */
	*poption=option;
}
/*}}}*/
/*FUNCTION IoModel::FetchData(int num,...){{{*/
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);

		/*Some checks in debugging mode*/
		/*{{{*/
		#ifdef _ISSM_DEBUG_
		_assert_(dataenum<MaximumNumberOfEnums);
		if(this->data[dataenum]){
			_error2_("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);

}
/*}}}*/
/*FUNCTION IoModel::FetchDataToInput{{{*/
void IoModel::FetchDataToInput(Elements* elements,int vector_enum,int default_vector_enum,IssmDouble default_value){

	/*intermediary: */
	int     i;
	bool    defaulting    = false;
	bool    transient     = false;

	FILE   *fid           = NULL;
	int     code          = 0;
	int     vector_layout = 0;
	int     counter;
	int     nods;
	int     nel;
	int     numberofelements;


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

	/*Fetch parameters: */
	this->constants->FindParam(&numberofelements,MeshNumberofelementsEnum);

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

	switch(code){
		case 1: //boolean constant.  {{{
			this->FetchData(&boolean,vector_enum);

			/*Add boolean constant input to all elements: */
			counter=0;
			for (i=0;i<numberofelements;i++){
				if(this->my_elements[i]){
					Element* element=(Element*)elements->GetObjectByOffset(counter);
					element->InputCreate((IssmDouble)boolean,vector_enum,code);
					counter++;
				}
			}
			break; /*}}}*/
		case 2: //integer constant.  {{{
			this->FetchData(&integer,vector_enum);

			/*Add integer constant input to all elements: */
			counter=0;
			for (i=0;i<numberofelements;i++){
				if(this->my_elements[i]){
					Element* element=(Element*)elements->GetObjectByOffset(counter);
					element->InputCreate((IssmDouble)integer,vector_enum,code);
					counter++;
				}
			}
			break; /*}}}*/
		case 3: //IssmDouble constant.  {{{
			this->FetchData(&scalar,vector_enum);

			/*Add IssmDouble constant input to all elements: */
			counter=0;
			for (i=0;i<numberofelements;i++){
				if(this->my_elements[i]){
					Element* element=(Element*)elements->GetObjectByOffset(counter);
					element->InputCreate(scalar,vector_enum,code); 
					counter++;
				}
			}
			break; /*}}}*/
		case 5: //boolean vector.  {{{

			/*Fetch vector:*/
			this->FetchData(&IssmDoublevector,&M,&N,vector_enum); //we still have a IssmDoublevector, because it might include times in transient mode
			/*Check we got something, otherwise fetch default: */
			if(IssmDoublevector){
				defaulting=false;  //we are not defaulting, because  we do have a vector
			}
			else{
				/*do we have a default enum to go fetch? */
				if(default_vector_enum!=NoneEnum){
					/*yes. fetch it: */
					this->FetchData(&IssmDoublevector,&M,&N,default_vector_enum);
					if(IssmDoublevector){
						defaulting=false;  //we are not defaulting, because  we do have a vector
					}
					else{
						/*even the default vector is non existent. we are defaulting to the default value: */
						defaulting=true;
					}
				}
				else{
					/*we are left with the default value: */
					defaulting=true;
				}
			}

			/*Create inputs:*/
			counter=0;
			for (i=0;i<numberofelements;i++){
				if(this->my_elements[i]){
					Element* element=(Element*)elements->GetObjectByOffset(counter);
					if(defaulting) element->InputCreate(default_value,vector_enum,code); 
					else           element->InputCreate(IssmDoublevector,i,this,M,N,vector_layout,vector_enum,code);//we need i to index into elements.
					counter++;
				}
			}
			break; /*}}}*/
		case 6: //int vector{{{

			/*Fetch vector:*/
			this->FetchData(&IssmDoublevector,&M,&N,vector_enum); //we still have a IssmDoublevector, because it might include times in transient mode
			/*Check we got something, otherwise fetch default: */
			if(IssmDoublevector){
				defaulting=false;  //we are not defaulting, because  we do have a vector
			}
			else{
				/*do we have a default enum to go fetch? */
				if(default_vector_enum!=NoneEnum){
					/*yes. fetch it: */
					this->FetchData(&IssmDoublevector,&M,&N,default_vector_enum);
					if(IssmDoublevector){
						defaulting=false;  //we are not defaulting, because  we do have a vector
					}
					else{
						/*even the default vector is non existent. we are defaulting to the default value: */
						defaulting=true;
					}
				}
				else{
					/*we are left with the default value: */
					defaulting=true;
				}
			}

			/*Create inputs:*/
			counter=0;
			for (i=0;i<numberofelements;i++){
				if(this->my_elements[i]){
					Element* element=(Element*)elements->GetObjectByOffset(counter);
					if(defaulting) element->InputCreate(default_value,vector_enum,code); 
					else           element->InputCreate(IssmDoublevector,i,this,M,N,vector_layout,vector_enum,code);//we need i to index into elements.
					counter++;
				}
			}
			break; /*}}}*/
		case 7: //IssmDouble vector{{{

			/*Fetch vector:*/
			this->FetchData(&IssmDoublevector,&M,&N,vector_enum);
			/*Check we got something, otherwise fetch default: */
			if(IssmDoublevector){
				defaulting=false;  //we are not defaulting, because  we do have a vector
			}
			else{
				/*do we have a default enum to go fetch? */
				if(default_vector_enum!=NoneEnum){
					/*yes. fetch it: */
					this->FetchData(&IssmDoublevector,&M,&N,default_vector_enum);
					if(IssmDoublevector){
						defaulting=false;  //we are not defaulting, because  we do have a vector
					}
					else{
						/*even the default vector is non existent. we are defaulting to the default value: */
						defaulting=true;
					}
				}
				else{
					/*we are left with the default value: */
					defaulting=true;
				}
			}

			/*Create inputs:*/
			counter=0;
			for (i=0;i<numberofelements;i++){
				if(this->my_elements[i]){
					Element* element=(Element*)elements->GetObjectByOffset(counter);
					if(defaulting) element->InputCreate(default_value,vector_enum,code); 
					else           element->InputCreate(IssmDoublevector,i,this,M,N,vector_layout,vector_enum,code);//we need i to index into elements.
					counter++;
				}
			}

			break; /*}}}*/
		default: /*{{{*/
			_error2_("data code " << code << " not supported yet!");
			break;
			/*}}}*/
	}
	/*Free ressources:*/
	xDelete<IssmDouble>(IssmDoublevector);
	xDelete<char>(string);
}
/*FUNCTION IoModel::LastIndex{{{*/
void IoModel::LastIndex(int *pindex){

	extern int my_rank;
	int        lastindex,index;
	int        record_length;

	/*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: */
			fread(&record_length,sizeof(int),1,fid);
			fseek(fid,record_length,SEEK_CUR);
			lastindex=index;
		}
	}
	/*Broadcast code and vector type: */
#ifdef _HAVE_MPI_
	MPI_Bcast(&lastindex,1,MPI_INT,0,MPI_COMM_WORLD); 
#endif

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

	extern int my_rank;
	extern int num_procs;

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

	/*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);
				fread(&record_code,sizeof(int),1,fid);

				/*if record_code points to a vector, get its type (nodal or elementary): */
				if(5<=record_code && record_code<=7)fread(&vector_type,sizeof(int),1,fid);
				found=1;
				break;
			}
			else{
				/*This is not the correct string, read the record length, and use it to skip this record: */
				fread(&record_length,sizeof(int),1,fid);
				/*skip: */
				fseek(fid,record_length,SEEK_CUR);
			}
		}
	}
#ifdef _HAVE_MPI_
	MPI_Bcast(&found,1,MPI_INT,0,MPI_COMM_WORLD); 
	if(!found)_error2_("could not find data with name" << " " << EnumToStringx(data_enum) << " ");
#endif

	/*Broadcast code and vector type: */
#ifdef _HAVE_MPI_
	MPI_Bcast(&record_code,1,MPI_INT,0,MPI_COMM_WORLD); 
	MPI_Bcast(&vector_type,1,MPI_INT,0,MPI_COMM_WORLD); 
	if(record_code==5) MPI_Bcast(&vector_type,1,MPI_INT,0,MPI_COMM_WORLD); 
#endif

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

	return fid;
}
/*}}}*/
