Index: /issm/branches/trunk-larour-SLPS2020/src/c/main/issm_post.cpp
===================================================================
--- /issm/branches/trunk-larour-SLPS2020/src/c/main/issm_post.cpp	(revision 25495)
+++ /issm/branches/trunk-larour-SLPS2020/src/c/main/issm_post.cpp	(revision 25496)
@@ -156,4 +156,9 @@
 	int*         xsize=NULL;
 
+	IssmDouble** meanx=NULL;
+	IssmDouble** meanx2=NULL;
+	int*         meantype=NULL;
+	int*         meansize=NULL;
+
 	/*Retrieve parameters:*/
 	parameters->FindParam(&nsamples,QmuNsampleEnum);
@@ -175,4 +180,9 @@
 	xtype=xNew<int>(nfields*nsteps);
 	xsize=xNew<int>(nfields*nsteps);
+
+	meantype=xNew<int>(nfields);
+	meansize=xNew<int>(nfields);
+	meanx=xNew<IssmDouble*>(nfields);
+	meanx2=xNew<IssmDouble*>(nfields);
 
 	/*Start opening files:*/
@@ -247,4 +257,57 @@
 			}
 		}
+		
+		/*Deal with time mean: */
+		for (int f=0;f<nfields;f++){
+			char* field=fields[f];
+			fseek(fid,0,SEEK_SET);
+			meantype[f]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[0]);
+			if(i==(lower_row+1)){
+				if(meantype[f]==1){
+					meanx[f]=xNewZeroInit<IssmDouble>(1);
+					meanx2[f]=xNewZeroInit<IssmDouble>(1);
+					meansize[f]=1;
+				}
+				else{
+					meanx[f]=xNewZeroInit<IssmDouble>(doublematsize);
+					meanx2[f]=xNewZeroInit<IssmDouble>(doublematsize);
+					meansize[f]=doublematsize;
+				}
+			}
+			fseek(fid,0,SEEK_SET);
+			if(meantype[f]==1){
+				IssmDouble sc=0;
+				IssmDouble sc2=0;
+				for(int j=0;j<nsteps;j++){
+					readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
+					sc+=scalar/nsteps;
+				}
+				sc2+=pow(sc,2.0);
+				*meanx[f]+=sc;
+				*meanx2[f]+=sc2;
+			}
+			else{
+				IssmDouble* sc=meanx[f];
+				IssmDouble* sc2=meanx2[f];
+				IssmDouble* timemean=xNewZeroInit<IssmDouble>(doublematsize);
+				IssmDouble* timemean2=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;
+					}
+				}
+				for (int k=0;k<doublematsize;k++){
+					timemean2[k]=pow(timemean[k],2.0);
+				}
+				for (int k=0;k<doublematsize;k++){
+					sc[k]+=timemean[k];
+					sc2[k]+=timemean2[k];
+				}
+
+			}
+
+		}
 		fclose(fid);
 
@@ -258,10 +321,4 @@
 	 *cpu0 and then compute statistics:*/
 	for (int f=0;f<nfields;f++){
-
-		IssmDouble sumscalar_alltimes=0;
-		IssmDouble sumscalar2_alltimes=0;
-		IssmDouble* sumarray_alltimes=NULL;
-		IssmDouble* sumarray2_alltimes=NULL;
-			
 		int counter0=f*nsteps+0;
 		if (xtype[counter0]==1){ /*deal with scalars {{{*/
@@ -293,20 +350,4 @@
 				}
 
-				/*We are also looking for the mean and stddev  over all time steps, 
-				 *so accumulate sums again:*/
-				sumscalar_alltimes+=sumscalar;
-				sumscalar2_alltimes+=sumscalar2;
-			}
-			/*Build averae and standard deviation for all time steps:*/
-			mean=sumscalar_alltimes/nsamples/nsteps;
-			stddev=sqrt(sumscalar2_alltimes/nsamples/nsteps-pow(mean,2.0));
-			/*add to results at step 1:*/
-			if(my_rank==0){
-				char fieldname[1000];
-
-				sprintf(fieldname,"%s%s",fields[f],"TimeMean");
-				results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,mean,1,0));
-				sprintf(fieldname,"%s%s",fields[f],"TimeStddev");
-				results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,stddev,1,0));
 			}
 		} /*}}}*/
@@ -314,6 +355,4 @@
 
 			int size=xsize[counter0];
-			sumarray_alltimes=xNewZeroInit<IssmDouble>(size);
-			sumarray2_alltimes=xNewZeroInit<IssmDouble>(size);
 
 			IssmDouble*  mean=xNew<IssmDouble>(size);
@@ -349,19 +388,59 @@
 					results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,stddev,size,1,j+1,0));
 				}
-
-				/*We want mean and stddev across time too, so accumulate through time: */
-				for (int k=0;k<size;k++){
-					sumarray_alltimes[k]+=sumx[k];
-					sumarray2_alltimes[k]+=sumx2[k];
-				}
-			}
-
-			/*build time mean and stddev:*/
+			}
+		} /*}}}*/
+	}
+	/*Do the same but for the time mean:*/
+	for (int f=0;f<nfields;f++){
+		if (meantype[f]==1){ /*deal with scalars {{{*/
+			IssmDouble mean,stddev;
+
+			/*we are broadcasting doubles:*/
+			IssmDouble scalar=*meanx[f];
+			IssmDouble scalar2=*meanx2[f];
+			IssmDouble sumscalar;
+			IssmDouble sumscalar2;
+
+			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],"TimeMean");
+				results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,mean,1,0));
+				sprintf(fieldname,"%s%s",fields[f],"TimeStddev");
+				results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,stddev,1,0));
+			}
+		} /*}}}*/
+		else{ /*deal with arrays:{{{*/
+
+			int size=meansize[f];
+			IssmDouble*  mean=xNew<IssmDouble>(size);
+			IssmDouble*  stddev=xNew<IssmDouble>(size);
+
+			/*we are broadcasting double arrays:*/
+			x=meanx[f];
+			x2=meanx2[f];
+
+			IssmDouble*  sumx=xNew<IssmDouble>(size);
+			IssmDouble*  sumx2=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]=sumarray_alltimes[k]/nsamples/nsteps;
-				stddev[k]=sqrt(sumarray2_alltimes[k]/nsamples/nsteps-pow(mean[k],2.0));
-			}
-
-			/*add to results at step 1:*/
+				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];
@@ -373,6 +452,6 @@
 			}
 		} /*}}}*/
-
-	}
+	}
+
 
 } /*}}}*/