READS: Coordinates::coords(everywhere)

    READS: Coordinates::vertex_coords(everywhere)

      CCTK_VarIndex("Coordinates::coordx"),
      CCTK_VarIndex("Coordinates::coordy"),
      CCTK_VarIndex("Coordinates::coordz"),

      CCTK_VarIndex("Coordinates::vcoordx"),
      CCTK_VarIndex("Coordinates::vcoordy"),
      CCTK_VarIndex("Coordinates::vcoordz"),

CCTK_REAL vertex_coords TYPE=gf TAGS='index={0 0 0}' { vcoordx vcoordy vcoordz } "Vertex coordinates"

CCTK_REAL coords TYPE=gf TAGS='index={0 0 0}' { coordx coordy coordz } "Coordinates"

CCTK_REAL cell_coords TYPE=gf TAGS='index={1 1 1}' { ccoordx ccoordy ccoordz } "Cell centre coordinates"
CCTK_REAL cell_volume TYPE=gf TAGS='index={1 1 1}' { cvol } "Cell volume"

  WRITES: coords(everywhere)

  WRITES: vertex_coords(everywhere)
  WRITES: cell_coords(everywhere)
  WRITES: cell_volume(everywhere)

  const Loop::GF3D<CCTK_REAL, 0, 0, 0> vcoordx_(cctkGH, vcoordx);
  const Loop::GF3D<CCTK_REAL, 0, 0, 0> vcoordy_(cctkGH, vcoordy);
  const Loop::GF3D<CCTK_REAL, 0, 0, 0> vcoordz_(cctkGH, vcoordz);
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> ccoordx_(cctkGH, ccoordx);
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> ccoordy_(cctkGH, ccoordy);
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> ccoordz_(cctkGH, ccoordz);

  const Loop::GF3D<CCTK_REAL, 0, 0, 0> coordx_(cctkGH, coordx);
  const Loop::GF3D<CCTK_REAL, 0, 0, 0> coordy_(cctkGH, coordy);
  const Loop::GF3D<CCTK_REAL, 0, 0, 0> coordz_(cctkGH, coordz);

  const Loop::GF3D<CCTK_REAL, 1, 1, 1> cvol_(cctkGH, cvol);

      cctkGH, [&](const Loop::PointDesc &p) { coordx_(p.I) = p.x; });

      cctkGH, [&](const Loop::PointDesc &p) { vcoordx_(p.I) = p.x; });

      cctkGH, [&](const Loop::PointDesc &p) { coordy_(p.I) = p.y; });

      cctkGH, [&](const Loop::PointDesc &p) { vcoordy_(p.I) = p.y; });

      cctkGH, [&](const Loop::PointDesc &p) { coordz_(p.I) = p.z; });

      cctkGH, [&](const Loop::PointDesc &p) { vcoordz_(p.I) = p.z; });
  Loop::loop_all<1, 1, 1>(
      cctkGH, [&](const Loop::PointDesc &p) { ccoordx_(p.I) = p.x; });
  Loop::loop_all<1, 1, 1>(
      cctkGH, [&](const Loop::PointDesc &p) { ccoordy_(p.I) = p.y; });
  Loop::loop_all<1, 1, 1>(
      cctkGH, [&](const Loop::PointDesc &p) { ccoordz_(p.I) = p.z; });
  Loop::loop_all<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
    cvol_(p.I) = p.dx * p.dy * p.dz;
  });

    READS: Coordinates::coords(everywhere)

    READS: Coordinates::vertex_coords(everywhere)

        xx = coordx[ind] - center_offset[0];
        yy = coordy[ind] - center_offset[1];
        zz = coordz[ind] - center_offset[2];

        xx = vcoordx[ind] - center_offset[0];
        yy = vcoordy[ind] - center_offset[1];
        zz = vcoordz[ind] - center_offset[2];

                    coordx, coordy, coordz,

                    vcoordx, vcoordy, vcoordz,