Index: /issm/branches/trunk-larour-SLPS2020/src/c/classes/FemModel.cpp
===================================================================
--- /issm/branches/trunk-larour-SLPS2020/src/c/classes/FemModel.cpp	(revision 25360)
+++ /issm/branches/trunk-larour-SLPS2020/src/c/classes/FemModel.cpp	(revision 25361)
@@ -54,4 +54,7 @@
 
 /*Object constructors and destructor*/
+FemModel::FemModel(void){ /*{{{*/
+	/*do nothing:*/
+} /*}}}*/
 FemModel::FemModel(int argc,char** argv,ISSM_MPI_Comm incomm,bool trace){/*{{{*/
 
Index: /issm/branches/trunk-larour-SLPS2020/src/c/classes/FemModel.h
===================================================================
--- /issm/branches/trunk-larour-SLPS2020/src/c/classes/FemModel.h	(revision 25360)
+++ /issm/branches/trunk-larour-SLPS2020/src/c/classes/FemModel.h	(revision 25361)
@@ -67,4 +67,5 @@
 
 		/*constructors, destructors: */
+		FemModel(void);
 		FemModel(int argc,char** argv,ISSM_MPI_Comm comm_init,bool trace=false);
 		FemModel(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, char* restartfilename, ISSM_MPI_Comm incomm, int solution_type,IssmPDouble* X);
Index: /issm/branches/trunk-larour-SLPS2020/src/c/main/issm_post.cpp
===================================================================
--- /issm/branches/trunk-larour-SLPS2020/src/c/main/issm_post.cpp	(revision 25360)
+++ /issm/branches/trunk-larour-SLPS2020/src/c/main/issm_post.cpp	(revision 25361)
@@ -3,43 +3,579 @@
  */ 
 
+/*includes and prototypes:*/
 #include "./issm.h"
-
-int main(int argc,char **argv){
+int ComputeMeanVariance(Parameters* parameters,Results* results);
+int ComputeSampleSeries(Parameters* parameters,Results* results);
+int readdata(IssmDouble** pdoublemat, int* pdoublematsize, 
+	   	     IssmDouble* pdouble, FILE* fid,char* field,int step);
+int OutputStatistics(Parameters* parameters,Results* results);
+
+int main(int argc,char **argv){ /*{{{*/
 	
 	char* method=NULL;
 
+	int nfields;
+	char*  string=NULL;
+	char** fields=NULL;
+	char*  field=NULL;
+	int    offset;
+	char*  stepstring=NULL;
+	int    step1,step2;
+	char*  pattern=NULL;
+	int    nsteps; 
+	int*   steps=NULL;
+	int    nindices;
+	int*   indices=NULL;
+	
+	/*datasets*/
+	Parameters  *parameters   = NULL;
+	Results  *results   = NULL;
+
 	/*Initialize environment (MPI, PETSC, MUMPS, etc ...)*/
-	ISSM_MPI_Comm comm_init=EnvironmentInit(argc,argv);
+	ISSM_MPI_Comm comm=EnvironmentInit(argc,argv);
 
 	/*First things first, store the communicator, and set it as a global variable: */
-	IssmComm::SetComm(comm_init);
+	IssmComm::SetComm(comm);
+
+	/*Initialize and retrieve parameters:{{{*/
+	parameters   = new Parameters();
+	results   = new Results();
+
+	/*Method and Solution:*/
+	method=argv[1];
+	parameters->AddObject(new IntParam(SolutionTypeEnum,StatisticsSolutionEnum));
+	parameters->AddObject(new StringParam(AnalysisTypeEnum,method));
+
+	/*Directory:*/
+	parameters->AddObject(new StringParam(DirectoryNameEnum,argv[2]));
+
+	/*Model name:*/
+	parameters->AddObject(new StringParam(InputFileNameEnum,argv[3]));
+
+	/*Number of samples:*/
+	parameters->AddObject(new IntParam(QmuNsampleEnum,atoi(argv[4])));
+
+	/*Retrieve fields:*/
+	nfields=atoi(argv[5]);
+	for(int i=0;i<nfields;i++){
+		fields=xNew<char*>(nfields);
+		string= argv[6+i];
+		field=xNew<char>(strlen(string)+1);
+		xMemCpy<char>(field,string,(strlen(string)+1));
+		fields[i]=field;
+	}
+	parameters->AddObject(new StringArrayParam(FieldsEnum,fields,nfields));
+	offset=6+nfields;
+
+	/*Retrieve time steps:*/
+	stepstring=argv[offset]; 
+	pattern=strstr(stepstring,":");
+	if (pattern==NULL){
+		step1=atoi(stepstring);
+		step2=step1;
+	}
+	else{
+		step2=atoi(pattern+1);
+		stepstring[pattern-stepstring]='\0';
+		step1=atoi(stepstring);
+	}
+	nsteps=step2-step1+1; 
+	steps=xNew<int>(nsteps);
+	for (int i=step1;i<=step2;i++)steps[i-step1]=i;
+	parameters->AddObject(new IntVecParam(StepsEnum,steps,nsteps));
+
+	offset++;
+	/*}}}*/
+	
+	/*Key off method:*/
+	if(strcmp(method,"MeanVariance")==0){
+
+		ComputeMeanVariance(parameters,results);
+	}
+	else if(strcmp(method,"SampleSeries")==0){
+
+		/*Retrieve the vertex indices where we'll retrieve our sample time series:*/
+		nindices=atoi(argv[offset]); offset++;
+		indices=xNew<int>(nindices);
+		for (int i=0;i<nindices;i++){
+			indices[i]=atoi(argv[offset+i]);
+		}
+		parameters->AddObject(new IntVecParam(IndicesEnum,indices,nindices));
+		
+		ComputeSampleSeries(parameters,results);
+	}
+
+
+
+	else _error_("unknown method: " << method << "\n");
+
+	/*output results:*/
+	ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD); _printf0_("Output file.\n");
+	OutputStatistics(parameters,results);
+
+	/*Finalize ISSM:*/
+	ISSM_MPI_Finalize();
+
+	/*Return unix success: */
+	return 0; 
+
+} /*}}}*/
+int ComputeMeanVariance(Parameters* parameters,Results* results){ /*{{{*/
+
+	int nsamples; 
+	char* directory=NULL;
+	char* model=NULL;
+	char** fields=NULL;
+	int* steps=NULL;
+	int nsteps;
+	int nfields;
+	int range,lower_row,upper_row;
+	
+	/*intermediary:*/
+	IssmDouble* doublemat=NULL;
+	int         doublematsize;
+	IssmDouble scalar;
+
+	/*computation of average and variance itself:*/
+	IssmDouble*  x = NULL;
+	IssmDouble*  x2 = NULL;
+	IssmDouble** xs = NULL;
+	IssmDouble** xs2 = NULL;
+	int*         xtype=NULL;
+	int*         xsize=NULL;
+
+	/*Retrieve parameters:*/
+	parameters->FindParam(&nsamples,QmuNsampleEnum);
+	parameters->FindParam(&directory,DirectoryNameEnum);
+	parameters->FindParam(&model,InputFileNameEnum);
+	parameters->FindParam(&fields,&nfields,FieldsEnum);
+	parameters->FindParam(&steps,&nsteps,StepsEnum);
 
 	/*Get rank:*/
 	int my_rank=IssmComm::GetRank();
 
-	/*Initialize parameters:*/
-	Parameters  *parameters   = new Parameters();
-
-	/*Method:*/
-	method=argv[1];
-	parameters->AddObject(new StringParam(SolutionTypeEnum,argv[1]));
-
-	/*Directory:*/
-	parameters->AddObject(new StringParam(DirectoryEnum,argv[2]));
-
-	/*Model name:*/
-	parameters->AddObject(new StringParam(InputFileNameEnum,argv[3]));
-
-	/*Key off method:*/
-	if(strcmp(model,"MeanVariance")==0)
-	if 
-	switch(method){
-		case 
-
-	/*Finalize ISSM:*/
-	ISSM_MPI_Finalize();
-
-	/*Return unix success: */
-	return 0; 
+	/*Open files and read them complelety, in a distributed way:*/
+	range=DetermineLocalSize(nsamples,IssmComm::GetComm());
+	GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm());
+
+	/*Initialize arrays:*/
+	xs=xNew<IssmDouble*>(nfields*nsteps);
+	xs2=xNew<IssmDouble*>(nfields*nsteps);
+	xtype=xNew<int>(nfields*nsteps);
+	xsize=xNew<int>(nfields*nsteps);
+
+	/*Start opening files:*/
+	for (int i=(lower_row+1);i<=upper_row;i++){
+		_printf0_("reading file #: " << i << "\n");
+		char file[1000];
+		int  length;
+		char* buffer=NULL;
+
+		/*string:*/
+		sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i);
+
+		/*open file: */
+		_printf0_("    opening file:\n");
+		FILE* fid=fopen(file,"rb");
+
+		/*figure out size of file, and read the whole thing:*/
+		_printf0_("    reading file:\n");
+		fseek (fid, 0, SEEK_END);
+		length = ftell (fid);
+		fseek (fid, 0, SEEK_SET);
+		buffer = xNew<char>(length);
+		fread (buffer, sizeof(char), length, fid);
+
+		/*close file:*/
+		fclose (fid);
+
+		/*create a memory stream with this buffer:*/
+		_printf0_("    processing file:\n");
+		fid=fmemopen(buffer, length, "rb");
+
+		/*start reading data from the buffer directly:*/
+		for (int f=0;f<nfields;f++){
+			char* field=fields[f];
+			for (int j=0;j<nsteps;j++){
+				int counter=f*nsteps+j;
+				xtype[counter]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
+				if(i==(lower_row+1)){
+					if(xtype[counter]==1){
+						xs[counter]=xNew<IssmDouble>(1); 
+						xs2[counter]=xNew<IssmDouble>(1); 
+						*xs[counter]=scalar;
+						*xs2[counter]=pow(scalar,2.0);
+						xsize[counter]=1;
+					}
+					else if (xtype[counter]==3){
+						IssmDouble* doublemat2=xNew<IssmDouble>(doublematsize);
+						for(int k=0;k<doublematsize;k++)doublemat2[k]=pow(doublemat[k],2.0);
+						xs[counter]=doublemat;
+						xs2[counter]=doublemat2;
+						xsize[counter]=doublematsize;
+					}
+					else _error_("cannot carry out statistics on type " << xtype[counter]); 
+				}
+				else{
+					if(xtype[counter]==1){
+						*xs[counter]+=scalar;
+						*xs2[counter]+=pow(scalar,2.0);
+					}
+					else if (xtype[counter]==3){
+						IssmDouble* newdoublemat=xs[counter];
+						IssmDouble* newdoublemat2=xs2[counter];
+						for(int k=0;k<doublematsize;k++){
+							newdoublemat[k]+=doublemat[k];
+							newdoublemat2[k]+=pow(doublemat[k],2.0);
+						}
+						xs[counter]=newdoublemat;
+						xs2[counter]=newdoublemat2;
+					}
+					else _error_("cannot carry out statistics on type " << xtype[counter]); 
+				}
+			}
+		}
+		fclose(fid);
+
+		/*delete buffer:*/
+		xDelete<char>(buffer);
+	}
+	ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD); 
+	_printf0_("Done reading files, now computing mean and variance.\n"); 
+
+	/*We have agregated x and x^2 across the cluster, now consolidate: */
+	for (int f=0;f<nfields;f++){
+		for (int j=0;j<nsteps;j++){
+			int counter=f*nsteps+j;
+			if (xtype[counter]==1){
+				/*we are broadcasting doubles:*/
+				IssmDouble scalar=*xs[counter];
+				IssmDouble scalar2=*xs2[counter];
+				IssmDouble sumscalar;
+				IssmDouble sumscalar2;
+				IssmDouble mean,stddev;
+
+				ISSM_MPI_Reduce(&scalar,&sumscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
+				ISSM_MPI_Reduce(&scalar2,&sumscalar2,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
+				/*Build average and standard deviation. For standard deviation, use the 
+				 *following formula: sigma^2=E(x^2)-mu^2:*/
+				mean=sumscalar/nsamples;
+				stddev=sqrt(sumscalar2/nsamples-pow(mean,2.0));
+
+				/*add to results:*/
+				if(my_rank==0){
+					char fieldname[1000];
+					
+					sprintf(fieldname,"%s%s",fields[f],"Mean");
+					results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,mean,j+1,0));
+					sprintf(fieldname,"%s%s",fields[f],"Stddev");
+					results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,stddev,j+1,0));
+				}
+			}
+			else{
+				/*we are broadcasting double arrays:*/
+				x=xs[counter];
+				x2=xs2[counter];
+				int size=xsize[counter];
+				
+				IssmDouble*  sumx=xNew<IssmDouble>(size);
+				IssmDouble*  sumx2=xNew<IssmDouble>(size);
+				IssmDouble*  mean=xNew<IssmDouble>(size);
+				IssmDouble*  stddev=xNew<IssmDouble>(size);
+
+				ISSM_MPI_Reduce(x,sumx,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
+				ISSM_MPI_Reduce(x2,sumx2,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
+
+				/*Build average and standard deviation. For standard deviation, use the 
+				 *following formula: sigma^2=E(x^2)-mu^2:*/
+				for (int k=0;k<size;k++){
+					mean[k]=sumx[k]/nsamples;
+					stddev[k]=sqrt(sumx2[k]/nsamples-pow(mean[k],2.0));
+				}
+
+				/*add to results:*/
+				if(my_rank==0){
+					char fieldname[1000];
+
+					sprintf(fieldname,"%s%s",fields[f],"Mean");
+					results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,mean,size,1,j+1,0));
+					sprintf(fieldname,"%s%s",fields[f],"Stddev");
+					results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,stddev,size,1,j+1,0));
+				}
+			}
+		}
+	}
+
+
+} /*}}}*/
+int ComputeSampleSeries(Parameters* parameters,Results* results){ /*{{{*/
+
+	int nsamples; 
+	char* directory=NULL;
+	char* model=NULL;
+	char** fields=NULL;
+	int* steps=NULL;
+	int nsteps;
+	int nfields;
+	int range,lower_row,upper_row;
+	int* indices=NULL;
+	int  nindices;
+	
+	/*intermediary:*/
+	IssmDouble* doublemat=NULL;
+	int         doublematsize;
+	IssmDouble scalar;
+
+	/*computation of average and variance itself:*/
+	IssmDouble*  x = NULL;
+	IssmDouble*  allx=NULL;
+	IssmDouble** xs = NULL;
+	int*         xtype=NULL;
+	int*         xsize=NULL;
+
+	/*Retrieve parameters:*/
+	parameters->FindParam(&nsamples,QmuNsampleEnum);
+	parameters->FindParam(&directory,DirectoryNameEnum);
+	parameters->FindParam(&model,InputFileNameEnum);
+	parameters->FindParam(&fields,&nfields,FieldsEnum);
+	parameters->FindParam(&steps,&nsteps,StepsEnum);
+	parameters->FindParam(&indices,&nindices,IndicesEnum);
+
+	/*Get rank:*/
+	int my_rank=IssmComm::GetRank();
+
+	/*Open files and read them complelety, in a distributed way:*/
+	range=DetermineLocalSize(nsamples,IssmComm::GetComm());
+	GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm());
+
+	/*Initialize arrays:*/
+	xs=xNew<IssmDouble*>(nfields*nsteps);
+	xtype=xNew<int>(nfields*nsteps);
+	xsize=xNew<int>(nfields*nsteps);
+
+	/*Start opening files:*/
+	for (int i=(lower_row+1);i<=upper_row;i++){
+		_printf0_("reading file #: " << i << "\n");
+		char file[1000];
+		int  length;
+		char* buffer=NULL;
+
+		/*string:*/
+		sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i);
+
+		/*open file: */
+		_printf0_("    opening file:\n");
+		FILE* fid=fopen(file,"rb");
+
+		/*figure out size of file, and read the whole thing:*/
+		_printf0_("    reading file:\n");
+		fseek (fid, 0, SEEK_END);
+		length = ftell (fid);
+		fseek (fid, 0, SEEK_SET);
+		buffer = xNew<char>(length);
+		fread (buffer, sizeof(char), length, fid);
+
+		/*close file:*/
+		fclose (fid);
+
+		/*create a memory stream with this buffer:*/
+		_printf0_("    processing file:\n");
+		fid=fmemopen(buffer, length, "rb");
+
+		/*start reading data from the buffer directly:*/
+		for (int f=0;f<nfields;f++){
+			char* field=fields[f];
+			for (int j=0;j<nsteps;j++){
+				int counter=f*nsteps+j;
+				xtype[counter]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
+				if(i==(lower_row+1)){
+					if(xtype[counter]==1){
+						x=xNew<IssmDouble>(range);
+						x[0]=scalar;
+						xs[counter]=x;
+						xsize[counter]=range;
+					}
+					else if (xtype[counter]==3){
+						x=xNew<IssmDouble>(nindices*range);
+						for(int k=0;k<nindices;k++)x[(i-(lower_row+1))*nindices+k]=doublemat[indices[k]-1];
+						xs[counter]=x;
+						xsize[counter]=range*nindices;
+					}
+					else _error_("cannot carry out statistics on type " << xtype[counter]); 
+				}
+				else{
+					if(xtype[counter]==1){
+						x=xs[counter]; 
+						x[i-(lower_row+1)]=scalar;
+						xs[counter]=x;
+					}
+					else if (xtype[counter]==3){
+						x=xs[counter];
+						for(int k=0;k<nindices;k++)x[(i-(lower_row+1))*nindices+k]=doublemat[indices[k]-1];
+						xs[counter]=x;
+					}
+					else _error_("cannot carry out statistics on type " << xtype[counter]); 
+				}
+			}
+		}
+		fclose(fid);
+
+		/*delete buffer:*/
+		xDelete<char>(buffer);
+	}
+	ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD); 
+	_printf0_("Done reading files, now assembling time series.\n");
+
+	for (int f=0;f<nfields;f++){
+		for (int j=0;j<nsteps;j++){
+			int counter=f*nsteps+j;
+			if (xtype[counter]==1){
+				/*we are broadcasting range times doubles:*/
+				x=xs[counter]; 
+				allx=xNew<IssmDouble>(nsamples);
+				MPI_Gather(x, range, ISSM_MPI_PDOUBLE,allx, 1, ISSM_MPI_PDOUBLE, 0, IssmComm::GetComm());
+				/*add to results:*/
+				if(my_rank==0){
+					char fieldname[1000];
+					
+					sprintf(fieldname,"%s%s",fields[f],"Samples");
+					results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allx,nsamples,1,j+1,0));
+				}
+			}
+			else{
+				/*we are broadcasting double arrays:*/
+				x=xs[counter];
+				allx=xNew<IssmDouble>(nsamples*nindices);
+				
+				MPI_Gather(x, range*nindices, ISSM_MPI_PDOUBLE,allx, range*nindices, ISSM_MPI_PDOUBLE, 0, IssmComm::GetComm());
+				
+				/*add to results:*/
+				if(my_rank==0){
+					char fieldname[1000];
+					sprintf(fieldname,"%s%s",fields[f],"Samples");
+					results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allx,nsamples,nindices,j+1,0));
+				}
+			}
+		}
+	}
+
+
+} /*}}}*/
+int OutputStatistics(Parameters* parameters,Results* results){ /*{{{*/
+	
+	char   outputfilename[1000];
+	char* directory=NULL;
+	char* model=NULL;
+	char* method=NULL;
+	int   nsamples;
+	FemModel* femmodel=new FemModel();
+
+	
+	/*Some parameters that will allow us to use the OutputResultsx module:*/
+	parameters->AddObject(new BoolParam(QmuIsdakotaEnum,false));
+	parameters->AddObject(new BoolParam(SettingsIoGatherEnum,true));
+
+	parameters->FindParam(&directory,DirectoryNameEnum);
+	parameters->FindParam(&model,InputFileNameEnum);
+	parameters->FindParam(&method,AnalysisTypeEnum);
+	parameters->FindParam(&nsamples,QmuNsampleEnum);
+	sprintf(outputfilename,"%s/%s-%s.stats-%i",directory,model,method,nsamples);
+	parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename));
+
+	/*Call OutputResults module:*/
+	femmodel->parameters=parameters;
+	femmodel->results=results;
+
+	OutputResultsx(femmodel);
+} /*}}}*/
+int readdata(IssmDouble** pdoublemat, int* pdoublematsize, IssmDouble* pdouble, FILE* fid,char* field,int step){ /*{{{*/
+
+	int length;
+	char fieldname[1000];
+	int   fieldname_size;
+	IssmDouble   rtime;
+	int          rstep;
+	int M,N;
+
+	//fields that we retrive: 
+	IssmDouble  dfield; 
+	char*       sfield    = NULL;
+	IssmDouble* dmatfield = NULL; 
+	int*        imatfield = NULL; 
+			
+	//type of the returned field: 
+	int type;
+	int found=0;
+
+	while(1){
+
+		size_t ret_code = fread(&fieldname_size, sizeof(int), 1, fid); 
+		if(ret_code != 1) break; //we are done.
+
+		fread(fieldname, sizeof(char), fieldname_size, fid); 
+		//_printf0_("fieldname: " << fieldname << "\n");
+
+		fread(&rtime, sizeof(IssmDouble), 1, fid); 
+		fread(&rstep, sizeof(int), 1, fid); 
+
+		//check on field: 
+		if ((step==rstep) && (strcmp(field,fieldname)==0)){
+
+			//ok, go read the result really: 
+			fread(&type,sizeof(int),1,fid);
+			fread(&M,sizeof(int),1,fid);
+			if (type==1){
+				fread(&dfield,sizeof(IssmDouble),1,fid);
+			}
+			else if (type==2){
+				fread(&M,sizeof(int),1,fid);
+				sfield=xNew<char>(M);
+				fread(sfield,sizeof(char),M,fid);
+			}
+			else if (type==3){
+				fread(&N,sizeof(int),1,fid);
+				dmatfield=xNew<IssmDouble>(M*N);
+				fread(dmatfield,sizeof(IssmDouble),M*N,fid);
+			}
+			else if (type==4){
+				fread(&N,sizeof(int),1,fid);
+				imatfield=xNew<int>(M*N);
+				fread(imatfield,sizeof(int),M*N,fid);
+			}
+			else _error_("cannot read data of type " << type << "\n");
+			found=1;
+			break;
+		}
+		else{
+			//just skim to next results.
+			fread(&type,sizeof(int),1,fid);
+			fread(&M,sizeof(int),1,fid);
+			if (type==1){
+				fseek(fid,sizeof(IssmDouble),SEEK_CUR);
+			}
+			else if(type==2){
+				fseek(fid,M*sizeof(char),SEEK_CUR);
+			}
+			else if(type==3){
+				fread(&N,sizeof(int),1,fid);
+				fseek(fid,M*N*sizeof(IssmDouble),SEEK_CUR);
+			}
+			else if(type==4){
+				fread(&N,sizeof(int),1,fid);
+				fseek(fid,M*N*sizeof(int),SEEK_CUR);
+			}
+			else _error_("cannot read data of type " << type << "\n");
+		}
+	}
+	if(found==0)_error_("cound not find " << field << " at step " << step  << "\n");
+
+	/*assign output pointers:*/
+	*pdoublemat=dmatfield;
+	*pdoublematsize=M*N;
+	*pdouble=dfield;
+
+	/*return:*/
+	return type;
 
 }
+/*}}}*/