diff --git a/lib/gpu/lal_neighbor.cpp b/lib/gpu/lal_neighbor.cpp
index 10816e2fa6..051f55f0a3 100644
--- a/lib/gpu/lal_neighbor.cpp
+++ b/lib/gpu/lal_neighbor.cpp
@@ -586,8 +586,25 @@ void Neighbor::build_nbor_list(double **x, const int inum, const int host_inum,
     const int b2y=_block_cell_2d;
     const int g2x=static_cast<int>(ceil(static_cast<double>(_maxspecial)/b2x));
     const int g2y=static_cast<int>(ceil(static_cast<double>(nt)/b2y));
-    _shared->k_transpose.set_size(g2x,g2y,b2x,b2y);
-    _shared->k_transpose.run(&dev_special,&dev_special_t,&_maxspecial,&nt);
+    // the maximum number of blocks on the device is typically 65535
+    // in principle we can use a lower number to have more resource per block 32768
+    const int max_num_blocks = 65535;
+    int shift = 0;
+    if (g2y < max_num_blocks) {
+      _shared->k_transpose.set_size(g2x,g2y,b2x,b2y);
+      _shared->k_transpose.run(&dev_special,&dev_special_t,&_maxspecial,&nt,&shift);
+    } else {
+      // using a fixed number of blocks
+      int g2y_m = max_num_blocks;
+      _shared->k_transpose.set_size(g2x,g2y_m,b2x,b2y);
+      // number of chunks needed for the whole transpose
+      const int num_chunks = ceil(static_cast<double>(g2y) / g2y_m);
+      for (int i = 0; i < num_chunks; i++) {        
+        _shared->k_transpose.run(&dev_special,&dev_special_t,&_maxspecial,&nt,&shift);
+        shift += g2y_m*b2y;
+      }
+    }
+
     time_transpose.stop();
   }
 
diff --git a/lib/gpu/lal_neighbor_gpu.cu b/lib/gpu/lal_neighbor_gpu.cu
index a7506fc5c3..7d0941ccd5 100644
--- a/lib/gpu/lal_neighbor_gpu.cu
+++ b/lib/gpu/lal_neighbor_gpu.cu
@@ -147,7 +147,7 @@ __kernel void kernel_calc_cell_counts(const unsigned *restrict cell_id,
 
 __kernel void transpose(__global tagint *restrict out,
                         const __global tagint *restrict in,
-                        int columns_in, int rows_in)
+                        int columns_in, int rows_in, int shift)
 {
   __local tagint block[BLOCK_CELL_2D][BLOCK_CELL_2D+1];
 
@@ -158,15 +158,15 @@ __kernel void transpose(__global tagint *restrict out,
 
   unsigned i=bi*BLOCK_CELL_2D+ti;
   unsigned j=bj*BLOCK_CELL_2D+tj;
-  if ((i<columns_in) && (j<rows_in))
-    block[tj][ti]=in[j*columns_in+i];
+  if ((i<columns_in) && (j+shift<rows_in))
+    block[tj][ti]=in[(j+shift)*columns_in+i];
 
    __syncthreads();
 
   i=bj*BLOCK_CELL_2D+ti;
   j=bi*BLOCK_CELL_2D+tj;
-  if ((i<rows_in) && (j<columns_in))
-    out[j*rows_in+i] = block[ti][tj];
+  if ((i+shift<rows_in) && (j<columns_in))
+    out[j*rows_in+i+shift] = block[ti][tj];
 }
 
 #ifndef LAL_USE_OLD_NEIGHBOR