diff --git a/src/exo3/cs_find_block_id.cpp b/src/exo3/cs_find_block_id.cpp
index 5f7a58e2..e5f7d45b 100644
--- a/src/exo3/cs_find_block_id.cpp
+++ b/src/exo3/cs_find_block_id.cpp
@@ -9,8 +9,9 @@
 #include "cubed_sphere.hpp"
 
 int CubedSphere::FindBlockID(LogicalLocation const& loc) {
-  int lv2_lx2 = loc.lx2 >> (loc.level - 2);
-  int lv2_lx3 = loc.lx3 >> (loc.level - 2);
+
+  int lv2_lx2, lv2_lx3, local_lx2, local_lx3, bound_lim;
+  GetLocalIndex(&lv2_lx2, &lv2_lx3, &local_lx2, &local_lx3, &bound_lim, loc);
 
   // Determine the block number
   int block_id;
@@ -27,6 +28,7 @@ int CubedSphere::FindBlockID(LogicalLocation const& loc) {
           std::stringstream msg;
           msg << "Error: something wrong, check the geometry setup of the "
                  "cubed sphere. \n";
+          msg << "lv2_lx2: " << lv2_lx2 << " lv2_lx3: " << lv2_lx3 << "bound_lim: " << bound_lim << "loc.lx2: " << loc.lx2 << "loc.lx3: " << loc.lx3 << std::endl;
           msg << "----------------------------------" << std::endl;
           ATHENA_ERROR(msg);
       }
@@ -73,3 +75,24 @@ int CubedSphere::FindBlockID(LogicalLocation const& loc) {
 
   return block_id;
 }
+
+void CubedSphere::GetLocalIndex(int *lv2_lx2, int *lv2_lx3, int *local_lx2, int *local_lx3, int *bound_lim, LogicalLocation const& loc) {
+#ifdef USE_NBLOCKS
+  // Updated method, need to manually setup in configure.hpp, allow 6*n^2 blocks
+  *bound_lim = (int)(sqrt(NBLOCKS / 6) - 0.5);
+  // Find relative location within block
+  
+  *lv2_lx2 = loc.lx2 / (*bound_lim + 1);
+  *lv2_lx3 = loc.lx3 / (*bound_lim + 1);
+  *local_lx2 = loc.lx2 - (*lv2_lx2 * (*bound_lim + 1));
+  *local_lx3 = loc.lx3 - (*lv2_lx3 * (*bound_lim + 1));
+#else
+  // Old method, suitable for 6*4^n blocks
+  *bound_lim = (1 << (loc.level - 2)) - 1;
+    // Find relative location within block
+  *lv2_lx2 = loc.lx2 >> (loc.level - 2);
+  *lv2_lx3 = loc.lx3 >> (loc.level - 2);
+  *local_lx2 = loc.lx2 - (*lv2_lx2 << (loc.level - 2));
+  *local_lx3 = loc.lx3 - (*lv2_lx3 << (loc.level - 2));
+#endif
+}
\ No newline at end of file
diff --git a/src/exo3/cs_transform_ox.cpp b/src/exo3/cs_transform_ox.cpp
index 477b1059..0154f464 100644
--- a/src/exo3/cs_transform_ox.cpp
+++ b/src/exo3/cs_transform_ox.cpp
@@ -12,16 +12,12 @@ void CubedSphere::TransformOX(int *ox2, int *ox3, int *tox2, int *tox3,
                               LogicalLocation const &loc) {
   // Find the block ID
   int block_id = FindBlockID(loc);
+  int lv2_lx2, lv2_lx3, local_lx2, local_lx3, bound_lim;
+  GetLocalIndex(&lv2_lx2, &lv2_lx3, &local_lx2, &local_lx3, &bound_lim, loc);
 
-  // Find relative location within block
-  int lv2_lx2 = loc.lx2 >> (loc.level - 2);
-  int lv2_lx3 = loc.lx3 >> (loc.level - 2);
-  int local_lx2 = loc.lx2 - (lv2_lx2 << (loc.level - 2));
-  int local_lx3 = loc.lx3 - (lv2_lx3 << (loc.level - 2));
-  int bound_lim = (1 << (loc.level - 2)) - 1;
 
-  // Hard code the cases...
-  // No need to consider the corner cases, abandon in reading buffers.
+  // Hard code the cases... 
+  // The corner cases are not used, abandoned after reading buffers.
   int target_block = -1;           // Block id of target
   int target_loc_2, target_loc_3;  // local x2 and x3 in target block
 
@@ -168,6 +164,7 @@ void CubedSphere::TransformOX(int *ox2, int *ox3, int *tox2, int *tox3,
       std::stringstream msg;
       msg << "Error: something wrong, check the geometry setup of the cubed "
              "sphere. \n";
+      msg << "Block ID: " << block_id << std::endl;
       msg << "----------------------------------" << std::endl;
       ATHENA_ERROR(msg);
   }
@@ -210,8 +207,13 @@ void CubedSphere::TransformOX(int *ox2, int *ox3, int *tox2, int *tox3,
         msg << "----------------------------------" << std::endl;
         ATHENA_ERROR(msg);
     }
+#ifdef USE_NBLOCKS
+    lx3_0 = (lx3_0 * (bound_lim + 1));
+    lx2_0 = (lx2_0 * (bound_lim + 1));
+#else
     lx3_0 = (lx3_0 << (loc.level - 2));
     lx2_0 = (lx2_0 << (loc.level - 2));
+#endif
     // Add up first block and local positions
     int lx3_t = lx3_0 + target_loc_3;
     int lx2_t = lx2_0 + target_loc_2;
diff --git a/src/exo3/cubed_sphere.cpp b/src/exo3/cubed_sphere.cpp
index b9a6933e..90795f80 100644
--- a/src/exo3/cubed_sphere.cpp
+++ b/src/exo3/cubed_sphere.cpp
@@ -39,9 +39,10 @@ CubedSphere::CubedSphere(MeshBlock *pmb) : pmy_block_(pmb) {
 
 Real CubedSphere::GenerateMeshX2(Real x, LogicalLocation const &loc) {
   Real x_l, x_u;
-  int lx2_lv2 = loc.lx2 >> (loc.level - 2);
+  int lv2_lx2, lv2_lx3, local_lx2, local_lx3, bound_lim;
+  GetLocalIndex(&lv2_lx2, &lv2_lx3, &local_lx2, &local_lx3, &bound_lim, loc);
 
-  switch (lx2_lv2) {
+  switch (lv2_lx2) {
     case 0:
       x_l = -0.5;
       x_u = 0.0;
@@ -56,9 +57,10 @@ Real CubedSphere::GenerateMeshX2(Real x, LogicalLocation const &loc) {
 
 Real CubedSphere::GenerateMeshX3(Real x, LogicalLocation const &loc) {
   Real x_l, x_u;
-  int lx3_lv2 = loc.lx3 >> (loc.level - 2);
+  int lv2_lx2, lv2_lx3, local_lx2, local_lx3, bound_lim;
+  GetLocalIndex(&lv2_lx2, &lv2_lx3, &local_lx2, &local_lx3, &bound_lim, loc);
 
-  switch (lx3_lv2) {
+  switch (lv2_lx3) {
     case 0:
       x_l = -0.5;
       x_u = -1.0 / 6.0;
@@ -81,10 +83,8 @@ Real CubedSphere::GenerateMeshX3(Real x, LogicalLocation const &loc) {
 void CubedSphere::GetLatLon(Real *lat, Real *lon, int k, int j, int i) const {
   auto pcoord = pmy_block_->pcoord;
   auto &loc = pmy_block_->loc;
+  int blockID = FindBlockID(loc);
 
-  int lv2_lx2 = loc.lx2 >> (loc.level - 2);
-  int lv2_lx3 = loc.lx3 >> (loc.level - 2);
-  int blockID = lv2_lx2 + lv2_lx3 * 2 + 1;
   // Calculate the needed parameters
   Real dX = tan(pcoord->x2v(j));
   Real dY = tan(pcoord->x3v(k));
@@ -99,10 +99,7 @@ void CubedSphere::GetLatLonFace2(Real *lat, Real *lon, int k, int j,
                                  int i) const {
   auto pcoord = pmy_block_->pcoord;
   auto &loc = pmy_block_->loc;
-
-  int lv2_lx2 = loc.lx2 >> (loc.level - 2);
-  int lv2_lx3 = loc.lx3 >> (loc.level - 2);
-  int blockID = lv2_lx2 + lv2_lx3 * 2 + 1;
+  int blockID = FindBlockID(loc);
 
   // Calculate the needed parameters
   Real dX = tan(pcoord->x2f(j));
@@ -118,10 +115,7 @@ void CubedSphere::GetLatLonFace3(Real *lat, Real *lon, int k, int j,
                                  int i) const {
   auto pcoord = pmy_block_->pcoord;
   auto &loc = pcoord->pmy_block->loc;
-
-  int lv2_lx2 = loc.lx2 >> (loc.level - 2);
-  int lv2_lx3 = loc.lx3 >> (loc.level - 2);
-  int blockID = lv2_lx2 + lv2_lx3 * 2 + 1;
+  int blockID = FindBlockID(loc);
 
   // Calculate the needed parameters
   Real dX = tan(pcoord->x2v(j));
@@ -136,10 +130,7 @@ void CubedSphere::GetUV(Real *U, Real *V, Real V2, Real V3, int k, int j,
                         int i) const {
   auto pcoord = pmy_block_->pcoord;
   auto &loc = pmy_block_->loc;
-
-  int lv2_lx2 = loc.lx2 >> (loc.level - 2);
-  int lv2_lx3 = loc.lx3 >> (loc.level - 2);
-  int blockID = lv2_lx2 + lv2_lx3 * 2 + 1;
+  int blockID = FindBlockID(loc);
 
   Real X = tan(pcoord->x2v(j));
   Real Y = tan(pcoord->x3v(k));
@@ -153,10 +144,7 @@ void CubedSphere::GetVyVz(Real *V2, Real *V3, Real U, Real V, int k, int j,
                           int i) const {
   auto pcoord = pmy_block_->pcoord;
   auto &loc = pmy_block_->loc;
-
-  int lv2_lx2 = loc.lx2 >> (loc.level - 2);
-  int lv2_lx3 = loc.lx3 >> (loc.level - 2);
-  int blockID = lv2_lx2 + lv2_lx3 * 2 + 1;
+  int blockID = FindBlockID(loc);
 
   // Calculate the needed parameters
   Real X = tan(pcoord->x2v(j));
@@ -169,10 +157,7 @@ void CubedSphere::CalculateCoriolisForce2(int i2, int i3, Real v2, Real v3,
                                           Real *cF3) const {
   auto pcoord = pmy_block_->pcoord;
   auto &loc = pmy_block_->loc;
-
-  int lv2_lx2 = loc.lx2 >> (loc.level - 2);
-  int lv2_lx3 = loc.lx3 >> (loc.level - 2);
-  int blockID = lv2_lx2 + lv2_lx3 * 2 + 1;
+  int blockID = FindBlockID(loc);
 
   Real x = tan(pcoord->x2v(i2));
   Real y = tan(pcoord->x3v(i3));
@@ -203,10 +188,7 @@ void CubedSphere::CalculateCoriolisForce3(int i2, int i3, Real v1, Real v2,
                                           Real *cF3) const {
   auto pcoord = pmy_block_->pcoord;
   auto &loc = pmy_block_->loc;
-
-  int lv2_lx2 = loc.lx2 >> (loc.level - 2);
-  int lv2_lx3 = loc.lx3 >> (loc.level - 2);
-  int blockID = lv2_lx2 + lv2_lx3 * 2 + 1;
+  int blockID = FindBlockID(loc);
 
   Real x = tan(pcoord->x2v(i2));
   Real y = tan(pcoord->x3v(i3));
@@ -333,14 +315,10 @@ void CubedSphere::sendNeighborBlocks(int ox2, int ox3, int tg_rank,
                                      int tg_gid) {
   MeshBlock *pmb = pmy_block_;
   auto &loc = pmb->loc;
-
-  int lv2_lx2 = loc.lx2 >> (loc.level - 2);
-  int lv2_lx3 = loc.lx3 >> (loc.level - 2);
-  int blockID = lv2_lx2 + lv2_lx3 * 2 + 1;
+  int lv2_lx2, lv2_lx3, local_lx2, local_lx3, bound_lim;
+  GetLocalIndex(&lv2_lx2, &lv2_lx3, &local_lx2, &local_lx3, &bound_lim, loc);
+  int blockID = FindBlockID(loc);
   // Calculate local ID
-  int local_lx2 = loc.lx2 - (lv2_lx2 << (loc.level - 2));
-  int local_lx3 = loc.lx3 - (lv2_lx3 << (loc.level - 2));
-  int bound_lim = (1 << (loc.level - 2)) - 1;
   int tox2, tox3;
   int DirTag, DirNum, ownTag;  // Tag for target, and the numbering of axis
   int ox2_bkp = ox2;
@@ -560,6 +538,16 @@ void CubedSphere::sendNeighborBlocks(int ox2, int ox3, int tg_rank,
   }
 
   // Calculate the tag of destination
+#ifdef USE_NBLOCKS
+  if (tox2 == -1)
+    DirTag = 0 + 4 * pmb->gid + 24 * (bound_lim + 1) * tg_gid;
+  if (tox2 == 1)
+    DirTag = 1 + 4 * pmb->gid + 24 * (bound_lim + 1) * tg_gid;
+  if (tox3 == -1)
+    DirTag = 2 + 4 * pmb->gid + 24 * (bound_lim + 1) * tg_gid;
+  if (tox3 == 1)
+    DirTag = 3 + 4 * pmb->gid + 24 * (bound_lim + 1) * tg_gid;
+#else
   if (tox2 == -1)
     DirTag = 0 + 4 * pmb->gid + 24 * (1 << (loc.level - 2)) * tg_gid;
   if (tox2 == 1)
@@ -568,6 +556,7 @@ void CubedSphere::sendNeighborBlocks(int ox2, int ox3, int tg_rank,
     DirTag = 2 + 4 * pmb->gid + 24 * (1 << (loc.level - 2)) * tg_gid;
   if (tox3 == 1)
     DirTag = 3 + 4 * pmb->gid + 24 * (1 << (loc.level - 2)) * tg_gid;
+#endif
   // Send by MPI: we don't care whether it is in the same process for now
   if (ox2 == -1) ownTag = 0;
   if (ox2 == 1) ownTag = 1;
@@ -586,13 +575,10 @@ void CubedSphere::recvNeighborBlocks(int ox2, int ox3, int tg_rank,
   MeshBlock *pmb = pmy_block_;
   auto &loc = pmb->loc;
 
-  int lv2_lx2 = loc.lx2 >> (loc.level - 2);
-  int lv2_lx3 = loc.lx3 >> (loc.level - 2);
+  int lv2_lx2, lv2_lx3, local_lx2, local_lx3, bound_lim;
+  GetLocalIndex(&lv2_lx2, &lv2_lx3, &local_lx2, &local_lx3, &bound_lim, loc);
   int blockID = FindBlockID(loc);
-  // Calculate local ID
-  int local_lx2 = loc.lx2 - (lv2_lx2 << (loc.level - 2));
-  int local_lx3 = loc.lx3 - (lv2_lx3 << (loc.level - 2));
-  int bound_lim = (1 << (loc.level - 2)) - 1;
+
   int tox2, tox3;
   int DirTag, DirNum, ownTag;  // Tag for receiving, and the numbering of axis
                                // to place the values
@@ -658,6 +644,16 @@ void CubedSphere::recvNeighborBlocks(int ox2, int ox3, int tg_rank,
   int dsize = ((kb2 - kb1 + 1) * (jb2 - jb1 + 1) * (ib2 - ib1 + 1) * NWAVE);
   Real *data = new Real[dsize];
   // Calculate the tag for receiving
+#ifdef USE_NBLOCKS
+  if (ox2 == -1)
+    DirTag = 0 + 4 * tg_gid + 24 * (bound_lim + 1) * pmb->gid;
+  if (ox2 == 1)
+    DirTag = 1 + 4 * tg_gid + 24 * (bound_lim + 1) * pmb->gid;
+  if (ox3 == -1)
+    DirTag = 2 + 4 * tg_gid + 24 * (bound_lim + 1) * pmb->gid;
+  if (ox3 == 1)
+    DirTag = 3 + 4 * tg_gid + 24 * (bound_lim + 1) * pmb->gid;
+#else
   if (ox2 == -1)
     DirTag = 0 + 4 * tg_gid + 24 * (1 << (loc.level - 2)) * pmb->gid;
   if (ox2 == 1)
@@ -666,6 +662,7 @@ void CubedSphere::recvNeighborBlocks(int ox2, int ox3, int tg_rank,
     DirTag = 2 + 4 * tg_gid + 24 * (1 << (loc.level - 2)) * pmb->gid;
   if (ox3 == 1)
     DirTag = 3 + 4 * tg_gid + 24 * (1 << (loc.level - 2)) * pmb->gid;
+#endif
 
 #ifdef MPI_PARALLEL
   MPI_Recv(data, dsize, MPI_DOUBLE, tg_rank, DirTag, MPI_COMM_WORLD,
diff --git a/src/exo3/cubed_sphere.hpp b/src/exo3/cubed_sphere.hpp
index 8995f23a..9c8649d1 100644
--- a/src/exo3/cubed_sphere.hpp
+++ b/src/exo3/cubed_sphere.hpp
@@ -17,6 +17,9 @@ class CubedSphere {
   ~CubedSphere() {}
 
   static int FindBlockID(LogicalLocation const &loc);
+  static void GetLocalIndex(int *lv2_lx2, int *lv2_lx3, int *local_lx2,
+                            int *local_lx3, int *bound_lim,
+                            LogicalLocation const &loc);
   static void TransformOX(int *ox2, int *ox3, int *tox2, int *tox3,
                           LogicalLocation const &loc);
 
@@ -66,7 +69,6 @@ class CubedSphere {
 #endif
 
   std::vector<Real> LRDataBuffer[4];
-
   MeshBlock *pmy_block_;
 };
 
diff --git a/src/exo3/cubed_sphere_utility.cpp b/src/exo3/cubed_sphere_utility.cpp
index b6180877..933f3231 100644
--- a/src/exo3/cubed_sphere_utility.cpp
+++ b/src/exo3/cubed_sphere_utility.cpp
@@ -102,10 +102,9 @@ void InteprolateX2(const AthenaArray<Real> &src, AthenaArray<Real> &tgt,
                    int sk, int ek, int TypeFlag) {
   // Interpolation along X2 (j) axis, used before sending data to X3 (k) axis
   // Get the local indices
-  int lv2_lx2 = loc.lx2 >> (loc.level - 2);
-  int lv2_lx3 = loc.lx3 >> (loc.level - 2);
-  int local_lx2 = loc.lx2 - (lv2_lx2 << (loc.level - 2));
-  int bound_lim = (1 << (loc.level - 2)) - 1;
+  int lv2_lx2, lv2_lx3, local_lx2, local_lx3, bound_lim;
+  CubedSphere::GetLocalIndex(&lv2_lx2, &lv2_lx3, &local_lx2, &local_lx3, &bound_lim, loc);
+
   int meshblock_size = ej - sj + 1;
   int N_blk = meshblock_size *
               (bound_lim + 1);  // N in X2 direction for each panel. This value
@@ -233,10 +232,9 @@ void InteprolateX3(const AthenaArray<Real> &src, AthenaArray<Real> &tgt,
                    int sk, int ek, int TypeFlag) {
   // Interpolation along X3 (k) axis, used before sending data to ghost zone in
   // X2 (j) direction Get the local indices
-  int lv2_lx2 = loc.lx2 >> (loc.level - 2);
-  int lv2_lx3 = loc.lx3 >> (loc.level - 2);
-  int local_lx3 = loc.lx3 - (lv2_lx3 << (loc.level - 2));
-  int bound_lim = (1 << (loc.level - 2)) - 1;
+  int lv2_lx2, lv2_lx3, local_lx2, local_lx3, bound_lim;
+  CubedSphere::GetLocalIndex(&lv2_lx2, &lv2_lx3, &local_lx2, &local_lx3, &bound_lim, loc);
+
   int meshblock_size = ek - sk + 1;
   int N_blk = meshblock_size *
               (bound_lim + 1);  // N in X2 direction for each panel. This value
@@ -403,11 +401,8 @@ void PackData(const AthenaArray<Real> &src, Real *buf, int sn, int en, int si,
   }
 
   // Get the local indices
-  int lv2_lx2 = loc.lx2 >> (loc.level - 2);
-  int lv2_lx3 = loc.lx3 >> (loc.level - 2);
-  int local_lx2 = loc.lx2 - (lv2_lx2 << (loc.level - 2));
-  int local_lx3 = loc.lx3 - (lv2_lx3 << (loc.level - 2));
-  int bound_lim = (1 << (loc.level - 2)) - 1;
+  int lv2_lx2, lv2_lx3, local_lx2, local_lx3, bound_lim;
+  CubedSphere::GetLocalIndex(&lv2_lx2, &lv2_lx3, &local_lx2, &local_lx3, &bound_lim, loc);
 
   // Work on interpolation
   AthenaArray<Real> interpolatedSrc;