diff --git a/doc/fix_temp_rescale.html b/doc/fix_temp_rescale.html
index 665f6a7e7..fc8c0ba24 100644
--- a/doc/fix_temp_rescale.html
+++ b/doc/fix_temp_rescale.html
@@ -15,17 +15,28 @@ <H3>fix temp/rescale/kk command
 </H3>
 <P><B>Syntax:</B>
 </P>
-<PRE>fix ID temp/rescale N Tstart Tstop 
+<PRE>fix ID temp/rescale N Tstart Tstop keyword value ... 
 </PRE>
-<UL><LI>ID is documented in <A HREF = "fix.html">fix</A> command
-<LI>temp/rescale = style name of this fix command
-<LI>N = thermostat every N timesteps
+<UL><LI>ID is documented in <A HREF = "fix.html">fix</A> command 
+
+<LI>temp/rescale = style name of this fix command 
+
+<LI>N = thermostat every N timesteps 
+
 <LI>Tstart,Tstop = desired temperature at start/end of run (temperature units) 
+
+<LI>zero or more keyword/args pairs may be appended 
+
+<LI>keyword = <I>ave</I> 
+
+<PRE>  ave values = <I>yes</I> or <I>no</I> 
+</PRE>
+
 </UL>
 <P><B>Examples:</B>
 </P>
 <PRE>fix 1 temp/rescale 100 300.0 300.0
-fix 5 temp/rescale 10 300.0 10.0 
+fix 5 temp/rescale 10 300.0 10.0 ave yes 
 </PRE>
 <P><B>Description:</B>
 </P>
@@ -40,9 +51,9 @@ <H3>fix temp/rescale/kk command
 for the particles.  Their rotational or vibrational degrees of freedom
 are not altered.
 </P>
-<P>Rescaling is performed every N timesteps. The target temperature is a
-ramped value between the Tstart and Tstop temperatures at the
-beginning and end of the run.
+<P>Rescaling is performed every N timesteps. The target temperature
+(Ttarget) is a ramped value between the Tstart and Tstop temperatures
+at the beginning and end of the run.
 </P>
 <P>This fix performs thermostatting on a per grid cell basis.  For each
 grid cell, the center-of-mass velocity and thermal temperature of the
@@ -55,13 +66,36 @@ <H3>fix temp/rescale/kk command
 temp/global/rescale</A> doc page for a
 command that thermostats the temperature of the global system.
 </P>
-<P>From the current thermal temperature and the current target
-temperature, a velocity scale factor is calculated.  That factor is
-applied to each of the components of the thermal velocity for each
-particle in the grid cell.
+<P>How the rescaling of particle velocities is done depends on the value
+of the <I>ave</I> keyword.
+</P>
+<P>For <I>ave</I> with a value <I>no</I> (the default), the thermal temperature
+(Tthermal) of each cell is used to compute a velocity scale factor for
+that cell, which is
+</P>
+<PRE>vscale = sqrt(Ttarget/Tthermal) 
+</PRE>
+<P>The vscale factor is applied to each of the components of the thermal
+velocity for each particle in the grid cell.  Only cells with 2 or
+more particles have their particle velocities rescaled.
+</P>
+<P>For <I>ave</I> with a value <I>yes</I>, the thermal temperatures of all the
+cells are averaged.  The average thermal temperature is simply the sum
+of cell thermal temperatures divided by the number of cells.  Cells
+with less than 2 particles or whose thermal temperature = 0.0
+contribute a thermal temperaure = Ttarget to the average.  The average
+thermal temperature (Tthermal_ave) for all cells is used to compute a
+velocity scale factor for all cells, which is
+</P>
+<PRE>vscale = sqrt(Ttarget/Tthermal_ave) 
+</PRE>
+<P>This single vscale factor is applied to each of the components of the
+thermal velocity for each particle in all the grid cells, including
+the particles in single-particle cells.
 </P>
-<P>After this rescaling, if the thermal temperature were re-computed for
-the grid cell, it would be exactly the target temperature.
+<P>After rescaling, for either <I>ave</I> = <I>no</I> or <I>yes</I>, if the thermal
+temperature were re-computed for any grid cell with more than one
+particle, it would be exactly the target temperature.
 </P>
 <HR>
 
@@ -107,6 +141,8 @@ <H3>fix temp/rescale/kk command
 </P>
 <P><A HREF = "fix_temp_global_rescale.html">fix temp/global/rescale</A>
 </P>
-<P><B>Default:</B> none
+<P><B>Default:</B>
+</P>
+<P>The default is ave = no.
 </P>
 </HTML>
diff --git a/doc/fix_temp_rescale.txt b/doc/fix_temp_rescale.txt
index 1602ff2bb..5f7a2ea99 100644
--- a/doc/fix_temp_rescale.txt
+++ b/doc/fix_temp_rescale.txt
@@ -11,17 +11,21 @@ fix temp/rescale/kk command :h3
 
 [Syntax:]
 
-fix ID temp/rescale N Tstart Tstop :pre
+fix ID temp/rescale N Tstart Tstop keyword value ... :pre
 
-ID is documented in "fix"_fix.html command
-temp/rescale = style name of this fix command
-N = thermostat every N timesteps
-Tstart,Tstop = desired temperature at start/end of run (temperature units) :ul
+ID is documented in "fix"_fix.html command :ulb,l
+temp/rescale = style name of this fix command :l
+N = thermostat every N timesteps :l
+Tstart,Tstop = desired temperature at start/end of run (temperature units) :l
+zero or more keyword/args pairs may be appended :l
+keyword = {ave} :l
+  ave values = {yes} or {no} :pre
+:ule
 
 [Examples:]
 
 fix 1 temp/rescale 100 300.0 300.0
-fix 5 temp/rescale 10 300.0 10.0 :pre
+fix 5 temp/rescale 10 300.0 10.0 ave yes :pre
 
 [Description:]
 
@@ -36,9 +40,9 @@ The rescaling is applied to only the translational degrees of freedom
 for the particles.  Their rotational or vibrational degrees of freedom
 are not altered.
 
-Rescaling is performed every N timesteps. The target temperature is a
-ramped value between the Tstart and Tstop temperatures at the
-beginning and end of the run.
+Rescaling is performed every N timesteps. The target temperature
+(Ttarget) is a ramped value between the Tstart and Tstop temperatures
+at the beginning and end of the run.
 
 This fix performs thermostatting on a per grid cell basis.  For each
 grid cell, the center-of-mass velocity and thermal temperature of the
@@ -51,13 +55,36 @@ doc page for the equations.  See the "fix
 temp/global/rescale"_fix_temp_globalrescale.html doc page for a
 command that thermostats the temperature of the global system.
 
