   int s2b = dom->Gfx._s2b;
 
   if((ti < dom->Gfx._inb) && (tj < dom->Gfx._jnb)) {
-    u[ti + tj*s1b + (dom->Gfx._keb-1)*s2b] = bc_b[ti + tj*dom->Gfx.inb];
+    u[ti + tj*s1b + (dom->Gfx._ke-1)*s2b] = bc_b[ti + tj*dom->Gfx.inb];
     // velocity on ghost cell
-    u[ti + tj*s1b + dom->Gfx._keb*s2b] = bc_t[ti + tj*dom->Gfx.inb];
+    u[ti + tj*s1b + (dom->Gfx._keb-1)*s2b] = bc_t[ti + tj*dom->Gfx.inb];
   }
 }
 

  int s2b = dom->Gfx._s2b;

  if((ti < dom->Gfx._inb) && (tk < dom->Gfx._knb)) {
    // velocity with computational domain, near the boundary
    u[ti + (dom->Gfx._jeb-1)*s1b + tk*s2b] = bc_s[tk + ti*dom->Gfx.knb];
    //velocity on ghost cells
    u[ti + dom->Gfx._jeb*s1b + tk*s2b] = bc_n[tk + ti*dom->Gfx.knb];
    
  }
}

  int s2b = dom->Gfy._s2b;

  if((tj < dom->Gfy._jnb) && (tk < dom->Gfy._knb)) {
    v[(dom->Gfy._ie-1) + tj*s1b + tk*s2b] = bc_w[tj + tk*dom->Gfy.jnb];
    //velocity on ghost cell
    v[dom->Gfy._ieb-1 + tj*s1b + tk*s2b] = bc_e[tj + tk*dom->Gfy.jnb];
  }
}

