diff --git a/inputs/z4c/awa/z4c_gauge_wave.athinput b/inputs/z4c/awa/z4c_gauge_wave.athinput
new file mode 100644
index 00000000..5e05b99d
--- /dev/null
+++ b/inputs/z4c/awa/z4c_gauge_wave.athinput
@@ -0,0 +1,76 @@
+# AthenaXXX input file for Z4c linear wave tests
+
+<comment>
+problem   = z4c linear waves
+reference = e.g. Daverio et al. arxiv:1810.12346 (2018)
+
+<job>
+basename  = gauge    # problem ID: basename of output filenames
+
+<mesh>
+nghost    = 2          # Number of ghost cells
+nx1       = 50         # Number of zones in X1-direction
+x1min     = 0.0        # minimum value of X1
+x1max     = 1.0        # maximum value of X1
+ix1_bc    = periodic   # inner-X1 boundary flag
+ox1_bc    = periodic   # outer-X1 boundary flag
+
+nx2       = 4         # Number of zones in X2-direction
+x2min     = 0.0        # minimum value of X2
+x2max     = 1.0        # maximum value of X2
+ix2_bc    = periodic   # inner-X2 boundary flag
+ox2_bc    = periodic   # outer-X2 boundary flag
+
+nx3       = 4         # Number of zones in X3-direction
+x3min     = 0.0        # minimum value of X3
+x3max     = 1.0        # maximum value of X3
+ix3_bc    = periodic   # inner-X3 boundary flag
+ox3_bc    = periodic   # outer-X3 boundary flag
+
+<meshblock>
+nx1       = 50         # Number of cells in each MeshBlock, X1-dir
+nx2       = 4         # Number of cells in each MeshBlock, X2-dir
+nx3       = 4         # Number of cells in each MeshBlock, X3-dir
+
+<time>
+evolution  = dynamic   # dynamic/kinematic/static
+integrator = rk4       # time integration algorithm
+cfl_number = 0.5       # The Courant, Friedrichs, & Lewy (CFL) Number
+nlim       = -1        # cycle limit (no limit if <0)
+tlim       = 1000.0       # time limit
+ndiag      = 1         # cycles between diagostic output
+
+<z4c>
+lapse_harmonic  = 0.0       # Harmonic lapse parameter mu_L
+lapse_oplog     = 2.0       # 1+log lapse parameter
+shift_eta       = 2.0       # Shift damping term
+diss            = 0.02       # Kreiss-Oliger dissipation parameter
+chi_div_floor   = 0.00001   # Floor on the conformal factor
+damp_kappa1     = 0.02       # Constraint damping factor 1
+damp_kappa2     = 0.0       # Constraint damping factor 2
+
+<problem>
+pgen_name = z4c_gauge_wave # problem generator name
+amp = 0.01
+
+<output1>
+file_type   = tab       # legacy VTK output
+variable    = z4c   # variables to be output
+dt          = 1      # time increment between outputs
+slice_x2    = 0.5
+slice_x3    = 0.5
+ghost_zones = true     # switch to output ghost cells
+
+<output2>
+file_type   = hst       # history data dump
+user_hist   = false
+data_format = %12.5e    # Optional data format string
+dt          = 1       # time increment between outputs
+
+#<output3>
+#file_type   = tab       # legacy VTK output
+#variable    = z4c   # variables to be output
+#dt          = 0.025      # time increment between outputs
+#slice_x2    = 0.5
+#slice_x3    = 0.5
+#ghost_zones = true     # switch to output ghost cells
diff --git a/inputs/z4c/awa/z4c_stability.athinput b/inputs/z4c/awa/z4c_stability.athinput
new file mode 100644
index 00000000..853cad95
--- /dev/null
+++ b/inputs/z4c/awa/z4c_stability.athinput
@@ -0,0 +1,76 @@
+# AthenaXXX input file for Z4c linear wave tests
+
+<comment>
+problem   = z4c linear waves
+reference = e.g. Daverio et al. arxiv:1810.12346 (2018)
+
+<job>
+basename  = stability    # problem ID: basename of output filenames
+
+<mesh>
+nghost    = 2          # Number of ghost cells
+nx1       = 4         # Number of zones in X1-direction
+x1min     = 0.0        # minimum value of X1
+x1max     = 1.0        # maximum value of X1
+ix1_bc    = periodic   # inner-X1 boundary flag
+ox1_bc    = periodic   # outer-X1 boundary flag
+
+nx2       = 4         # Number of zones in X2-direction
+x2min     = 0.0        # minimum value of X2
+x2max     = 1.0        # maximum value of X2
+ix2_bc    = periodic   # inner-X2 boundary flag
+ox2_bc    = periodic   # outer-X2 boundary flag
+
+nx3       = 50         # Number of zones in X3-direction
+x3min     = 0.0        # minimum value of X3
+x3max     = 1.0        # maximum value of X3
+ix3_bc    = periodic   # inner-X3 boundary flag
+ox3_bc    = periodic   # outer-X3 boundary flag
+
+<meshblock>
+nx1       = 4         # Number of cells in each MeshBlock, X1-dir
+nx2       = 4         # Number of cells in each MeshBlock, X2-dir
+nx3       = 50         # Number of cells in each MeshBlock, X3-dir
+
+<time>
+evolution  = dynamic   # dynamic/kinematic/static
+integrator = rk4       # time integration algorithm
+cfl_number = 0.5       # The Courant, Friedrichs, & Lewy (CFL) Number
+nlim       = -1        # cycle limit (no limit if <0)
+tlim       = 1000.0       # time limit
+ndiag      = 1         # cycles between diagostic output
+
+<z4c>
+diss       = 0.02
+lapse_harmonic  = 0.0                 # Harmonic lapse parameter mu_L
+lapse_oplog     = 1.0                 # 1+log lapse parameter
+shift_eta       = 2.0                 # Shift damping term
+chi_div_floor   = 0.00001
+damp_kappa1     = 0.00      # Constraint damping factor 1
+damp_kappa2     = 0.02
+
+
+<problem>
+pgen_name = z4c_stability # problem generator name
+user_hist = true      # enroll user-defined history function
+rho = 1
+
+<output1>
+file_type   = tab       # legacy VTK output
+variable    = z4c   # variables to be output
+dt          = 100      # time increment between outputs
+slice_x1    = 0.5
+slice_x2    = 0.5
+ghost_zones = true     # switch to output ghost cells
+
+<output2>
+file_type   = hst       # history data dump
+user_hist   = true
+data_format = %12.5e    # Optional data format string
+dt          = 1       # time increment between outputs
+
+<output2>
+file_type   = hst       # history data dump
+user_hist   = false
+data_format = %12.5e    # Optional data format string
+dt          = 1       # time increment between outputs
diff --git a/src/outputs/history.cpp b/src/outputs/history.cpp
index ea96bfca..fe77a589 100644
--- a/src/outputs/history.cpp
+++ b/src/outputs/history.cpp
@@ -161,14 +161,14 @@ void HistoryOutput::LoadHydroHistoryData(HistoryData *pdata, Mesh *pm) {
 void HistoryOutput::LoadZ4cHistoryData(HistoryData *pdata, Mesh *pm) {
   // set number of and names of history variables for z4c
   pdata->nhist = 8;
-  pdata->label[0] = "H-norm2";
-  pdata->label[1] = "M-norm2";
-  pdata->label[2] = "Mx-norm2";
-  pdata->label[3] = "My-norm2";
-  pdata->label[4] = "Mz-norm2";
-  pdata->label[5] = "Z-norm2";
-  pdata->label[6] = "Theta-norm2";
-  pdata->label[7] = "C-norm2";
+  pdata->label[0] = "C-norm2";
+  pdata->label[1] = "H-norm2";
+  pdata->label[2] = "M-norm2";
+  pdata->label[3] = "Z-norm2";
+  pdata->label[4] = "Mx-norm2";
+  pdata->label[5] = "My-norm2";
+  pdata->label[6] = "Mz-norm2";
+  pdata->label[7] = "Theta-norm2";
 
   // capture class variabels for kernel
   auto &u0_ = pm->pmb_pack->pz4c->u0;
@@ -201,14 +201,14 @@ void HistoryOutput::LoadZ4cHistoryData(HistoryData *pdata, Mesh *pm) {
 
     // Hydro conserved variables:
     array_sum::GlobalSum hvars;
-    hvars.the_array[0] = vol*SQR(u_con_(m,0,k,j,i)); // ||H||^2
-    hvars.the_array[1] = vol*u_con_(m,1,k,j,i);      // ||M||^2 (comes already squared)
-    hvars.the_array[2] = vol*SQR(u_con_(m,2,k,j,i)); // ||Mx||^2
-    hvars.the_array[3] = vol*SQR(u_con_(m,3,k,j,i)); // ||My||^2
-    hvars.the_array[4] = vol*u_con_(m,4,k,j,i);      // ||Mz||^2
-    hvars.the_array[5] = vol*u_con_(m,5,k,j,i);      // ||Z||^2 (comes already squared)
-    hvars.the_array[6] = vol*SQR(u0_(m,I_Z4c_Theta_,k,j,i)); // ||Theta||^2
-    hvars.the_array[7] = vol*u_con_(m,6,k,j,i);      // ||C||^2 (comes already squared)
+    hvars.the_array[0] = vol*u_con_(m,0,k,j,i); // ||C||^2 (comes already squared)
+    hvars.the_array[1] = vol*SQR(u_con_(m,1,k,j,i)); //||H||^2
+    hvars.the_array[2] = vol*u_con_(m,2,k,j,i); // ||M||^2 (comes already squared)
+    hvars.the_array[3] = vol*u_con_(m,3,k,j,i); // ||Z||^2 (comes already squared)
+    hvars.the_array[4] = vol*SQR(u_con_(m,4,k,j,i));      // ||Mx||^2
+    hvars.the_array[5] = vol*SQR(u_con_(m,5,k,j,i));      // ||My||^2
+    hvars.the_array[6] = vol*SQR(u_con_(m,6,k,j,i));      // ||Mz||^2
+    hvars.the_array[7] = vol*SQR(u0_(m,I_Z4c_Theta_,k,j,i)); // ||Theta||^2
 
     // fill rest of the_array with zeros, if nhist < NHISTORY_VARIABLES
     for (int n=nhist_; n<NHISTORY_VARIABLES; ++n) {
diff --git a/src/pgen/z4c_gauge_wave.cpp b/src/pgen/z4c_gauge_wave.cpp
new file mode 100644
index 00000000..e3c7e294
--- /dev/null
+++ b/src/pgen/z4c_gauge_wave.cpp
@@ -0,0 +1,119 @@
+//========================================================================================
+// AthenaXXX astrophysical plasma code
+// Copyright(C) 2020 James M. Stone <jmstone@ias.edu> and the Athena code team
+// Licensed under the 3-clause BSD License (the "LICENSE")
+//========================================================================================
+//! \file z4c_gauge_wave.cpp
+//! \brief z4c gauge wave test
+
+
+// C/C++ headers
+#include <algorithm> // min, max
+#include <cmath>     // sqrt()
+#include <cstdio>    // fopen(), fprintf(), freopen()
+#include <iostream>  // endl
+#include <sstream>   // stringstream
+#include <string>    // c_str()
+#include <limits>
+
+// Athena++ headers
+#include "athena.hpp"
+#include "globals.hpp"
+#include "parameter_input.hpp"
+#include "coordinates/cell_locations.hpp"
+#include "coordinates/adm.hpp"
+#include "mesh/mesh.hpp"
+#include "z4c/z4c.hpp"
+#include "driver/driver.hpp"
+#include "pgen/pgen.hpp"
+
+// function to compute errors in solution at end of run
+void Z4cGaugeWaveErrors(ParameterInput *pin, Mesh *pm);
+
+//----------------------------------------------------------------------------------------
+//! \fn void ProblemGenerator::UserProblem()
+//! \brief Sets initial conditions for gw linear wave tests
+
+void ProblemGenerator::UserProblem(ParameterInput *pin, const bool restart) {
+  if (restart)
+    return;
+
+  MeshBlockPack *pmbp = pmy_mesh_->pmb_pack;
+  if (pmbp->pz4c == nullptr) {
+    std::cout << "### FATAL ERROR in " << __FILE__ << " at line " << __LINE__ << std::endl
+              << "Z4c Wave test can only be run in Z4c, but no <z4c> block "
+              << "in input file" << std::endl;
+    exit(EXIT_FAILURE);
+  }
+
+  // Prepare Initial Data
+
+  // capture variables for the kernel
+  auto &indcs = pmbp->pmesh->mb_indcs;
+  auto &size = pmbp->pmb->mb_size;
+  int &is = indcs.is;
+  int &ie = indcs.ie;
+  int &js = indcs.js;
+  int &je = indcs.je;
+  int &ks = indcs.ks;
+  int &ke = indcs.ke;
+
+  int isg = is-indcs.ng; int ieg = ie+indcs.ng;
+  int jsg = js-indcs.ng; int jeg = je+indcs.ng;
+  int ksg = ks-indcs.ng; int keg = ke+indcs.ng;
+
+  auto &pz4c = pmbp->pz4c;
+  auto &adm = pmbp->padm->adm;
+  Real x1size = pmy_mesh_->mesh_size.x1max - pmy_mesh_->mesh_size.x1min;
+  std::cout << x1size << std::endl;
+
+  // Wave amplitude
+  Real amp = pin->GetOrAddReal("problem", "amp", 0.001);
+
+  par_for("pgen_gauge_wave", DevExeSpace(), 0, (pmbp->nmb_thispack - 1),
+      ksg, keg, jsg, jeg, isg, ieg, KOKKOS_LAMBDA(int m, int k, int j, int i) {
+    Real &x1min = size.d_view(m).x1min;
+    Real &x1max = size.d_view(m).x1max;
+
+    int nx1 = indcs.nx1;
+
+    Real x1v = CellCenterX(i - is, nx1, x1min, x1max);
+    Real H = amp * sin(2 * M_PI * x1v / x1size);
+    Real dH_dt = - amp * 2 * M_PI / x1size *
+                    cos(2 * M_PI * x1v / x1size);
+
+    adm.g_dd(m,0,0,k,j,i) = 1 - H;
+    adm.g_dd(m,0,1,k,j,i) = 0;
+    adm.g_dd(m,0,2,k,j,i) = 0;
+    adm.g_dd(m,1,1,k,j,i) = 1;
+    adm.g_dd(m,1,2,k,j,i) = 0;
+    adm.g_dd(m,2,2,k,j,i) = 1;
+
+    adm.vK_dd(m,0,0,k,j,i) = 0.5 * dH_dt / sqrt(adm.g_dd(m,0,0,k,j,i));
+    adm.vK_dd(m,0,1,k,j,i) = 0;
+    adm.vK_dd(m,0,2,k,j,i) = 0;
+    adm.vK_dd(m,1,1,k,j,i) = 0;
+    adm.vK_dd(m,1,2,k,j,i) = 0;
+    adm.vK_dd(m,2,2,k,j,i) = 0;
+
+    adm.alpha(m,k,j,i) = sqrt(1 - H);
+  });
+  switch (indcs.ng) {
+    case 2: pmbp->pz4c->ADMToZ4c<2>(pmbp, pin);
+            break;
+    case 3: pmbp->pz4c->ADMToZ4c<3>(pmbp, pin);
+            break;
+    case 4: pmbp->pz4c->ADMToZ4c<4>(pmbp, pin);
+            break;
+  }
+  pmbp->pz4c->Z4cToADM(pmbp);
+  switch (indcs.ng) {
+    case 2: pmbp->pz4c->ADMConstraints<2>(pmbp);
+            break;
+    case 3: pmbp->pz4c->ADMConstraints<3>(pmbp);
+            break;
+    case 4: pmbp->pz4c->ADMConstraints<4>(pmbp);
+            break;
+  }
+  return;
+}
diff --git a/src/pgen/z4c_stability.cpp b/src/pgen/z4c_stability.cpp
new file mode 100644
index 00000000..966e66a4
--- /dev/null
+++ b/src/pgen/z4c_stability.cpp
@@ -0,0 +1,242 @@
+//========================================================================================
+// AthenaXXX astrophysical plasma code
+// Copyright(C) 2020 James M. Stone <jmstone@ias.edu> and the Athena code team
+// Licensed under the 3-clause BSD License (the "LICENSE")
+//========================================================================================
+//! \file z4c_stability.cpp
+//! \brief z4c stability test
+
+
+// C/C++ headers
+#include <algorithm> // min, max
+#include <cmath>     // sqrt()
+#include <cstdio>    // fopen(), fprintf(), freopen()
+#include <iostream>  // endl
+#include <sstream>   // stringstream
+#include <string>    // c_str()
+#include <limits>
+
+// Athena++ headers
+#include "athena.hpp"
+#include "globals.hpp"
+#include "parameter_input.hpp"
+#include "coordinates/cell_locations.hpp"
+#include "coordinates/adm.hpp"
+#include "mesh/mesh.hpp"
+#include "z4c/z4c.hpp"
+#include "driver/driver.hpp"
+#include "pgen/pgen.hpp"
+
+#include <Kokkos_Random.hpp>
+
+
+// history function to output error as a function of time
+void Z4cStabilityErrors(HistoryData *pdata, Mesh *pm);
+
+//----------------------------------------------------------------------------------------
+//! \fn void ProblemGenerator::UserProblem()
+//! \brief Sets initial conditions for robust stability test
+
+void ProblemGenerator::UserProblem(ParameterInput *pin, const bool restart) {
+  user_hist_func = &Z4cStabilityErrors;
+  if (restart)
+    return;
+
+  MeshBlockPack *pmbp = pmy_mesh_->pmb_pack;
+  if (pmbp->pz4c == nullptr) {
+    std::cout << "### FATAL ERROR in " << __FILE__ << " at line " << __LINE__ << std::endl
+              << "Z4c stability test can only be run in Z4c, but no <z4c> block "
+              << "in input file" << std::endl;
+    exit(EXIT_FAILURE);
+  }
+
+  // Prepare Initial Data
+
+  // capture variables for the kernel
+  auto &indcs = pmbp->pmesh->mb_indcs;
+  int &is = indcs.is; int &ie = indcs.ie;
+  int &js = indcs.js; int &je = indcs.je;
+  int &ks = indcs.ks; int &ke = indcs.ke;
+  // For GLOOPS
+  int ksg = ks-indcs.ng; int keg = ke+indcs.ng;
+  int nmb = pmbp->nmb_thispack;
+  int ncell1 = 2*indcs.ng + indcs.nx1;
+  int ncell2 = 2*indcs.ng + indcs.nx2;
+  int ncell12 = ncell1 * ncell2;
+
+  auto &adm = pmbp->padm->adm;
+  auto &z4c = pmbp->pz4c->z4c;
+
+  // perturbation amplitude 1e-10 dx^2
+  Real rho = pin->GetOrAddReal("problem", "rho", 1);
+  Real amp = 1e-10/rho/rho;
+
+  // random number
+  Kokkos::Random_XorShift64_Pool<> rand_pool64(pmbp->gids);
+  size_t scr_size = ScrArray1D<Real>::shmem_size(12);
+  int scr_level = 0;
+  par_for_outer("pgen_stability",DevExeSpace(), scr_size, scr_level, 0,nmb-1,ksg,keg,
+  KOKKOS_LAMBDA(TeamMember_t member, const int m, const int k) {
+    // get random number state this thread
+    auto rand_gen = rand_pool64.get_state();
+
+    // scratch array holding random numbers for 12 variables.
+    // only perturbaion in the z direction is added
+    // in the x and y directions the values of perturbation
+    // are always the same
+    ScrArray1D<Real> rnum(member.team_scratch(scr_level), 12);
+    for(int a = 0; a < 12; ++a) {
+      rnum(a) = amp*2*(rand_gen.drand() - 0.5);
+    }
+
+    par_for_inner(member, 0, ncell12-1, [&](const int ij) {
+      const int j = ij/ncell1;
+      const int i = ij%ncell1;
+
+      // Minkowski spacetime
+      for(int a = 0; a < 3; ++a)
+      for(int b = a; b < 3; ++b) {
+        adm.g_dd(m,a,b,k,j,i) = (a == b ? 1. : 0.);
+      }
+      adm.alpha(m,k,j,i) = 1.0;
+
+      // g_ab perturbation
+      int c = 0;
+      for(int a = 0; a < 3; ++a)
+      for(int b = a; b < 3; ++b) {
+        adm.g_dd(m,a,b,k,j,i) += rnum(c);
+        c++;
+      }
+
+      // K_ab perturbation
+      for(int a = 0; a < 3; ++a)
+      for(int b = a; b < 3; ++b) {
+        adm.vK_dd(m,a,b,k,j,i) += rnum(c);
+        c++;
+      }
+
+      rand_pool64.free_state(rand_gen);  // free state for use by other threads
+
+      z4c.alpha(m,k,j,i) = 1;
+      z4c.chi(m,k,j,i) = 1;
+      z4c.vKhat(m,k,j,i) = 0;
+      z4c.vTheta(m,k,j,i) = 0;
+    });
+  });
+
+  switch (indcs.ng) {
+    case 2: pmbp->pz4c->ADMToZ4c<2>(pmbp, pin);
+            break;
+    case 3: pmbp->pz4c->ADMToZ4c<3>(pmbp, pin);
+            break;
+    case 4: pmbp->pz4c->ADMToZ4c<4>(pmbp, pin);
+            break;
+  }
+  pmbp->pz4c->Z4cToADM(pmbp);
+  switch (indcs.ng) {
+    case 2: pmbp->pz4c->ADMConstraints<2>(pmbp);
+            break;
+    case 3: pmbp->pz4c->ADMConstraints<3>(pmbp);
+            break;
+    case 4: pmbp->pz4c->ADMConstraints<4>(pmbp);
+            break;
+  }
+  return;
+}
+
+void Z4cStabilityErrors(HistoryData *pdata, Mesh *pm) {
+  pdata->nhist = 12;
+  pdata->label[0] = "LINF-Err";
+  pdata->label[1] = "RMS-Err";
+
+  Real l1_err[6];
+  Real linfty_err=0.0;
+  int nvars=0;
+
+  // capture class variables for kernel
+  auto &indcs = pm->mb_indcs;
+  int &nx1 = indcs.nx1;
+  int &nx2 = indcs.nx2;
+  int &nx3 = indcs.nx3;
+  int &is = indcs.is;
+  int &js = indcs.js;
+  int &ks = indcs.ks;
+  MeshBlockPack *pmbp = pm->pmb_pack;
+  auto &size = pmbp->pmb->mb_size;
+
+  // compute errors for Hydro  -----------------------------------------------------------
+  if (pmbp->pz4c != nullptr) {
+    nvars = 6; // 6 metric components
+    auto &pz4c = pmbp->pz4c;
+    auto &u0_ = pmbp->pz4c->u0;
+
+    const int nmkji = (pmbp->nmb_thispack)*nx3*nx2*nx1;
+    const int nkji = nx3*nx2*nx1;
+    const int nji  = nx2*nx1;
+    array_sum::GlobalSum sum_this_mb;
+    Kokkos::parallel_reduce("LW-err",Kokkos::RangePolicy<>(DevExeSpace(), 0, nmkji),
+    KOKKOS_LAMBDA(const int &idx, array_sum::GlobalSum &mb_sum, Real &max_err) {
+      // compute n,k,j,i indices of thread
+      int m = (idx)/nkji;
+      int k = (idx - m*nkji)/nji;
+      int j = (idx - m*nkji - k*nji)/nx1;
+      int i = (idx - m*nkji - k*nji - j*nx1) + is;
+      k += ks;
+      j += js;
+
+      Real vol = size.d_view(m).dx1*size.d_view(m).dx2*size.d_view(m).dx3;
+
+      // g_ij's:
+      array_sum::GlobalSum evars;
+      evars.the_array[0] = vol*fabs(1.
+                          - u0_(m,pz4c->I_Z4C_GXX,k,j,i));
+      max_err = fmax(max_err, evars.the_array[0]);
+      evars.the_array[1] = vol*fabs(0.
+                          - u0_(m,pz4c->I_Z4C_GXY,k,j,i));
+      max_err = fmax(max_err, evars.the_array[1]);
+      evars.the_array[2] = vol*fabs(0.
+                          - u0_(m,pz4c->I_Z4C_GXZ,k,j,i));
+      max_err = fmax(max_err, evars.the_array[2]);
+      evars.the_array[3] = vol*fabs(1.
+                          - u0_(m,pz4c->I_Z4C_GYY,k,j,i));
+      max_err = fmax(max_err, evars.the_array[3]);
+      evars.the_array[4] = vol*fabs(0.
+                          - u0_(m,pz4c->I_Z4C_GYZ,k,j,i));
+      max_err = fmax(max_err, evars.the_array[4]);
+      evars.the_array[5] = vol*fabs(1.
+                          - u0_(m,pz4c->I_Z4C_GZZ,k,j,i));
+      max_err = fmax(max_err, evars.the_array[5]);
+
+      // fill rest of the_array with zeros, if narray < NREDUCTION_VARIABLES
+      for (int n=nvars; n<NREDUCTION_VARIABLES; ++n) {
+        evars.the_array[n] = 0.0;
+      }
+
+      // sum into parallel reduce
+      mb_sum += evars;
+    }, Kokkos::Sum<array_sum::GlobalSum>(sum_this_mb), Kokkos::Max<Real>(linfty_err));
+
+    // store data into l1_err array
+    for (int n=0; n<nvars; ++n) {
+      l1_err[n] = sum_this_mb.the_array[n];
+    }
+  }
+
+  // normalize errors by number of cells
+  Real vol=  (pmbp->pmesh->mesh_size.x1max - pmbp->pmesh->mesh_size.x1min)
+            *(pmbp->pmesh->mesh_size.x2max - pmbp->pmesh->mesh_size.x2min)
+            *(pmbp->pmesh->mesh_size.x3max - pmbp->pmesh->mesh_size.x3min);
+  for (int i=0; i<nvars; ++i) l1_err[i] = l1_err[i]/vol;
+  linfty_err /= vol;
+
+  // compute rms error
+  Real rms_err = 0.0;
+  for (int i=0; i<nvars; ++i) {
+    rms_err += SQR(l1_err[i]);
+  }
+  rms_err = std::sqrt(rms_err);
+
+  pdata->hdata[0] = linfty_err;
+  pdata->hdata[1] = rms_err;
+  return;
+}