-From the current thermal temperature and the current target
-temperature, a velocity scale factor is calculated.  That factor is
-applied to each of the components of the thermal velocity for each
-particle in the grid cell.
+How the rescaling of particle velocities is done depends on the value
+of the {ave} keyword.
 
-After this rescaling, if the thermal temperature were re-computed for
-the grid cell, it would be exactly the target temperature.
+For {ave} with a value {no} (the default), the thermal temperature
+(Tthermal) of each cell is used to compute a velocity scale factor for
+that cell, which is
+
+vscale = sqrt(Ttarget/Tthermal) :pre
+
+The vscale factor is applied to each of the components of the thermal
+velocity for each particle in the grid cell.  Only cells with 2 or
+more particles have their particle velocities rescaled.
+
+For {ave} with a value {yes}, the thermal temperatures of all the
+cells are averaged.  The average thermal temperature is simply the sum
+of cell thermal temperatures divided by the number of cells.  Cells
+with less than 2 particles or whose thermal temperature = 0.0
+contribute a thermal temperaure = Ttarget to the average.  The average
+thermal temperature (Tthermal_ave) for all cells is used to compute a
+velocity scale factor for all cells, which is
+
+vscale = sqrt(Ttarget/Tthermal_ave) :pre
+
+This single vscale factor is applied to each of the components of the
+thermal velocity for each particle in all the grid cells, including
+the particles in single-particle cells.
+
+After rescaling, for either {ave} = {no} or {yes}, if the thermal
+temperature were re-computed for any grid cell with more than one
+particle, it would be exactly the target temperature.
 
 :line
 
@@ -103,4 +130,6 @@ effectively.
 
 "fix temp/global/rescale"_fix_temp_global_rescale.html
 
-[Default:] none
+[Default:]
+
+The default is ave = no.
diff --git a/doc/surf_collide.html b/doc/surf_collide.html
index eff6725d7..a825c58ad 100644
--- a/doc/surf_collide.html
+++ b/doc/surf_collide.html
@@ -21,7 +21,8 @@ <H3>surf_collide command
 
 <LI>args = arguments for specific style 
 
-<PRE>  <I>specular</I> or <I>specular/kk</I> args = none
+<PRE>  <I>specular</I> or <I>specular/kk</I> args = noslip (optional)
+    noslip = reflect all velocity components off surface (not just normal component)
   <I>diffuse</I> or <I>diffuse/kk</I> args = Tsurf acc
     Tsurf = temperature of surface (temperature units)
             Tsurf can be a variable (see below)
@@ -129,6 +130,13 @@ <H3>surf_collide command
 reversed with respect to the outward normal of the surface element.
 The particle's speed is unchanged.
 </P>
+<P>Specular reflection means that a particle bounce off a surface element
+reverses only the component of its velocity normal to the surface. If
+the optional <I>noslip</I> keyword is used, the particle bounce flips the sign
+of all 3 xyz components of the particle's incident velocity, so that it
+now moves in the opposite direction, creating a no slip boundary condition.
+In either case, the particle's speed is unchanged.
+</P>
 <HR>
 
 <P>The <I>diffuse</I> style computes a simple diffusive reflection model.
