-    *out << O815->paraQ->getParaVals();
-    *out << "\t" << ( Sim->nequi + (nthmeas+1)*Sim->nskip );
-    *out << "\t" << kpimode;
-    *out << "\t" << 0 + int( kpimode/4. * O815->comargs.lsize[0] );
-    *out << "\t" << 0 + int( kpimode/4. * O815->comargs.lsize[1] );
-    *out << "\t" << abs( OM->phip2[ kpimode ] - 1./(ppt*ppt + ppx*ppx + Sim->m*Sim->m) ) << endl;
+    avgphip2[kpimode] += OM->phip2[ kpimode ];
+
+    if ( (nthmeas+1)%avgmany == 0 ) {
+      *out << O815->paraQ->getParaVals();
+      *out << "\t" << ( Sim->nequi + (nthmeas+1)*Sim->nskip );
+      *out << "\t" << kpimode;
+      *out << "\t" << 0 + int( kpimode/4. * O815->comargs.lsize[0] );
+      *out << "\t" << 0 + int( kpimode/4. * O815->comargs.lsize[1] );
+      *out << "\t" << abs( avgphip2[kpimode] / (double)avgmany - 1./(ppt*ppt + ppx*ppx + Sim->m*Sim->m) ) << endl;
+    }