   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 within computational domain, near the boundary
+    u[ti + (dom->Gfx._je-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];
+    u[ti + (dom->Gfx._jeb-1)*s1b + tk*s2b] = bc_n[tk + ti*dom->Gfx.knb];
     
   }
 }

  int s2b = dom->Gfy._s2b;

  if((ti < dom->Gfy._inb) && (tj < dom->Gfy._jnb)) {
    v[ti + tj*s1b + (dom->Gfy._keb-1)*s2b] = bc_s[ti + tj*dom->Gfy.inb];
    v[ti + tj*s1b + dom->Gfy._keb*s2b] = bc_n[ti + tj*dom->Gfy.inb];
  }
}


  int s2b = dom->Gfz._s2b;

  if((ti < dom->Gfz._inb) && (tk < dom->Gfz._knb)) {
    w[ti + (dom->Gfz._je-1)*s1b + tk*s2b] = bc_s[tk + ti*dom->Gfz.knb];
    //velocity on ghost cell
    w[ti + (dom->Gfz._jeb-1)*s1b + tk*s2b] = bc_n[tk + ti*dom->Gfz.knb];
  }
}