diff --git a/examples/thermostat/in.thermostat_ave b/examples/thermostat/in.thermostat_ave
new file mode 100644
index 000000000..b5159a432
--- /dev/null
+++ b/examples/thermostat/in.thermostat_ave
@@ -0,0 +1,46 @@
+################################################################################
+# thermal gas in a 3d box with collisions
+# the fix temp/rescale command is used to ramp and cool the temperature
+# the compute thermal/grid command is used to measure the temperature
+
+# Note:
+#  - The "comm/sort” option to the “global” command is used to match MPI runs.
+#  - The “twopass” option is used to match Kokkos runs.
+# The "comm/sort" and "twopass" options should not be used for production runs.
+################################################################################
+
+seed	    	    12345
+dimension   	    3
+global              gridcut 1.0e-5 comm/sort yes
+
+boundary	    rr rr rr
+
+create_box  	    0 0.0001 0 0.0001 0 0.0001
+create_grid 	    10 10 10
+
+balance_grid        rcb part
+
+species		    ar.species Ar
+mixture		    air Ar vstream 0.0 0.0 0.0 temp 300.0
+
+global              nrho 7.07043E22
+global              fnum 7.07043E6
+
+collide		    vss air ar.vss
+
+create_particles    air n 10000 twopass
+
+compute             1 thermal/grid all air temp
+compute             2 reduce ave c_1[1]
+
+stats		    100
+compute             temp temp
+stats_style	    step cpu np nattempt ncoll c_temp c_2
+
+timestep 	    7.00E-9
+
+fix                 tstat temp/rescale 200 300 600 ave yes
+run 		    1000
+
+fix                 tstat temp/rescale 200 600 100 ave yes
+run 		    1000
diff --git a/examples/thermostat/log.05Apr22.mpi_1.thermostat_ave b/examples/thermostat/log.05Apr22.mpi_1.thermostat_ave
new file mode 100644
index 000000000..efc47d7c9
--- /dev/null
+++ b/examples/thermostat/log.05Apr22.mpi_1.thermostat_ave
@@ -0,0 +1,188 @@
+SPARTA (7 Jan 2022)
+################################################################################
+# thermal gas in a 3d box with collisions
+# the fix temp/rescale command is used to ramp and cool the temperature
+# the compute thermal/grid command is used to measure the temperature
+
+# Note:
+#  - The "comm/sort” option to the “global” command is used to match MPI runs.
+#  - The “twopass” option is used to match Kokkos runs.
+# The "comm/sort" and "twopass" options should not be used for production runs.
+################################################################################
+
+seed	    	    12345
+dimension   	    3
+global              gridcut 1.0e-5 comm/sort yes
+
+boundary	    rr rr rr
+
+create_box  	    0 0.0001 0 0.0001 0 0.0001
+Created orthogonal box = (0 0 0) to (0.0001 0.0001 0.0001)
+create_grid 	    10 10 10
+Created 1000 child grid cells
+  CPU time = 0.0025177 secs
+  create/ghost percent = 37.2633 62.7367
+
+balance_grid        rcb part
+Balance grid migrated 0 cells
+  CPU time = 0.00149417 secs
+  reassign/sort/migrate/ghost percent = 17.0895 0.622307 9.73352 72.5547
+
+species		    ar.species Ar
+mixture		    air Ar vstream 0.0 0.0 0.0 temp 300.0
+
+global              nrho 7.07043E22
+global              fnum 7.07043E6
+
+collide		    vss air ar.vss
+
+create_particles    air n 10000 twopass
+Created 10000 particles
+  CPU time = 0.00697708 secs
+
+compute             1 thermal/grid all air temp
+compute             2 reduce ave c_1[1]
+
+stats		    100
+compute             temp temp
+stats_style	    step cpu np nattempt ncoll c_temp c_2
+
+timestep 	    7.00E-9
+
+fix                 tstat temp/rescale 200 300 600 ave yes
+run 		    1000
+Memory usage per proc in Mbytes:
+  particles (ave,min,max) = 1.6875 1.6875 1.6875
+  grid      (ave,min,max) = 1.51379 1.51379 1.51379
+  surf      (ave,min,max) = 0 0 0
+  total     (ave,min,max) = 3.24707 3.24707 3.24707
+Step CPU Np Natt Ncoll c_temp c_2 
+       0            0    10000        0        0    299.96218    269.72175 
+     100   0.17400479    10000      987      706    299.96218    266.41225 
+     200   0.35235786    10000      989      713    395.48591    359.02435 
+     300   0.53674841    10000     1056      756    395.48591    351.65308 
+     400   0.72269797    10000     1050      724    463.79879    419.49979 
+     500   0.91654754    10000     1124      807    463.79879    415.83757 
+     600    1.1086497    10000     1133      798    531.18925          480 
+     700    1.3034542    10000     1139      806    531.18925    471.33265 
+     800    1.4995565    10000     1151      822    598.61029          540 
+     900    1.6989169    10000     1136      806    598.61029    535.97582 
+    1000    1.9032514    10000     1183      850     668.9891    599.32148 
+Loop time of 1.90329 on 1 procs for 1000 steps with 10000 particles
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Move    | 0.87913    | 0.87913    | 0.87913    |   0.0 | 46.19
+Coll    | 0.88581    | 0.88581    | 0.88581    |   0.0 | 46.54
+Sort    | 0.12321    | 0.12321    | 0.12321    |   0.0 |  6.47
+Comm    | 0.0012875  | 0.0012875  | 0.0012875  |   0.0 |  0.07
+Modify  | 0.004735   | 0.004735   | 0.004735   |   0.0 |  0.25
+Output  | 0.0069201  | 0.0069201  | 0.0069201  |   0.0 |  0.36
+Other   |            | 0.002201   |            |       |  0.12
+
+Particle moves    = 10000000 (10M)
+Cells touched     = 14624452 (14.6M)
+Particle comms    = 0 (0K)
+Boundary collides = 513895 (0.514M)
+Boundary exits    = 0 (0K)
+SurfColl checks   = 0 (0K)
+SurfColl occurs   = 0 (0K)
+Surf reactions    = 0 (0K)
+Collide attempts  = 1086520 (1.09M)
+Collide occurs    = 775143 (0.775M)
+Reactions         = 0 (0K)
+Particles stuck   = 0
+Axisymm bad moves = 0
+
+Particle-moves/CPUsec/proc: 5.25406e+06
+Particle-moves/step: 10000
+Cell-touches/particle/step: 1.46245
+Particle comm iterations/step: 1
+Particle fraction communicated: 0
+Particle fraction colliding with boundary: 0.0513895
+Particle fraction exiting boundary: 0
+Surface-checks/particle/step: 0
+Surface-collisions/particle/step: 0
+Surf-reactions/particle/step: 0
+Collision-attempts/particle/step: 0.108652
+Collisions/particle/step: 0.0775143
+Reactions/particle/step: 0
+
+Particles: 10000 ave 10000 max 10000 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Cells:      1000 ave 1000 max 1000 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+GhostCell: 0 ave 0 max 0 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+EmptyCell: 0 ave 0 max 0 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+
+fix                 tstat temp/rescale 200 600 100 ave yes
+run 		    1000
+Memory usage per proc in Mbytes:
+  particles (ave,min,max) = 1.6875 1.6875 1.6875
+  grid      (ave,min,max) = 1.51379 1.51379 1.51379
+  surf      (ave,min,max) = 0 0 0
+  total     (ave,min,max) = 3.24707 3.24707 3.24707
+Step CPU Np Natt Ncoll c_temp c_2 
+    1000            0    10000     1183      850     668.9891    599.32148 
+    1100   0.20216942    10000     1173      850     668.9891    597.96635 
+    1200   0.40663767    10000     1185      864    572.60781          500 
+    1300   0.60565662    10000     1170      797    572.60781    509.60331 
+    1400    0.8055687    10000     1202      841    464.77997    399.68397 
+    1500    1.0000267    10000     1222      802    464.77997    413.66814 
+    1600    1.1970019    10000     1186      783    350.92242     299.7815 
+    1700    1.3824208    10000     1204      744    350.92242    311.93353 
+    1800    1.5686469    10000     1190      721    238.99705          200 
+    1900    1.7456334    10000     1222      689    238.99705    211.86189 
+    2000    1.9235685    10000     1180      674    125.31819          100 
+Loop time of 1.92358 on 1 procs for 1000 steps with 10000 particles
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Move    | 0.87614    | 0.87614    | 0.87614    |   0.0 | 45.55
+Coll    | 0.90974    | 0.90974    | 0.90974    |   0.0 | 47.29
+Sort    | 0.123      | 0.123      | 0.123      |   0.0 |  6.39
+Comm    | 0.0012643  | 0.0012643  | 0.0012643  |   0.0 |  0.07
+Modify  | 0.0043652  | 0.0043652  | 0.0043652  |   0.0 |  0.23
+Output  | 0.0069406  | 0.0069406  | 0.0069406  |   0.0 |  0.36
+Other   |            | 0.002133   |            |       |  0.11
+
+Particle moves    = 10000000 (10M)
+Cells touched     = 14593538 (14.6M)
+Particle comms    = 0 (0K)
+Boundary collides = 510638 (0.511M)
+Boundary exits    = 0 (0K)
+SurfColl checks   = 0 (0K)
+SurfColl occurs   = 0 (0K)
+Surf reactions    = 0 (0K)
+Collide attempts  = 1183103 (1.18M)
+Collide occurs    = 771983 (0.772M)
+Reactions         = 0 (0K)
+Particles stuck   = 0
+Axisymm bad moves = 0
+
+Particle-moves/CPUsec/proc: 5.19863e+06
+Particle-moves/step: 10000
+Cell-touches/particle/step: 1.45935
+Particle comm iterations/step: 1
+Particle fraction communicated: 0
+Particle fraction colliding with boundary: 0.0510638
+Particle fraction exiting boundary: 0
+Surface-checks/particle/step: 0
+Surface-collisions/particle/step: 0
+Surf-reactions/particle/step: 0
+Collision-attempts/particle/step: 0.11831
+Collisions/particle/step: 0.0771983
+Reactions/particle/step: 0
+
+Particles: 10000 ave 10000 max 10000 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Cells:      1000 ave 1000 max 1000 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+GhostCell: 0 ave 0 max 0 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+EmptyCell: 0 ave 0 max 0 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
diff --git a/examples/thermostat/log.05Apr22.mpi_4.thermostat_ave b/examples/thermostat/log.05Apr22.mpi_4.thermostat_ave
new file mode 100644
index 000000000..762e61992
--- /dev/null
+++ b/examples/thermostat/log.05Apr22.mpi_4.thermostat_ave
@@ -0,0 +1,189 @@
+SPARTA (7 Jan 2022)
+################################################################################
+# thermal gas in a 3d box with collisions
+# the fix temp/rescale command is used to ramp and cool the temperature
+# the compute thermal/grid command is used to measure the temperature
+
+# Note:
+#  - The "comm/sort” option to the “global” command is used to match MPI runs.
+#  - The “twopass” option is used to match Kokkos runs.
+# The "comm/sort" and "twopass" options should not be used for production runs.
+################################################################################
+
+seed	    	    12345
+dimension   	    3
+global              gridcut 1.0e-5 comm/sort yes
+
+boundary	    rr rr rr
+
+create_box  	    0 0.0001 0 0.0001 0 0.0001
+Created orthogonal box = (0 0 0) to (0.0001 0.0001 0.0001)
+create_grid 	    10 10 10
+WARNING: Could not acquire nearby ghost cells b/c grid partition is not clumped (/ascldap/users/stamoor/sparta_master/src/grid.cpp:410)
+Created 1000 child grid cells
+  CPU time = 0.00113797 secs
+  create/ghost percent = 84.4123 15.5877
+
+balance_grid        rcb part
+Balance grid migrated 740 cells
+  CPU time = 0.00202799 secs
+  reassign/sort/migrate/ghost percent = 37.5147 0.329179 12.6617 49.4945
+
+species		    ar.species Ar
+mixture		    air Ar vstream 0.0 0.0 0.0 temp 300.0
+
+global              nrho 7.07043E22
+global              fnum 7.07043E6
+
+collide		    vss air ar.vss
+
+create_particles    air n 10000 twopass
+Created 10000 particles
+  CPU time = 0.00269508 secs
+
+compute             1 thermal/grid all air temp
+compute             2 reduce ave c_1[1]
+
+stats		    100
+compute             temp temp
+stats_style	    step cpu np nattempt ncoll c_temp c_2
+
+timestep 	    7.00E-9
+
+fix                 tstat temp/rescale 200 300 600 ave yes
+run 		    1000
+Memory usage per proc in Mbytes:
+  particles (ave,min,max) = 1.6875 1.6875 1.6875
+  grid      (ave,min,max) = 1.51379 1.51379 1.51379
+  surf      (ave,min,max) = 0 0 0
+  total     (ave,min,max) = 3.21274 3.21274 3.21274
+Step CPU Np Natt Ncoll c_temp c_2 
+       0            0    10000        0        0    297.86736    267.57161 
+     100  0.050527096    10000     1002      733    297.86736    263.36221 
+     200   0.10065508    10000     1044      747    395.35883          360 
+     300   0.15379119    10000     1055      747    395.35883    350.41188 
+     400   0.20689249    10000     1066      733    462.42646    419.50124 
+     500   0.26080036    10000     1115      802    462.42646    410.67507 
+     600   0.31535935    10000     1116      797    532.91753          480 
+     700   0.37099457    10000     1135      793    532.91753    470.83463 
+     800   0.42727995    10000     1132      798    600.49043          540 
+     900   0.48454762    10000     1160      824    600.49043    532.64919 
+    1000   0.54235458    10000     1197      863    671.24686          600 
+Loop time of 0.542403 on 4 procs for 1000 steps with 10000 particles
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Move    | 0.22895    | 0.22929    | 0.22966    |   0.1 | 42.27
+Coll    | 0.22996    | 0.23054    | 0.23153    |   0.1 | 42.50
+Sort    | 0.030813   | 0.030974   | 0.031146   |   0.1 |  5.71
+Comm    | 0.022657   | 0.023259   | 0.023987   |   0.3 |  4.29
+Modify  | 0.001287   | 0.0013166  | 0.0013523  |   0.1 |  0.24
+Output  | 0.0021853  | 0.0022685  | 0.0024858  |   0.3 |  0.42
+Other   |            | 0.02475    |            |       |  4.56
+
+Particle moves    = 10000000 (10M)
+Cells touched     = 14622757 (14.6M)
+Particle comms    = 339355 (0.339M)
+Boundary collides = 513928 (0.514M)
+Boundary exits    = 0 (0K)
+SurfColl checks   = 0 (0K)
+SurfColl occurs   = 0 (0K)
+Surf reactions    = 0 (0K)
+Collide attempts  = 1086754 (1.09M)
+Collide occurs    = 775604 (0.776M)
+Reactions         = 0 (0K)
+Particles stuck   = 0
+Axisymm bad moves = 0
+
+Particle-moves/CPUsec/proc: 4.60912e+06
+Particle-moves/step: 10000
+Cell-touches/particle/step: 1.46228
+Particle comm iterations/step: 1.002
+Particle fraction communicated: 0.0339355
+Particle fraction colliding with boundary: 0.0513928
+Particle fraction exiting boundary: 0
+Surface-checks/particle/step: 0
+Surface-collisions/particle/step: 0
+Surf-reactions/particle/step: 0
+Collision-attempts/particle/step: 0.108675
+Collisions/particle/step: 0.0775604
+Reactions/particle/step: 0
+
+Particles: 2500 ave 2616 max 2411 min
+Histogram: 1 0 0 2 0 0 0 0 0 1
+Cells:      250 ave 250 max 250 min
+Histogram: 4 0 0 0 0 0 0 0 0 0
+GhostCell: 172.5 ave 240 max 110 min
+Histogram: 1 0 0 0 2 0 0 0 0 1
+EmptyCell: 62.5 ave 130 max 0 min
+Histogram: 1 0 0 0 2 0 0 0 0 1
+
+fix                 tstat temp/rescale 200 600 100 ave yes
+run 		    1000
+Memory usage per proc in Mbytes:
+  particles (ave,min,max) = 1.6875 1.6875 1.6875
+  grid      (ave,min,max) = 1.51379 1.51379 1.51379
+  surf      (ave,min,max) = 0 0 0
+  total     (ave,min,max) = 3.21274 3.21274 3.21274
+Step CPU Np Natt Ncoll c_temp c_2 
+    1000            0    10000     1197      863    671.24686          600 
+    1100  0.058380365    10000     1158      813    671.24686    589.83045 
+    1200   0.11717296    10000     1199      807    579.32777          500 
+    1300   0.17433834    10000     1214      812    579.32777    510.95816 
+    1400   0.23125196    10000     1198      818    462.60677          400 
+    1500   0.28641748    10000     1202      793    462.60677    413.20264 
+    1600   0.34179258    10000     1206      800    351.85205    299.77984 
+    1700   0.39945364    10000     1189      723    351.85205    312.78497 
+    1800   0.49360871    10000     1176      722    237.83876          200 
+    1900   0.54414868    10000     1200      676    237.83876    210.98326 
+    2000   0.59484649    10000     1177      672    125.16707    99.952622 
+Loop time of 0.594903 on 4 procs for 1000 steps with 10000 particles
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Move    | 0.22819    | 0.23935    | 0.25537    |   2.2 | 40.23
+Coll    | 0.23572    | 0.2384     | 0.24061    |   0.4 | 40.07
+Sort    | 0.031091   | 0.031403   | 0.031751   |   0.1 |  5.28
+Comm    | 0.021785   | 0.022473   | 0.023264   |   0.4 |  3.78
+Modify  | 0.0012662  | 0.0012904  | 0.0013108  |   0.0 |  0.22
+Output  | 0.0021808  | 0.00231    | 0.0025542  |   0.3 |  0.39
+Other   |            | 0.05968    |            |       | 10.03
+
+Particle moves    = 10000000 (10M)
+Cells touched     = 14599168 (14.6M)
+Particle comms    = 337430 (0.337M)
+Boundary collides = 511295 (0.511M)
+Boundary exits    = 0 (0K)
+SurfColl checks   = 0 (0K)
+SurfColl occurs   = 0 (0K)
+Surf reactions    = 0 (0K)
+Collide attempts  = 1184112 (1.18M)
+Collide occurs    = 772234 (0.772M)
+Reactions         = 0 (0K)
+Particles stuck   = 0
+Axisymm bad moves = 0
+
+Particle-moves/CPUsec/proc: 4.20236e+06
+Particle-moves/step: 10000
+Cell-touches/particle/step: 1.45992
+Particle comm iterations/step: 1.002
+Particle fraction communicated: 0.033743
+Particle fraction colliding with boundary: 0.0511295
+Particle fraction exiting boundary: 0
+Surface-checks/particle/step: 0
+Surface-collisions/particle/step: 0
+Surf-reactions/particle/step: 0
+Collision-attempts/particle/step: 0.118411
+Collisions/particle/step: 0.0772234
+Reactions/particle/step: 0
+
+Particles: 2500 ave 2547 max 2432 min
+Histogram: 1 0 0 0 1 0 0 0 0 2
+Cells:      250 ave 250 max 250 min
+Histogram: 4 0 0 0 0 0 0 0 0 0
+GhostCell: 172.5 ave 240 max 110 min
+Histogram: 1 0 0 0 2 0 0 0 0 1
+EmptyCell: 62.5 ave 130 max 0 min
+Histogram: 1 0 0 0 2 0 0 0 0 1
diff --git a/src/KOKKOS/fix_temp_rescale_kokkos.cpp b/src/KOKKOS/fix_temp_rescale_kokkos.cpp
index e4c755022..2fa62213a 100644
--- a/src/KOKKOS/fix_temp_rescale_kokkos.cpp
+++ b/src/KOKKOS/fix_temp_rescale_kokkos.cpp
@@ -31,21 +31,13 @@ FixTempRescaleKokkos::FixTempRescaleKokkos(SPARTA *sparta, int narg, char **arg)
   datamask_modify = EMPTY_MASK;
 }
 
