Index: /issm/branches/trunk-larour-SLPS2020/src/c/main/issm_post.cpp
===================================================================
--- /issm/branches/trunk-larour-SLPS2020/src/c/main/issm_post.cpp	(revision 25494)
+++ /issm/branches/trunk-larour-SLPS2020/src/c/main/issm_post.cpp	(revision 25495)
@@ -631,6 +631,8 @@
 					}
 					else if (xtype[counter]==3){
-						maxxs[counter]=doublemat;
-						minxs[counter]=doublemat;
+						maxxs[counter]=xNew<IssmDouble>(doublematsize); 
+						xMemCpy<IssmDouble>(maxxs[counter],doublemat,doublematsize);
+						minxs[counter]=xNew<IssmDouble>(doublematsize); 
+						xMemCpy<IssmDouble>(minxs[counter],doublemat,doublematsize);
 						xsize[counter]=doublematsize;
 					}
@@ -643,12 +645,252 @@
 					}
 					else if (xtype[counter]==3){
-						IssmDouble* newmax=minxs[counter];
-						IssmDouble* newmin=maxxs[counter];
+						IssmDouble* newmax=maxxs[counter];
+						IssmDouble* newmin=minxs[counter];
 						for(int k=0;k<doublematsize;k++){
-							newmax[k]=max(newmax[k],doublemat[k]);
-							newmin[k]=min(newmin[k],doublemat[k]);
+							if(doublemat[k]>newmax[k])newmax[k]=doublemat[k];
+							if(doublemat[k]<newmin[k])newmin[k]=doublemat[k];
 						}
-						maxxs[counter]=newmax;
-						minxs[counter]=newmin;
+					}
+					else _error_("cannot carry out statistics on type " << xtype[counter]); 
+				}
+			}
+		}
+		_printf0_("    average in time:\n");
+
+		/*Deal with average in time: */
+		fseek(fid,0,SEEK_SET);
+		for (int f=0;f<nfields;f++){
+			char* field=fields[f];
+			meanxtype[f]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[0]);
+			
+			if(meanxtype[f]==1){
+				meanxsize[f]=1;
+				IssmDouble timemean=0;
+				fseek(fid,0,SEEK_SET);
+				for (int j=0;j<nsteps;j++){
+					readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
+					timemean+=scalar/nsteps;
+				}
+					
+				/*Figure out max and min of time means: */
+				if(i==(lower_row+1)){
+					maxmeans[f]=xNewZeroInit<IssmDouble>(1); 
+					minmeans[f]=xNewZeroInit<IssmDouble>(1); 
+					*maxmeans[f]=timemean;
+					*minmeans[f]=timemean;
+				}
+				else{
+					*maxmeans[f]=max(*maxmeans[f],timemean);
+					*minmeans[f]=min(*minmeans[f],timemean);
+				}
+			}
+			else{
+				meanxsize[f]=doublematsize;
+				fseek(fid,0,SEEK_SET);
+				IssmDouble* timemean=xNewZeroInit<IssmDouble>(doublematsize);
+				for (int j=0;j<nsteps;j++){
+					readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
+					for (int k=0;k<doublematsize;k++){
+						timemean[k]+=doublemat[k]/nsteps;
+					}
+				}
+
+				if(i==(lower_row+1)){
+					maxmeans[f]=xNew<IssmDouble>(doublematsize);
+					xMemCpy<IssmDouble>(maxmeans[f],timemean,doublematsize);
+					minmeans[f]=xNew<IssmDouble>(doublematsize);
+					xMemCpy<IssmDouble>(minmeans[f],timemean,doublematsize);
+				}
+				else{
+					IssmDouble* maxx=maxmeans[f];
+					IssmDouble* minx=minmeans[f];
+
+					for(int k=0;k<doublematsize;k++){
+						maxx[k]=max(maxx[k],timemean[k]);
+						minx[k]=min(minx[k],timemean[k]);
+					}
+					maxmeans[f]=maxx;
+					minmeans[f]=minx;
+				}
+			}
+		}
+		fclose(fid);
+
+		/*delete buffer:*/
+		xDelete<char>(buffer);
+	}
+	ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD); 
+	_printf0_("Done reading files, now computing min and max.\n"); 
+
+	/*We have agregated minx and max across the cluster, now gather across the cluster onto
+	 *cpu0 and then compute statistics:*/
+	for (int f=0;f<nfields;f++){
+		int counter0=f*nsteps+0;
+		if (xtype[counter0]==1){ /*deal with scalars {{{*/
+			for (int j=0;j<nsteps;j++){
+				int counter=f*nsteps+j;
+
+				/*we are broadcasting doubles:*/
+				IssmDouble maxscalar=*maxxs[counter];
+				IssmDouble minscalar=*minxs[counter];
+				IssmDouble allmaxscalar;
+				IssmDouble allminscalar;
+				IssmDouble sumscalar_alltimes=0;
+
+				ISSM_MPI_Allreduce(&maxscalar,&allmaxscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
+				ISSM_MPI_Allreduce(&minscalar,&allminscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
+
+				/*Store broadcasted value for later computation of histograms:*/
+				*maxxs[counter]=allmaxscalar;
+				*minxs[counter]=allminscalar;
+
+			}
+		} /*}}}*/
+		else{ /*deal with arrays:{{{*/
+
+			int size=xsize[counter0];
+			for (int j=0;j<nsteps;j++){
+				int counter=f*nsteps+j;
+
+				/*we are broadcasting double arrays:*/
+				IssmDouble* maxx=maxxs[counter];
+				IssmDouble* minx=minxs[counter];
+
+				IssmDouble*  allmax=xNew<IssmDouble>(size);
+				IssmDouble*  allmin=xNew<IssmDouble>(size);
+				
+				ISSM_MPI_Allreduce(maxx,allmax,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
+				ISSM_MPI_Allreduce(minx,allmin,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
+
+				/*Store broadcasted value for later computation of histograms:*/
+				maxxs[counter]=allmax;
+				minxs[counter]=allmin;
+			}
+		} /*}}}*/
+	}
+	
+	/*Now do the same for the time mean fields:*/
+	for (int f=0;f<nfields;f++){
+		if (meanxtype[f]==1){ /*deal with scalars {{{*/
+
+			/*we are broadcasting doubles:*/
+			IssmDouble maxscalar=*maxmeans[f];
+			IssmDouble minscalar=*minmeans[f];
+			IssmDouble allmaxscalar;
+			IssmDouble allminscalar;
+
+			ISSM_MPI_Allreduce(&maxscalar,&allmaxscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
+			ISSM_MPI_Allreduce(&minscalar,&allminscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
+
+			/*Store for later use in histogram computation:*/
+			*maxmeans[f]=allmaxscalar;
+			*minmeans[f]=allminscalar;
+
+		} /*}}}*/
+		else{ /*deal with arrays:{{{*/
+
+			int size=meanxsize[f];
+
+			/*we are broadcasting double arrays:*/
+			IssmDouble* maxx=maxmeans[f];
+			IssmDouble* minx=minmeans[f];
+
+			IssmDouble*  allmax=xNew<IssmDouble>(size);
+			IssmDouble*  allmin=xNew<IssmDouble>(size);
+
+			ISSM_MPI_Allreduce(maxx,allmax,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
+			ISSM_MPI_Allreduce(minx,allmin,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
+
+			
+			/*Store for later use in histogram computation:*/
+			maxmeans[f]=allmax;
+			minmeans[f]=allmin;
+
+		} /*}}}*/
+	}
+
+	/*Now that we have the min and max, we can start binning. First allocate 
+	 * histograms, then start filling them:*/
+	IssmDouble** histogram=xNew<IssmDouble*>(nfields*nsteps);
+	IssmDouble** timehistogram=xNew<IssmDouble*>(nfields);
+
+	_printf0_("Start reading files again, this time binning values in the histogram:\n");
+	/*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");
+		if(fid==NULL)_error_("cound not open file: " << file << "\n");
+
+		/*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];
+			fseek(fid,0,SEEK_SET);
+			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){
+						IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(nbins);
+						IssmDouble ma=*maxxs[counter];
+						IssmDouble mi=*minxs[counter];
+						int index=(scalar-mi)/(ma-mi)*nbins; if (index==nbins)index--;
+						localhistogram[index]++;
+						histogram[counter]=localhistogram;
+					}
+					else if (xtype[counter]==3){
+						IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(doublematsize*nbins);
+						IssmDouble* ma=maxxs[counter];
+						IssmDouble* mi=minxs[counter];
+						for (int k=0;k<doublematsize;k++){
+							IssmDouble scalar=doublemat[k];
+							int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index--;
+							_assert_(scalar<=ma[k]); _assert_(scalar>=mi[k]); _assert_(index<nbins);
+							localhistogram[k*nbins+index]++;
+						}
+						histogram[counter]=localhistogram;
+					}
+					else _error_("cannot carry out statistics on type " << xtype[counter]); 
+				}
+				else{
+					if(xtype[counter]==1){
+						IssmDouble* localhistogram=histogram[counter];
+						IssmDouble ma=*maxxs[counter];
+						IssmDouble mi=*minxs[counter];
+						int index=(scalar-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1;
+						localhistogram[index]++;
+					}
+					else if (xtype[counter]==3){
+						IssmDouble* localhistogram=histogram[counter];
+						IssmDouble* ma=maxxs[counter];
+						IssmDouble* mi=minxs[counter];
+						for (int k=0;k<doublematsize;k++){
+							IssmDouble scalar=doublemat[k];
+							int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index=nbins-1;
+							localhistogram[k*nbins+index]++;
+						}
 					}
 					else _error_("cannot carry out statistics on type " << xtype[counter]); 
@@ -674,278 +916,4 @@
 				/*Figure out max and min of time means: */
 				if(i==(lower_row+1)){
-					maxmeans[f]=xNewZeroInit<IssmDouble>(1); 
-					minmeans[f]=xNewZeroInit<IssmDouble>(1); 
-					*maxmeans[f]=timemean;
-					*minmeans[f]=timemean;
-				}
-				else{
-					*maxmeans[f]=max(*maxmeans[f],timemean);
-					*minmeans[f]=min(*minmeans[f],timemean);
-				}
-			}
-			else{
-				fseek(fid,0,SEEK_SET);
-				IssmDouble* timemean=xNewZeroInit<IssmDouble>(doublematsize);
-				for (int j=0;j<nsteps;j++){
-					readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
-					for (int k=0;i<doublematsize;k++){
-						timemean[k]+=doublemat[k]/nsteps;
-					}
-				}
-
-				if(i==(lower_row+1)){
-					maxmeans[f]=timemean;
-					minmeans[f]=timemean;
-				}
-				else{
-					IssmDouble* maxx=maxmeans[f];
-					IssmDouble* minx=minmeans[f];
-
-					for(int k=0;k<doublematsize;k++){
-						maxx[k]=max(maxx[k],timemean[k]);
-						minx[k]=min(minx[k],timemean[k]);
-					}
-					maxmeans[f]=maxx;
-					minmeans[f]=minx;
-				}
-			}
-		}
-		fclose(fid);
-
-		/*delete buffer:*/
-		xDelete<char>(buffer);
-	}
-	ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD); 
-	_printf0_("Done reading files, now computing min and max.\n"); 
-
-	/*We have agregated minx and max across the cluster, now gather across the cluster onto
-	 *cpu0 and then compute statistics:*/
-	for (int f=0;f<nfields;f++){
-		int counter0=f*nsteps+0;
-		if (xtype[counter0]==1){ /*deal with scalars {{{*/
-			for (int j=0;j<nsteps;j++){
-				int counter=f*nsteps+j;
-
-				/*we are broadcasting doubles:*/
-				IssmDouble maxscalar=*maxxs[counter];
-				IssmDouble minscalar=*minxs[counter];
-				IssmDouble allmaxscalar;
-				IssmDouble allminscalar;
-				IssmDouble sumscalar_alltimes=0;
-
-				ISSM_MPI_Allreduce(&maxscalar,&allmaxscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
-				ISSM_MPI_Allreduce(&minscalar,&allminscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
-
-				/*add to results:*/
-				if(my_rank==0){
-					char fieldname[1000];
-
-					sprintf(fieldname,"%s%s",fields[f],"Max");
-					results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,allmaxscalar,j+1,0));
-					sprintf(fieldname,"%s%s",fields[f],"Min");
-					results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,allminscalar,j+1,0));
-				}
-
-				/*Store broadcasted value for later computation of histograms:*/
-				*maxxs[counter]=allmaxscalar;
-				*minxs[counter]=allminscalar;
-
-			}
-		} /*}}}*/
-		else{ /*deal with arrays:{{{*/
-
-			int size=xsize[counter0];
-			for (int j=0;j<nsteps;j++){
-				int counter=f*nsteps+j;
-
-				/*we are broadcasting double arrays:*/
-				IssmDouble* maxx=maxxs[counter];
-				IssmDouble* minx=minxs[counter];
-
-				IssmDouble*  allmax=xNew<IssmDouble>(size);
-				IssmDouble*  allmin=xNew<IssmDouble>(size);
-				
-				ISSM_MPI_Allreduce(maxx,allmax,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
-				ISSM_MPI_Allreduce(minx,allmin,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
-
-				/*add to results:*/
-				if(my_rank==0){
-					char fieldname[1000];
-
-					sprintf(fieldname,"%s%s",fields[f],"Max");
-					results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allmax,size,1,j+1,0));
-					sprintf(fieldname,"%s%s",fields[f],"Min");
-					results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allmin,size,1,j+1,0));
-				}
-				/*Store broadcasted value for later computation of histograms:*/
-				maxxs[counter]=allmax;
-				minxs[counter]=allmin;
-			}
-		} /*}}}*/
-	}
-	
-	/*Now do the same for the time mean fields:*/
-	for (int f=0;f<nfields;f++){
-		if (meanxtype[f]==1){ /*deal with scalars {{{*/
-
-			/*we are broadcasting doubles:*/
-			IssmDouble maxscalar=*maxmeans[f];
-			IssmDouble minscalar=*minmeans[f];
-			IssmDouble allmaxscalar;
-			IssmDouble allminscalar;
-
-			ISSM_MPI_Allreduce(&maxscalar,&allmaxscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
-			ISSM_MPI_Allreduce(&minscalar,&allminscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
-
-			/*add to results at time step 1:*/
-			if(my_rank==0){
-				char fieldname[1000];
-
-				sprintf(fieldname,"%s%s",fields[f],"TimeMeanMax");
-				results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,allmaxscalar,1,0));
-				sprintf(fieldname,"%s%s",fields[f],"TimeMeaMin");
-				results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,allminscalar,1,0));
-			}
-			/*Store for later use in histogram computation:*/
-			*maxmeans[f]=allmaxscalar;
-			*minmeans[f]=allminscalar;
-
-		} /*}}}*/
-		else{ /*deal with arrays:{{{*/
-
-			int size=meanxsize[f];
-
-			/*we are broadcasting double arrays:*/
-			IssmDouble* maxx=maxmeans[f];
-			IssmDouble* minx=minmeans[f];
-
-			IssmDouble*  allmax=xNew<IssmDouble>(size);
-			IssmDouble*  allmin=xNew<IssmDouble>(size);
-
-			ISSM_MPI_Allreduce(maxx,allmax,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
-			ISSM_MPI_Allreduce(minx,allmin,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
-
-			/*add to results at step 1:*/
-			if(my_rank==0){
-				char fieldname[1000];
-
-				sprintf(fieldname,"%s%s",fields[f],"TimeMeanMax");
-				results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allmax,size,1,1,0));
-				sprintf(fieldname,"%s%s",fields[f],"TimeMeanMin");
-				results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allmin,size,1,1,0));
-			}
-			
-			/*Store for later use in histogram computation:*/
-			maxmeans[f]=allmax;
-			minmeans[f]=allmin;
-
-		} /*}}}*/
-	}
-
-	/*Now that we have the min and max, we can start binning. First allocate 
-	 * histograms, then start filling them:*/
-	IssmDouble** histogram=xNew<IssmDouble*>(nfields*nsteps);
-	IssmDouble** timehistogram=xNew<IssmDouble*>(nfields);
-
-	_printf0_("Start reading files again, this time binning values in the histogram:\n");
-	/*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");
-		if(fid==NULL)_error_("cound not open file: " << file << "\n");
-
-		/*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];
-			fseek(fid,0,SEEK_SET);
-			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){
-						IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(nbins);
-						IssmDouble ma=*maxxs[counter];
-						IssmDouble mi=*minxs[counter];
-						int index=(scalar-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1;
-						localhistogram[index]++;
-						histogram[counter]=localhistogram;
-					}
-					else if (xtype[counter]==3){
-						IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(doublematsize*nbins);
-						IssmDouble* ma=maxxs[counter];
-						IssmDouble* mi=minxs[counter];
-						for (int k=0;k<doublematsize;k++){
-							IssmDouble scalar=doublemat[k];
-							int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index=nbins-1;
-							localhistogram[k*nbins+index]++;
-						}
-						histogram[counter]=localhistogram;
-					}
-					else _error_("cannot carry out statistics on type " << xtype[counter]); 
-				}
-				else{
-					if(xtype[counter]==1){
-						IssmDouble* localhistogram=histogram[counter];
-						IssmDouble ma=*maxxs[counter];
-						IssmDouble mi=*minxs[counter];
-						int index=(scalar-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1;
-						localhistogram[index]++;
-					}
-					else if (xtype[counter]==3){
-						IssmDouble* localhistogram=histogram[counter];
-						IssmDouble* ma=maxxs[counter];
-						IssmDouble* mi=minxs[counter];
-						for (int k=0;k<doublematsize;k++){
-							IssmDouble scalar=doublemat[k];
-							int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index=nbins-1;
-							localhistogram[k*nbins+index]++;
-						}
-					}
-					else _error_("cannot carry out statistics on type " << xtype[counter]); 
-				}
-			}
-		}
-		_printf0_("    average in time:\n");
-
-		/*Deal with average in time: */
-		fseek(fid,0,SEEK_SET);
-		for (int f=0;f<nfields;f++){
-			char* field=fields[f];
-			meanxtype[f]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[0]);
-			
-			if(meanxtype[f]==1){
-				IssmDouble timemean=0;
-				fseek(fid,0,SEEK_SET);
-				for (int j=0;j<nsteps;j++){
-					readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
-					timemean+=scalar/nsteps;
-				}
-					
-				/*Figure out max and min of time means: */
-				if(i==(lower_row+1)){
 					IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(nbins); 
 					IssmDouble ma=*maxmeans[f];
@@ -968,5 +936,5 @@
 				for (int j=0;j<nsteps;j++){
 					readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
-					for (int k=0;i<doublematsize;k++){
+					for (int k=0;k<doublematsize;k++){
 						timemean[k]+=doublemat[k]/nsteps;
 					}
@@ -1025,4 +993,9 @@
 					sprintf(fieldname,"%s%s",fields[f],"Histogram");
 					results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,1,nbins,j+1,0));
+
+					sprintf(fieldname,"%s%s",fields[f],"Max");
+					results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,*maxxs[counter],j+1,0));
+					sprintf(fieldname,"%s%s",fields[f],"Min");
+					results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,*minxs[counter],j+1,0));
 				}
 			}
