// CCTK_REAL maxabs = 0.0;
      // for (auto &restrict leveldata : ghext->leveldata) {
      //   auto &restrict groupdata = leveldata.groupdata.at(gi);
      //   MultiFab &mfab = *groupdata.mfab.at(tl);
      //   maxabs = fmax(maxabs, mfab.norminf(vi));
      // }
      // CCTK_VINFO("            %g", maxabs);

#include <AMReX_FillPatchUtil.H>
#include <AMReX_Interpolater.H>

#include <AMReX_PhysBCFunct.H>

  for (int level = 0; level <numlevels; ++level) {

  for (int level = 0; level < numlevels; ++level) {

  // TODO: Synchronize only current time level. This probably requires
  // setting up different mfabs for each time level.

      for (int tl = 0; tl < sync_tl; ++tl)
        groupdata.mfab.at(tl)->FillBoundary(
            ghext->amrmesh->Geom(leveldata.level).periodicity());

      if (leveldata.level == 0) {
        // Coarsest level: Copy from adjacent boxes on same level
        for (int tl = 0; tl < sync_tl; ++tl)
          groupdata.mfab.at(tl)->FillBoundary(
              ghext->amrmesh->Geom(leveldata.level).periodicity());
      } else {
        // Refined level: Prolongate from next coarser level, and then
        // copy from adjacent boxes on same level
        const int level = leveldata.level;
        auto &restrict coarsegroupdata =
            ghext->leveldata.at(level - 1).groupdata.at(gi);
        assert(coarsegroupdata.numvars == groupdata.numvars);
        PhysBCFunctNoOp cphysbc;
        PhysBCFunctNoOp fphysbc;
        const IntVect reffact{2, 2, 2};
        CellBilinear interp;
        // periodic boundaries
        const BCRec bcrec(BCType::int_dir, BCType::int_dir, BCType::int_dir,
                          BCType::int_dir, BCType::int_dir, BCType::int_dir);
        const Vector<BCRec> bcs(groupdata.numvars, bcrec);
        for (int tl = 0; tl < sync_tl; ++tl)
          InterpFromCoarseLevel(
              *groupdata.mfab.at(tl), 0.0, *coarsegroupdata.mfab.at(tl), 0, 0,
              groupdata.numvars, ghext->amrmesh->Geom(level - 1),
              ghext->amrmesh->Geom(level), cphysbc, 0, fphysbc, 0, reffact,
              &interp, bcs, 0);
      }