Index: /issm/branches/trunk-larour-SLPS2022/src/c/modules/QmuStatisticsx/QmuStatisticsx.cpp
===================================================================
--- /issm/branches/trunk-larour-SLPS2022/src/c/modules/QmuStatisticsx/QmuStatisticsx.cpp	(revision 27272)
+++ /issm/branches/trunk-larour-SLPS2022/src/c/modules/QmuStatisticsx/QmuStatisticsx.cpp	(revision 27273)
@@ -657,6 +657,6 @@
 				xDelete<IssmDouble>(allhisto); 
 				if(my_rank==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));
+					sprintf(fieldname,"%s%s",fields[f],"Max"); results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,*maxxs[counter],steps[j],0));
+					sprintf(fieldname,"%s%s",fields[f],"Min"); results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,*minxs[counter],steps[j],0));
 				}
 				xDelete<IssmDouble>(maxxs[counter]); 
@@ -684,6 +684,6 @@
 				xDelete<IssmDouble>(allhisto);
 				if(my_rank==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));
+					sprintf(fieldname,"%s%s",fields[f],"Max"); results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,maxxs[counter],size,1,steps[j],0));
+					sprintf(fieldname,"%s%s",fields[f],"Min"); results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,minxs[counter],size,1,steps[j],0));
 				}
 				xDelete<IssmDouble>(maxxs[counter]);
@@ -763,4 +763,5 @@
 int ComputeMeanVariance(Parameters* parameters,Results* results,int color, ISSM_MPI_Comm statcomm){  /*{{{*/
 
+	/*variables: {{{*/
 	int nsamples; 
 	char* directory=NULL;
@@ -790,4 +791,5 @@
 	int*         meantype=NULL;
 	int*         meansize=NULL;
+	/*}}}*/
 
 	/*only work on the statistical communicator: */
@@ -810,5 +812,5 @@
 	GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm());
 
-	/*Initialize arrays:*/
+	/*Initialize arrays:{{{*/
 	xs=xNew<IssmDouble*>(nfields*nsteps);
 	xs2=xNew<IssmDouble*>(nfields*nsteps);
@@ -820,8 +822,10 @@
 	meanx=xNew<IssmDouble*>(nfields);
 	meanx2=xNew<IssmDouble*>(nfields);
+	/*}}}*/
 
 	/*Start opening files:*/
 	for (int i=(lower_row+1);i<=upper_row;i++){
 		_printf0_("reading file #: " << i << "\n");
+		/*open buffer linked to file: {{{*/
 		char file[1000];
 		long int  length;
@@ -850,6 +854,6 @@
 		_printf0_("    processing file:\n");
 		fid=fmemopen(buffer, length, "rb");
-
-		/*start reading data from the buffer directly:*/
+		/*}}}*/
+		/*read x and x^2 values for each time step:{{{*/
 		for (int f=0;f<nfields;f++){
 			char* field=fields[f];
@@ -893,7 +897,6 @@
 				}
 			}
-		}
-
-		/*Deal with time mean: */
+		}/*}}}*/
+		/*Same but for time mean: {{{*/
 		for (int f=0;f<nfields;f++){
 			char* field=fields[f];
@@ -910,7 +913,8 @@
 					meanx2[f]=xNewZeroInit<IssmDouble>(doublematsize);
 					meansize[f]=doublematsize;
-				}
-			}
-			fseek(fid,0,SEEK_SET);
+					xDelete<IssmDouble>(doublemat);
+				}
+			}
+			fseek(fid,0,SEEK_SET); 
 			if(meantype[f]==1){
 				IssmDouble sc=0;
@@ -935,4 +939,5 @@
 						timemean[k]+=doublemat[k]/nsteps;
 					}
+					xDelete<IssmDouble>(doublemat);
 				}
 				for (int k=0;k<doublematsize;k++){
@@ -946,9 +951,9 @@
 			}
 
-		}
+		} /*}}}*/
+		/*cleanup:{{{*/
 		fclose(fid);
-
-		/*delete buffer:*/
 		xDelete<char>(buffer);
+		/*}}}*/
 	}
 	ISSM_MPI_Barrier(IssmComm::GetComm());
@@ -972,4 +977,8 @@
 				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());
+				
+				xDelete<IssmDouble>(xs[counter]);
+				xDelete<IssmDouble>(xs2[counter]);
+
 				/*Build average and standard deviation. For standard deviation, use the 
 				 *following formula: sigma^2=E(x^2)-mu^2:*/
@@ -982,7 +991,7 @@
 
 					sprintf(fieldname,"%s%s",fields[f],"Mean");
-					results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,mean,j+1,0));
+					results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,mean,steps[j],0));
 					sprintf(fieldname,"%s%s",fields[f],"Stddev");
-					results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,stddev,j+1,0));
+					results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,stddev,steps[j],0));
 				}
 
@@ -1008,4 +1017,7 @@
 				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());
+				
+				xDelete<IssmDouble>(xs[counter]);
+				xDelete<IssmDouble>(xs2[counter]);
 
 				/*Build average and standard deviation. For standard deviation, use the 
@@ -1015,4 +1027,6 @@
 					stddev[k]=sqrt(sumx2[k]/nsamples-pow(mean[k],2.0));
 				}
+				xDelete<IssmDouble>(sumx);
+				xDelete<IssmDouble>(sumx2);
 
 				/*add to results:*/
@@ -1021,9 +1035,13 @@
 
 					sprintf(fieldname,"%s%s",fields[f],"Mean");
-					results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,mean,size,1,j+1,0));
+					results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,mean,size,1,steps[j],0));
 					sprintf(fieldname,"%s%s",fields[f],"Stddev");
-					results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,stddev,size,1,j+1,0));
-				}
-			}
+					results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,stddev,size,1,steps[j],0));
+				}
+			}
+			/*Deallocate:*/
+			xDelete<IssmDouble>(mean);
+			xDelete<IssmDouble>(stddev);
+
 		} /*}}}*/
 	}
@@ -1038,4 +1056,7 @@
 			IssmDouble sumscalar;
 			IssmDouble sumscalar2;
+
+			xDelete<IssmDouble>(meanx[f]);
+			xDelete<IssmDouble>(meanx2[f]);
 
 			ISSM_MPI_Reduce(&scalar,&sumscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
@@ -1072,4 +1093,7 @@
 			ISSM_MPI_Reduce(x2,sumx2,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
 
+			xDelete<IssmDouble>(meanx[f]);
+			xDelete<IssmDouble>(meanx2[f]);
+
 			/*Build average and standard deviation. For standard deviation, use the 
 			 *following formula: sigma^2=E(x^2)-mu^2:*/
@@ -1088,10 +1112,13 @@
 				results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,stddev,size,1,steps[0],0));
 			}
+			/*Deallocate:*/
+			xDelete<IssmDouble>(sumx);
+			xDelete<IssmDouble>(sumx2);
+			xDelete<IssmDouble>(mean);
+			xDelete<IssmDouble>(stddev);
 		} /*}}}*/
 	}
-
 	_printf0_("Done with MeanVariance:\n");
 	IssmComm::SetComm(ISSM_MPI_COMM_WORLD);
-
 	return 1;
 } /*}}}*/