-/* ---------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+   current thermal temperature is calculated on a per-cell basis
+---------------------------------------------------------------------- */
 
-void FixTempRescaleKokkos::end_of_step()
+void FixTempRescaleKokkos::end_of_step_no_average(double t_target_in)
 {
-  if (update->ntimestep % nevery) return;
-
-  // set current t_target
-
-  double delta = update->ntimestep - update->beginstep;
-  if (delta != 0.0) delta /= update->endstep - update->beginstep;
-  t_target = tstart + delta * (tstop-tstart);
-
-  // sort particles by grid cell if needed
-
-  if (!particle->sorted) particle->sort();
+  t_target = t_target_in;
 
   // loop over grid cells and twice over particles in each cell
   // 1st pass: calc thermal temp via same logic as in ComputeThermalGrid
@@ -59,13 +51,12 @@ void FixTempRescaleKokkos::end_of_step()
   GridKokkos* grid_kk = (GridKokkos*) grid;
   d_cellcount = grid_kk->d_cellcount;
   d_plist = grid_kk->d_plist;
-  grid_kk->sync(Device,CINFO_MASK);
 
   int nglocal = grid->nlocal;
 
   copymode = 1;
 
-  Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagFixTempRescale_end_of_step>(0,nglocal),*this);
+  Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagFixTempRescale_end_of_step_no_average>(0,nglocal),*this);
 
   copymode = 0;
 }