@@ -1038,5 +1011,5 @@
 				IssmDouble*  allhisto=xNew<IssmDouble>(size*nbins);
 				
-				ISSM_MPI_Allreduce(histo,allhisto,size*nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
+				ISSM_MPI_Allreduce(histo,allhisto,size*nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
 
 				/*add to results:*/
@@ -1046,4 +1019,9 @@
 					sprintf(fieldname,"%s%s",fields[f],"Histogram");
 					results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,size,nbins,j+1,0));
+
+					sprintf(fieldname,"%s%s",fields[f],"Max");
+					results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,maxxs[counter],size,1,j+1,0));
+					sprintf(fieldname,"%s%s",fields[f],"Min");
+					results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,minxs[counter],size,1,j+1,0));
 				}
 			}
@@ -1067,4 +1045,9 @@
 				sprintf(fieldname,"%s%s",fields[f],"TimeMeanHistogram");
 				results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,1,nbins,1,0));
+
+				sprintf(fieldname,"%s%s",fields[f],"TimeMeanMax");
+				results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,*maxmeans[f],1,0));
+				sprintf(fieldname,"%s%s",fields[f],"TimeMeaMin");
+				results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,*minmeans[f],1,0));
 			}
 		} /*}}}*/
@@ -1077,5 +1060,5 @@
 			IssmDouble* allhisto=xNewZeroInit<IssmDouble>(size*nbins);
 
-			ISSM_MPI_Allreduce(histo,allhisto,size*nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
+			ISSM_MPI_Allreduce(histo,allhisto,size*nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
 			/*add to results at step 1:*/
 			if(my_rank==0){
@@ -1084,4 +1067,8 @@
 				sprintf(fieldname,"%s%s",fields[f],"TimeMeanHistogram");
 				results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,size,nbins,1,0));
+				sprintf(fieldname,"%s%s",fields[f],"TimeMeanMax");
+				results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,maxmeans[f],size,1,1,0));
+				sprintf(fieldname,"%s%s",fields[f],"TimeMeanMin");
+				results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,minmeans[f],size,1,1,0));
 			}
 		} /*}}}*/