@@ -73,37 +64,30 @@ void FixTempRescaleKokkos::end_of_step()
 /* ---------------------------------------------------------------------- */
 
 KOKKOS_INLINE_FUNCTION
-void FixTempRescaleKokkos::operator()(TagFixTempRescale_end_of_step, const int &icell) const {
+void FixTempRescaleKokkos::operator()(TagFixTempRescale_end_of_step_no_average, const int &icell) const {
 
   // loop over grid cells with more than 1 particle
 
-  int ispecies;
-  double mass;
-  double count,totmass,mvx,mvy,mvz,mvsq;
-  double invtotmass,vscale;
-  double vxcom,vycom,vzcom;
-  double t_current;
+  double totmass,mvx,mvy,mvz,mvsq;
   double *v;
 
-  const int np = d_cellcount[icell];
-  if (np <= 1) return;
+  const int count = d_cellcount[icell];
+  if (count <= 1) return;
 
-  count = 0.0;
   totmass = 0.0;
   mvx = mvy = mvz = 0.0;
   mvsq = 0.0;
 
   // 1st pass: loop over particles in cell
-  // 6 tallies per particle: N, Mass, mVx, mVy, mVz, mV^2
+  // 5 tallies per particle: Mass, mVx, mVy, mVz, mV^2
 
-  for (int n = 0; n < np; n++) {
+  for (int n = 0; n < count; n++) {
     const int ip = d_plist(icell,n);
 
-    ispecies = d_particles[ip].ispecies;
-    mass = d_species[ispecies].mass;
+    const int ispecies = d_particles[ip].ispecies;
+    const double mass = d_species[ispecies].mass;
     v = d_particles[ip].v;
 
-    count += 1.0;
     totmass += mass;
     mvx += mass*v[0];
     mvy += mass*v[1];
@@ -113,27 +97,27 @@ void FixTempRescaleKokkos::operator()(TagFixTempRescale_end_of_step, const int &
 
   // COM velocity of particles in grid cell
 
-  invtotmass = 1.0/totmass;
-  vxcom = mvx * invtotmass;
-  vycom = mvy * invtotmass;
-  vzcom = mvz * invtotmass;
+  const double invtotmass = 1.0/totmass;
+  const double vxcom = mvx * invtotmass;
+  const double vycom = mvy * invtotmass;
+  const double vzcom = mvz * invtotmass;
 
   // t_current = thermal T of particles in grid cell
 
-  t_current = mvsq - (mvx*mvx + mvy*mvy + mvz*mvz)*invtotmass;
+  double t_current = mvsq - (mvx*mvx + mvy*mvy + mvz*mvz)*invtotmass;
   t_current *= tprefactor/count;
 
   // vscale = scale factor for thermal velocity components
 
-  vscale = sqrt(t_target/t_current);
+  const double vscale = sqrt(t_target/t_current);
 
   // 2nd pass: loop over particles in cell
   // rescale thermal velocity components
 
-  for (int n = 0; n < np; n++) {
+  for (int n = 0; n < count; n++) {
     const int ip = d_plist(icell,n);
-    ispecies = d_particles[ip].ispecies;
-    mass = d_species[ispecies].mass;
+    const int ispecies = d_particles[ip].ispecies;
+    const double mass = d_species[ispecies].mass;
     v = d_particles[ip].v;
 
     v[0] = vscale*(v[0]-vxcom) + vxcom;
@@ -141,3 +125,156 @@ void FixTempRescaleKokkos::operator()(TagFixTempRescale_end_of_step, const int &
     v[2] = vscale*(v[2]-vzcom) + vzcom;
   }
 }
+
+/* ----------------------------------------------------------------------
+   current thermal temperature is averaged over all per-cell temperatures
+---------------------------------------------------------------------- */
+
+void FixTempRescaleKokkos::end_of_step_average(double t_target_in)
+{
+  t_target = t_target_in;
+
+  // loop over grid cells and twice over particles in each cell
+  // 1st pass: calc thermal temp via same logic as in ComputeThermalGrid
+  // 2nd pass: rescale thermal velocity components
+
+  ParticleKokkos* particle_kk = (ParticleKokkos*) particle;
+  particle_kk->sync(Device,PARTICLE_MASK|SPECIES_MASK);
+  d_particles = particle_kk->k_particles.d_view;
+  d_species = particle_kk->k_species.d_view;
+
+  GridKokkos* grid_kk = (GridKokkos*) grid;
+  d_cellcount = grid_kk->d_cellcount;
+  d_plist = grid_kk->d_plist;
+  d_cells = grid_kk->k_cells.d_view;
+  grid_kk->sync(Device,CELL_MASK);
+
+  int nglocal = grid->nlocal;
+
+  // resize d_vcom if needed
+
+  if (nglocal > maxgrid) {
+    d_vcom = DAT::t_float_1d_3();
+    maxgrid = nglocal + grid->nghost;
+    d_vcom = DAT::t_float_1d_3("temp/rescale:d_vcom",maxgrid);
+  }
+
+  // loop over grid cells to compute thermal T of each
+  // only unsplit and sub cells, skip split cells
+
+  REDUCE current_mine;
+
+  copymode = 1;
+
+  Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagFixTempRescale_end_of_step_average1>(0,nglocal),*this,current_mine);
+
+  double t_current_mine = current_mine.t;
+  bigint n_current_mine = current_mine.n;
+
+  // t_current = average of cellwise thermal T across all cells
+  // n_current = total # of cells contributing to t_current
+
+  double t_current;
+  MPI_Allreduce(&t_current_mine,&t_current,1,MPI_DOUBLE,MPI_SUM,world);
+
+  bigint n_current;
+  MPI_Allreduce(&n_current_mine,&n_current,1,MPI_SPARTA_BIGINT,MPI_SUM,world);
+
+  // t_current = cellwise averaged thermal T
+  // scale all particles in all cells by vscale
+
+  t_current /= n_current;
+  vscale = sqrt(t_target/t_current);
+
+  // loop over grid cells to rescale velocity of particles in each
+  // single-particle cells are also rescaled, their d_vcom = 0.0
+
+  Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagFixTempRescale_end_of_step_average2>(0,nglocal),*this);
+
+  copymode = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+KOKKOS_INLINE_FUNCTION
+void FixTempRescaleKokkos::operator()(TagFixTempRescale_end_of_step_average1, const int &icell, REDUCE &current_mine) const {
+
+  double totmass,mvx,mvy,mvz,mvsq;
+  double t_one;
+  double *v;
+
+  if (d_cells[icell].nsplit > 1) return;
+
+  const int count = d_cellcount[icell];
+  totmass = 0.0;
+  mvx = mvy = mvz = 0.0;
+  mvsq = 0.0;
+
+  // loop over particles in cell
+  // 5 tallies per particle: Mass, mVx, mVy, mVz, mV^2
+
+  for (int n = 0; n < count; n++) {
+    const int ip = d_plist(icell,n);
+    const int ispecies = d_particles[ip].ispecies;
+    const double mass = d_species[ispecies].mass;
+    v = d_particles[ip].v;
+
+    totmass += mass;
+    mvx += mass*v[0];
+    mvy += mass*v[1];
+    mvz += mass*v[2];
+    mvsq += mass * (v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
+  }
+
+  // t_one = thermal T of particles in the cell
+  // d_vcom = COM velocity of particle in the cell
+  // if cell has <= 1 particle: set t_one = t_target, d_vcom = zero
+  // likewise if t_one = 0.0: set t_one = t_target, d_vcom = zero
+  //   corner case when all particles have same velocity
+
+  if (count > 1.0) {
+    const double invtotmass = 1.0/totmass;
+    t_one = mvsq - (mvx*mvx + mvy*mvy + mvz*mvz)*invtotmass;
+    t_one *= tprefactor/count;
+    d_vcom(icell,0) = mvx * invtotmass;
+    d_vcom(icell,1) = mvy * invtotmass;
+    d_vcom(icell,2) = mvz * invtotmass;
+
+  } else {
+    t_one = t_target;
+    d_vcom(icell,0) = d_vcom(icell,1) = d_vcom(icell,2) = 0.0;
+  }
+
+  if (t_one == 0.0) {
+    t_one = t_target;
+    d_vcom(icell,0) = d_vcom(icell,1) = d_vcom(icell,2) = 0.0;
+  }
+
+  // accumulate thermal T over my cells
+
+  current_mine.t += t_one;
+  current_mine.n++;
+}
+
+/* ---------------------------------------------------------------------- */
+
+KOKKOS_INLINE_FUNCTION
+void FixTempRescaleKokkos::operator()(TagFixTempRescale_end_of_step_average2, const int &icell) const {
+
+  const int count = d_cellcount[icell];
+  if (count == 0) return;
+
+  // loop over particles in cell
+  // rescale thermal velocity components
+
+  for (int n = 0; n < count; n++) {
+    const int ip = d_plist(icell,n);
+    const int ispecies = d_particles[ip].ispecies;
+    const double mass = d_species[ispecies].mass;
+    double* v = d_particles[ip].v;
+
+    v[0] = vscale*(v[0]-d_vcom(icell,0)) + d_vcom(icell,0);
+    v[1] = vscale*(v[1]-d_vcom(icell,1)) + d_vcom(icell,1);
+    v[2] = vscale*(v[2]-d_vcom(icell,2)) + d_vcom(icell,2);
+  }
+}
diff --git a/src/KOKKOS/fix_temp_rescale_kokkos.h b/src/KOKKOS/fix_temp_rescale_kokkos.h
index fa13411b3..3db4c0263 100644
--- a/src/KOKKOS/fix_temp_rescale_kokkos.h
+++ b/src/KOKKOS/fix_temp_rescale_kokkos.h
@@ -26,25 +26,61 @@ FixStyle(temp/rescale/kk,FixTempRescaleKokkos)
 
 namespace SPARTA_NS {
 
-struct TagFixTempRescale_end_of_step{};
+struct TagFixTempRescale_end_of_step_no_average{};
+struct TagFixTempRescale_end_of_step_average1{};
+struct TagFixTempRescale_end_of_step_average2{};
 
 class FixTempRescaleKokkos : public FixTempRescale {
  public:
+
+  struct REDUCE {
+    bigint n;
+    double t;
+    KOKKOS_INLINE_FUNCTION
+    REDUCE() {
+      n = 0;
+      t = 0.0;
+    }
+    KOKKOS_INLINE_FUNCTION
+    REDUCE& operator+=(const REDUCE &rhs) {
+      n += rhs.n;
+      t += rhs.t;
+      return *this;
+    }
+
+    KOKKOS_INLINE_FUNCTION
+    void operator+=(const volatile REDUCE &rhs) volatile {
+      n += rhs.n;
+      t += rhs.t;
+    }
+  };
+
   FixTempRescaleKokkos(class SPARTA *, int, char **);
   virtual ~FixTempRescaleKokkos() {}
-  void end_of_step();
 
   KOKKOS_INLINE_FUNCTION
-  void operator()(TagFixTempRescale_end_of_step, const int&) const;
+  void operator()(TagFixTempRescale_end_of_step_no_average, const int&) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(TagFixTempRescale_end_of_step_average1, const int&, REDUCE&) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(TagFixTempRescale_end_of_step_average2, const int&) const;
 
  private:
-  double t_target;
+  double t_target,vscale;
+
+  DAT::t_float_1d_3 d_vcom;
 
   t_particle_1d d_particles;
   t_species_1d d_species;
 
   DAT::t_int_1d d_cellcount;
   DAT::t_int_2d d_plist;
+  t_cell_1d d_cells;
+
+  void end_of_step_no_average(double);
+  void end_of_step_average(double);
 };
 
 }
diff --git a/src/KOKKOS/kokkos_type.h b/src/KOKKOS/kokkos_type.h
index d88f68780..8cb50d887 100644
--- a/src/KOKKOS/kokkos_type.h
+++ b/src/KOKKOS/kokkos_type.h
@@ -780,7 +780,7 @@ typedef tdual_bigint_1d::t_host_um t_bigint_1d_um;
 typedef tdual_bigint_1d::t_host_const_um t_bigint_1d_const_um;
 typedef tdual_bigint_1d::t_host_const_randomread t_bigint_1d_randomread;
 
-typedef Kokkos::DualView<int*[3], Kokkos::LayoutRight, DeviceType> tdual_int_1d_3;
+typedef Kokkos::DualView<int*[3], DeviceType::array_layout, DeviceType> tdual_int_1d_3;
 typedef tdual_int_1d_3::t_host t_int_1d_3;
 typedef tdual_int_1d_3::t_host_const t_int_1d_3_const;
 typedef tdual_int_1d_3::t_host_um t_int_1d_3_um;
diff --git a/src/fix_temp_rescale.cpp b/src/fix_temp_rescale.cpp
index d85bb55b2..7c5561e57 100644
--- a/src/fix_temp_rescale.cpp
+++ b/src/fix_temp_rescale.cpp
@@ -13,10 +13,12 @@
 ------------------------------------------------------------------------- */
 
 #include "stdlib.h"
+#include "string.h"
 #include "fix_temp_rescale.h"
 #include "update.h"
 #include "grid.h"
 #include "particle.h"
+#include "memory.h"
 #include "error.h"
 
 using namespace SPARTA_NS;
@@ -26,7 +28,7 @@ using namespace SPARTA_NS;
 FixTempRescale::FixTempRescale(SPARTA *sparta, int narg, char **arg) :
   Fix(sparta, narg, arg)
 {
-  if (narg != 5) error->all(FLERR,"Illegal fix temp/rescale command");
+  if (narg < 5) error->all(FLERR,"Illegal fix temp/rescale command");
 
   nevery = atoi(arg[2]);
   tstart = atof(arg[3]);
@@ -35,6 +37,35 @@ FixTempRescale::FixTempRescale(SPARTA *sparta, int narg, char **arg) :
   if (nevery <= 0) error->all(FLERR,"Illegal fix temp/rescale command");
   if (tstart < 0.0 || tstop < 0.0)
     error->all(FLERR,"Illegal fix temp/rescale command");
+
+  // optional keyword
+
+  aveflag = 0;
+
+  int iarg = 5;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"ave") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Invalid fix temp/rescale command");
+      if (strcmp(arg[iarg+1],"yes") == 0) aveflag = 1;
+      else if (strcmp(arg[iarg+1],"no") == 0) aveflag = 0;
+      else error->all(FLERR,"Invalid fix temp/rescale command");
+      iarg += 2;
+    } else error->all(FLERR,"Invalid fix temp/rescale command");
+  }
+
+  // per-cell array for aveflag = 1 case
+
+  maxgrid = 0;
+  vcom = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixTempRescale::~FixTempRescale()
+{
+  if (copymode) return;
+
+  memory->destroy(vcom);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -69,6 +100,18 @@ void FixTempRescale::end_of_step()
 
   if (!particle->sorted) particle->sort();
 
+  // 2 variants of thermostatting
+
+  if (!aveflag) end_of_step_no_average(t_target);
+  else end_of_step_average(t_target);
+}
+
+/* ----------------------------------------------------------------------
+   current thermal temperature is calculated on a per-cell basis
+---------------------------------------------------------------------- */
+
+void FixTempRescale::end_of_step_no_average(double t_target)
+{
   // loop over grid cells and twice over particles in each cell
   // 1st pass: calc thermal temp via same logic as in ComputeThermalGrid
   // 2nd pass: rescale thermal velocity components
@@ -81,9 +124,9 @@ void FixTempRescale::end_of_step()
 
   // loop over grid cells with more than 1 particle
 
-  int ip,ispecies;
+  int ip,ispecies,count;
   double mass;
-  double count,totmass,mvx,mvy,mvz,mvsq;
+  double totmass,mvx,mvy,mvz,mvsq;
   double invtotmass,vscale;
   double vxcom,vycom,vzcom;
   double t_current;
@@ -92,13 +135,13 @@ void FixTempRescale::end_of_step()
   for (int icell = 0; icell < nglocal; icell++) {
     if (cinfo[icell].count <= 1) continue;
 
-    count = 0.0;
+    count = cinfo[icell].count;
     totmass = 0.0;
     mvx = mvy = mvz = 0.0;
     mvsq = 0.0;
 
     // 1st pass: loop over particles in cell
-    // 6 tallies per particle: N, Mass, mVx, mVy, mVz, mV^2
+    // 5 tallies per particle: Mass, mVx, mVy, mVz, mV^2
 
     ip = cinfo[icell].first;
     while (ip >= 0) {
@@ -106,7 +149,6 @@ void FixTempRescale::end_of_step()
       mass = species[ispecies].mass;
       v = particles[ip].v;
 
-      count += 1.0;
       totmass += mass;
       mvx += mass*v[0];
       mvy += mass*v[1];
@@ -149,3 +191,142 @@ void FixTempRescale::end_of_step()
     }
   }
 }
+
+/* ----------------------------------------------------------------------
+   current thermal temperature is averaged over all per-cell temperatures
+---------------------------------------------------------------------- */
+
+void FixTempRescale::end_of_step_average(double t_target)
+{
+  Particle::OnePart *particles = particle->particles;
+  Particle::Species *species = particle->species;
+  int *next = particle->next;
+  Grid::ChildCell *cells = grid->cells;
+  Grid::ChildInfo *cinfo = grid->cinfo;
+  int nglocal = grid->nlocal;
+
+  // resize vcom if needed
+
+  if (nglocal > maxgrid) {
+    memory->destroy(vcom);
+    maxgrid = nglocal + grid->nghost;
+    memory->create(vcom,maxgrid,3,"temp/rescale:vcom");
+  }
+
+  // loop over grid cells to compute thermal T of each
+  // only unsplit and sub cells, skip split cells
+
+  int ip,ispecies,count;
+  double mass;
+  double totmass,mvx,mvy,mvz,mvsq;
+  double invtotmass,t_one;
+  double *v;
+
+  bigint n_current_mine = 0.0;
+  double t_current_mine = 0.0;
+  
+  for (int icell = 0; icell < nglocal; icell++) {
+    if (cells[icell].nsplit > 1) continue;
+
+    count = cinfo[icell].count;
+    totmass = 0.0;
+    mvx = mvy = mvz = 0.0;
+    mvsq = 0.0;
+
+    // loop over particles in cell
+    // 5 tallies per particle: Mass, mVx, mVy, mVz, mV^2
+
+    ip = cinfo[icell].first;
+    while (ip >= 0) {
+      ispecies = particles[ip].ispecies;
+      mass = species[ispecies].mass;
+      v = particles[ip].v;
+
+      totmass += mass;
+      mvx += mass*v[0];
+      mvy += mass*v[1];
+      mvz += mass*v[2];
+      mvsq += mass * (v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
+
+      ip = next[ip];
+    }
+
+    // t_one = thermal T of particles in the cell
+    // vcom = COM velocity of particle in the cell
+    // if cell has <= 1 particle: set t_one = t_target, vcom = zero
+    // likewise if t_one = 0.0: set t_one = t_target, vcom = zero
+    //   corner case when all particles have same velocity
+
+    if (count > 1) {
+      invtotmass = 1.0/totmass;
+      t_one = mvsq - (mvx*mvx + mvy*mvy + mvz*mvz)*invtotmass;
+      t_one *= tprefactor/count;
+      vcom[icell][0] = mvx * invtotmass;
+      vcom[icell][1] = mvy * invtotmass;
+      vcom[icell][2] = mvz * invtotmass;
+
+    } else {
+      t_one = t_target;
+      vcom[icell][0] = vcom[icell][1] = vcom[icell][2] = 0.0;
+    }
+
+    if (t_one == 0.0) {
+      t_one = t_target;
+      vcom[icell][0] = vcom[icell][1] = vcom[icell][2] = 0.0;
+    }
+
+    // accumulate thermal T over my cells
+
+    t_current_mine += t_one;
+    n_current_mine++;
+  }
+
+  // t_current = average of cellwise thermal T across all cells
+  // n_current = total # of cells contributing to t_current
+
+  double t_current;
+  MPI_Allreduce(&t_current_mine,&t_current,1,MPI_DOUBLE,MPI_SUM,world);
+
+  bigint n_current;
+  MPI_Allreduce(&n_current_mine,&n_current,1,MPI_SPARTA_BIGINT,MPI_SUM,world);
+
+  // t_current = cellwise averaged thermal T
+  // scale all particles in all cells by vscale
+
+  t_current /= n_current;
+  double vscale = sqrt(t_target/t_current);
+
+  // loop over grid cells to rescale velocity of particles in each
+  // single-particle cells are also rescaled, their vcom = 0.0
+
+  for (int icell = 0; icell < nglocal; icell++) {
+    if (cinfo[icell].count == 0) continue;
+
+    // loop over particles in cell
+    // rescale thermal velocity components
+
+    ip = cinfo[icell].first;
+    while (ip >= 0) {
+      ispecies = particles[ip].ispecies;
+      mass = species[ispecies].mass;
+      v = particles[ip].v;
+
+      v[0] = vscale*(v[0]-vcom[icell][0]) + vcom[icell][0];
+      v[1] = vscale*(v[1]-vcom[icell][1]) + vcom[icell][1];
+      v[2] = vscale*(v[2]-vcom[icell][2]) + vcom[icell][2];
+
+      ip = next[ip];
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   memory usage
+------------------------------------------------------------------------- */
+
+double FixTempRescale::memory_usage()
+{
+  double bytes = 0.0;
+  bytes += maxgrid*3 * sizeof(double);    // vcom
+  return bytes;
+}
diff --git a/src/fix_temp_rescale.h b/src/fix_temp_rescale.h
index 47a7783fa..b8750171f 100644
--- a/src/fix_temp_rescale.h
+++ b/src/fix_temp_rescale.h
@@ -29,14 +29,22 @@ namespace SPARTA_NS {
 class FixTempRescale : public Fix {
  public:
   FixTempRescale(class SPARTA *, int, char **);
-  virtual ~FixTempRescale() {}
+  virtual ~FixTempRescale();
   int setmask();
   void init();
-  virtual void end_of_step();
+  void end_of_step();
+  double memory_usage();
 
  protected:
+  int aveflag;
   double tstart,tstop;
   double tprefactor;
+
+  int maxgrid;
+  double **vcom;
+
+  virtual void end_of_step_no_average(double);
+  virtual void end_of_step_average(double);
 };
 
 }