diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp
index b013e81c6..0868b2c1e 100644
--- a/src/common/m_variables_conversion.fpp
+++ b/src/common/m_variables_conversion.fpp
@@ -81,9 +81,6 @@ module m_variables_conversion
 
     end interface ! ============================================================
 
-    integer, public :: ixb, ixe, iyb, iye, izb, ize
-    !$acc declare create(ixb, ixe, iyb, iye, izb, ize)
-
     !! In simulation, gammas, pi_infs, and qvs are already declared in m_global_variables
 #ifndef MFC_SIMULATION
     real(kind(0d0)), allocatable, public, dimension(:) :: gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps
@@ -628,26 +625,7 @@ contains
 
         integer :: i, j
 
-#ifdef MFC_PRE_PROCESS
-        ixb = 0; iyb = 0; izb = 0; 
-        ixe = m; iye = n; ize = p; 
-#else
-        ixb = -buff_size
-        ixe = m - ixb
-
-        iyb = 0; iye = 0; izb = 0; ize = 0; 
-        if (n > 0) then
-            iyb = -buff_size; iye = n - iyb
-
-            if (p > 0) then
-                izb = -buff_size; ize = p - izb
-            end if
-        end if
-#endif
-
-!$acc enter data copyin(ixb, ixe, iyb, iye, izb, ize)
 !$acc enter data copyin(is1b, is1e, is2b, is2e, is3b, is3e)
-!$acc update device(ixb, ixe, iyb, iye, izb, ize)
 
 #ifdef MFC_SIMULATION
         @:ALLOCATE_GLOBAL(gammas (1:num_fluids))
@@ -784,15 +762,16 @@ contains
 
     !Initialize mv at the quadrature nodes based on the initialized moments and sigma
     subroutine s_initialize_mv(qK_cons_vf, mv)
+
         type(scalar_field), dimension(sys_size), intent(in) :: qK_cons_vf
-        real(kind(0d0)), dimension(ixb:, iyb:, izb:, 1:, 1:), intent(inout) :: mv
+        real(kind(0d0)), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: mv
 
         integer :: i, j, k, l
         real(kind(0d0)) :: mu, sig, nbub_sc
 
-        do l = izb, ize
-            do k = iyb, iye
-                do j = ixb, ixe
+        do l = idwbuff(3)%beg, idwbuff(3)%end
+            do k = idwbuff(2)%beg, idwbuff(2)%end
+                do j = idwbuff(1)%beg, idwbuff(1)%end
 
                     nbub_sc = qK_cons_vf(bubxb)%sf(j, k, l)
 
@@ -816,15 +795,15 @@ contains
     !Initialize pb at the quadrature nodes using isothermal relations (Preston model)
     subroutine s_initialize_pb(qK_cons_vf, mv, pb)
         type(scalar_field), dimension(sys_size), intent(in) :: qK_cons_vf
-        real(kind(0d0)), dimension(ixb:, iyb:, izb:, 1:, 1:), intent(in) :: mv
-        real(kind(0d0)), dimension(ixb:, iyb:, izb:, 1:, 1:), intent(inout) :: pb
+        real(kind(0d0)), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(in) :: mv
+        real(kind(0d0)), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb
 
         integer :: i, j, k, l
         real(kind(0d0)) :: mu, sig, nbub_sc
 
-        do l = izb, ize
-            do k = iyb, iye
-                do j = ixb, ixe
+        do l = idwbuff(3)%beg, idwbuff(3)%end
+            do k = idwbuff(2)%beg, idwbuff(2)%end
+                do j = idwbuff(1)%beg, idwbuff(1)%end
 
                     nbub_sc = qK_cons_vf(bubxb)%sf(j, k, l)
 
@@ -855,19 +834,16 @@ contains
         !! @param iz Index bounds in third coordinate direction
     subroutine s_convert_conservative_to_primitive_variables(qK_cons_vf, &
                                                              qK_prim_vf, &
-                                                             gm_alphaK_vf, &
-                                                             ix, iy, iz)
+                                                             ibounds, &
+                                                             gm_alphaK_vf)
 
         type(scalar_field), dimension(sys_size), intent(in) :: qK_cons_vf
         type(scalar_field), dimension(sys_size), intent(inout) :: qK_prim_vf
+        type(int_bounds_info), dimension(1:3), intent(in) :: ibounds
         type(scalar_field), &
             allocatable, optional, dimension(:), &
             intent(in) :: gm_alphaK_vf
 
-        type(int_bounds_info), optional, intent(in) :: ix, iy, iz
-
-        type(int_bounds_info) :: aix, aiy, aiz
-
         real(kind(0d0)), dimension(num_fluids) :: alpha_K, alpha_rho_K
         real(kind(0d0)), dimension(2) :: Re_K
         real(kind(0d0)) :: rho_K, gamma_K, pi_inf_K, qv_K, dyn_pres_K
@@ -910,14 +886,10 @@ contains
             end if
         #:endif
 
-        if (present(ix)) then; aix = ix; else; aix%beg = ixb; aix%end = ixe; end if
-        if (present(iy)) then; aiy = iy; else; aiy%beg = iyb; aiy%end = iye; end if
-        if (present(iz)) then; aiz = iz; else; aiz%beg = izb; aiz%end = ize; end if
-
-        !$acc parallel loop collapse(3) gang vector default(present) copyin(aix, aiy, aiz) private(alpha_K, alpha_rho_K, Re_K, nRtmp, rho_K, gamma_K, pi_inf_K, qv_K, dyn_pres_K, R3tmp, rhoyks)
-        do l = aiz%beg, aiz%end
-            do k = aiy%beg, aiy%end
-                do j = aix%beg, aix%end
+        !$acc parallel loop collapse(3) gang vector default(present) private(alpha_K, alpha_rho_K, Re_K, nRtmp, rho_K, gamma_K, pi_inf_K, qv_K, dyn_pres_K, R3tmp, rhoyks)
+        do l = ibounds(3)%beg, ibounds(3)%end
+            do k = ibounds(2)%beg, ibounds(2)%end
+                do j = ibounds(1)%beg, ibounds(1)%end
                     dyn_pres_K = 0d0
 
                     !$acc loop seq
diff --git a/src/post_process/m_global_parameters.fpp b/src/post_process/m_global_parameters.fpp
index cfd9036c2..64623f958 100644
--- a/src/post_process/m_global_parameters.fpp
+++ b/src/post_process/m_global_parameters.fpp
@@ -130,6 +130,13 @@ module m_global_parameters
     type(int_bounds_info) :: temperature_idx       !< Indexes of first & last temperature eqns.
     !> @}
 
+    ! Cell Indices for the (local) interior points (O-m, O-n, 0-p).
+    type(int_bounds_info) :: idwint(1:3)
+
+    ! Cell Indices for the entire (local) domain. In simulation, this includes
+    ! the buffer region. idwbuff and idwint are the same otherwise.
+    type(int_bounds_info) :: idwbuff(1:3)
+
     !> @name Boundary conditions in the x-, y- and z-coordinate directions
     !> @{
     type(int_bounds_info) :: bc_x, bc_y, bc_z
@@ -259,7 +266,6 @@ module m_global_parameters
     real(kind(0d0)) :: poly_sigma
     real(kind(0d0)) :: sigR
     integer :: nmom
-
     !> @}
 
     !> @name surface tension coefficient
@@ -642,8 +648,6 @@ contains
         tempxb = temperature_idx%beg
         tempxe = temperature_idx%end
 
-        ! ==================================================================
-
 #ifdef MFC_MPI
         allocate (MPI_IO_DATA%view(1:sys_size))
         allocate (MPI_IO_DATA%var(1:sys_size))
@@ -691,6 +695,19 @@ contains
             buff_size = buff_size + fd_number
         end if
 
+        ! Configuring Coordinate Direction Indexes =========================
+        idwint(1)%beg = 0; idwint(2)%beg = 0; idwint(3)%beg = 0
+        idwint(1)%end = m; idwint(2)%end = n; idwint(3)%end = p
+
+        idwbuff(1)%beg = -buff_size
+        if (num_dims > 1) then; idwbuff(2)%beg = -buff_size; else; idwbuff(2)%beg = 0; end if
+        if (num_dims > 2) then; idwbuff(3)%beg = -buff_size; else; idwbuff(3)%beg = 0; end if
+
+        idwbuff(1)%end = idwint(1)%end - idwbuff(1)%beg
+        idwbuff(2)%end = idwint(2)%end - idwbuff(2)%beg
+        idwbuff(3)%end = idwint(3)%end - idwbuff(3)%beg
+        ! ==================================================================
+
         ! Allocating single precision grid variables if needed
         if (precision == 1) then
             allocate (x_cb_s(-1 - offset_x%beg:m + offset_x%end))
diff --git a/src/post_process/m_start_up.f90 b/src/post_process/m_start_up.f90
index 062fb614c..f1d48a831 100644
--- a/src/post_process/m_start_up.f90
+++ b/src/post_process/m_start_up.f90
@@ -180,7 +180,7 @@ subroutine s_perform_time_step(t_step)
         end if
 
         ! Converting the conservative variables to the primitive ones
-        call s_convert_conservative_to_primitive_variables(q_cons_vf, q_prim_vf)
+        call s_convert_conservative_to_primitive_variables(q_cons_vf, q_prim_vf, idwbuff)
 
     end subroutine s_perform_time_step
 
diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp
index d34cad914..1c032bba7 100644
--- a/src/pre_process/m_global_parameters.fpp
+++ b/src/pre_process/m_global_parameters.fpp
@@ -108,6 +108,13 @@ module m_global_parameters
     type(int_bounds_info) :: species_idx           !< Indexes of first & last concentration eqns.
     type(int_bounds_info) :: temperature_idx       !< Indexes of first & last temperature eqns.
 
+    ! Cell Indices for the (local) interior points (O-m, O-n, 0-p).
+    type(int_bounds_info) :: idwint(1:3)
+
+    ! Cell Indices for the entire (local) domain. In simulation and post_process,
+    ! this includes the buffer region. idwbuff and idwint are the same otherwise.
+    type(int_bounds_info) :: idwbuff(1:3)
+
     type(int_bounds_info) :: bc_x, bc_y, bc_z !<
     !! Boundary conditions in the x-, y- and z-coordinate directions
 
@@ -722,6 +729,12 @@ contains
         tempxb = temperature_idx%beg
         tempxe = temperature_idx%end
 
+        ! Configuring Coordinate Direction Indexes =========================
+        idwint(1)%beg = 0; idwint(2)%beg = 0; idwint(3)%beg = 0
+        idwint(1)%end = m; idwint(2)%end = n; idwint(3)%end = p
+
+        ! There is no buffer region in pre_process.
+        idwbuff(1) = idwint(1); idwbuff(2) = idwint(2); idwbuff(3) = idwint(3)
         ! ==================================================================
 
 #ifdef MFC_MPI
diff --git a/src/pre_process/m_initial_condition.fpp b/src/pre_process/m_initial_condition.fpp
index 2e5ecde00..7cb56e114 100644
--- a/src/pre_process/m_initial_condition.fpp
+++ b/src/pre_process/m_initial_condition.fpp
@@ -124,7 +124,8 @@ contains
         ! preexisting initial condition data files were read in on start-up
         if (old_ic) then
             call s_convert_conservative_to_primitive_variables(q_cons_vf, &
-                                                               q_prim_vf)
+                                                               q_prim_vf, &
+                                                               idwbuff)
         end if
 
         !  3D Patch Geometries =============================================
diff --git a/src/simulation/m_bubbles.fpp b/src/simulation/m_bubbles.fpp
index 21f7aaf58..804b473b8 100644
--- a/src/simulation/m_bubbles.fpp
+++ b/src/simulation/m_bubbles.fpp
@@ -58,15 +58,6 @@ contains
     subroutine s_initialize_bubbles_module
 
         integer :: i, j, k, l, q
-        type(int_bounds_info) :: ix, iy, iz
-
-        ! Configuring Coordinate Direction Indexes =========================
-        ix%beg = -buff_size; iy%beg = 0; iz%beg = 0
-
-        if (n > 0) iy%beg = -buff_size; if (p > 0) iz%beg = -buff_size
-
-        ix%end = m - ix%beg; iy%end = n - iy%beg; iz%end = p - iz%beg
-        ! ==================================================================
 
         @:ALLOCATE_GLOBAL(rs(1:nb))
         @:ALLOCATE_GLOBAL(vs(1:nb))
@@ -89,7 +80,7 @@ contains
             !$acc update device(ps, ms)
         end if
 
-        @:ALLOCATE(divu%sf(ix%beg:ix%end, iy%beg:iy%end, iz%beg:iz%end))
+        @:ALLOCATE(divu%sf(idwbuff(1)%beg:idwbuff(1)%end, idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end))
         @:ACC_SETUP_SFs(divu)
 
         @:ALLOCATE_GLOBAL(bub_adv_src(0:m, 0:n, 0:p))
diff --git a/src/simulation/m_chemistry.fpp b/src/simulation/m_chemistry.fpp
index 3a773fc1f..ae3aa3f83 100644
--- a/src/simulation/m_chemistry.fpp
+++ b/src/simulation/m_chemistry.fpp
@@ -17,10 +17,8 @@ module m_chemistry
 
     implicit none
 
-    type(int_bounds_info), private :: ix, iy, iz
     type(scalar_field), private :: grads(1:3)
 
-    !$acc declare create(ix, iy, iz)
     !$acc declare create(grads)
 
 contains
@@ -29,16 +27,11 @@ contains
 
         integer :: i
 
-        ix%beg = -buff_size
-        if (n > 0) then; iy%beg = -buff_size; else; iy%beg = 0; end if
-        if (p > 0) then; iz%beg = -buff_size; else; iz%beg = 0; end if
-
-        ix%end = m - ix%beg; iy%end = n - iy%beg; iz%end = p - iz%beg
-
-        !$acc update device(ix, iy, iz)
-
         do i = 1, 3
-            @:ALLOCATE(grads(i)%sf(ix%beg:ix%end, iy%beg:iy%end, iz%beg:iz%end))
+            @:ALLOCATE(grads(i)%sf(&
+                & idwbuff(1)%beg:idwbuff(1)%end, &
+                & idwbuff(2)%beg:idwbuff(2)%end, &
+                & idwbuff(3)%beg:idwbuff(3)%end))
         end do
 
         !$acc kernels
diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp
index 22698cfb5..e2c1ccf4b 100644
--- a/src/simulation/m_data_output.fpp
+++ b/src/simulation/m_data_output.fpp
@@ -510,7 +510,7 @@ contains
         if (.not. file_exist) call s_create_directory(trim(t_step_dir))
 
         if (prim_vars_wrt .or. (n == 0 .and. p == 0)) then
-            call s_convert_conservative_to_primitive_variables(q_cons_vf, q_prim_vf)
+            call s_convert_conservative_to_primitive_variables(q_cons_vf, q_prim_vf, idwbuff)
             do i = 1, sys_size
                 !$acc update host(q_prim_vf(i)%sf(:,:,:))
             end do
diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp
index 5b071117b..b8cf82cdf 100644
--- a/src/simulation/m_global_parameters.fpp
+++ b/src/simulation/m_global_parameters.fpp
@@ -238,6 +238,15 @@ module m_global_parameters
 
     !$acc declare create(bub_idx)
 
+    ! Cell Indices for the (local) interior points (O-m, O-n, 0-p).
+    type(int_bounds_info) :: idwint(1:3)
+    !$acc declare create(idwint)
+
+    ! Cell Indices for the entire (local) domain. In simulation and post_process,
+    ! this includes the buffer region. idwbuff and idwint are the same otherwise.
+    type(int_bounds_info) :: idwbuff(1:3)
+    !$acc declare create(idwbuff)
+
     !> @name The number of fluids, along with their identifying indexes, respectively,
     !! for which viscous effects, e.g. the shear and/or the volume Reynolds (Re)
     !! numbers, will be non-negligible.
@@ -702,8 +711,6 @@ contains
         integer :: i, j, k
         integer :: fac
 
-        type(int_bounds_info) :: ix, iy, iz
-
         #:if not MFC_CASE_OPTIMIZATION
             ! Determining the degree of the WENO polynomials
             weno_polyn = (weno_order - 1)/2
@@ -1037,18 +1044,25 @@ contains
         end if
 
         ! Configuring Coordinate Direction Indexes =========================
-        if (bubbles) then
-            ix%beg = -buff_size; iy%beg = 0; iz%beg = 0
-            if (n > 0) then
-                iy%beg = -buff_size
-                if (p > 0) then
-                    iz%beg = -buff_size
-                end if
-            end if
+        idwint(1)%beg = 0; idwint(2)%beg = 0; idwint(3)%beg = 0
+        idwint(1)%end = m; idwint(2)%end = n; idwint(3)%end = p
+
+        idwbuff(1)%beg = -buff_size
+        if (num_dims > 1) then; idwbuff(2)%beg = -buff_size; else; idwbuff(2)%beg = 0; end if
+        if (num_dims > 2) then; idwbuff(3)%beg = -buff_size; else; idwbuff(3)%beg = 0; end if
 
-            ix%end = m - ix%beg; iy%end = n - iy%beg; iz%end = p - iz%beg
+        idwbuff(1)%end = idwint(1)%end - idwbuff(1)%beg
+        idwbuff(2)%end = idwint(2)%end - idwbuff(2)%beg
+        idwbuff(3)%end = idwint(3)%end - idwbuff(3)%beg
+        !$acc update device(idwint, idwbuff)
+        ! ==================================================================
 
-            @:ALLOCATE_GLOBAL(ptil(ix%beg:ix%end, iy%beg:iy%end, iz%beg:iz%end))
+        ! Configuring Coordinate Direction Indexes =========================
+        if (bubbles) then
+            @:ALLOCATE_GLOBAL(ptil(&
+                & idwbuff(1)%beg:idwbuff(1)%end, &
+                & idwbuff(2)%beg:idwbuff(2)%end, &
+                & idwbuff(3)%beg:idwbuff(3)%end))
         end if
 
         if (probe_wrt) then
diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp
index 7ece8f5c5..360c484d0 100644
--- a/src/simulation/m_rhs.fpp
+++ b/src/simulation/m_rhs.fpp
@@ -169,21 +169,12 @@ module m_rhs
 
     !> @name Indical bounds in the x-, y- and z-directions
     !> @{
-    type(int_bounds_info) :: ix, iy, iz
-    !$acc declare create(ix, iy, iz)
-
-    type(int_bounds_info) :: ibounds_interior(1:3)
-    !$acc declare create(ibounds_interior)
-
-    type(int_bounds_info) :: ibounds_wbuffer(1:3)
-    !$acc declare create(ibounds_wbuffer)
+    type(int_bounds_info) :: irx, iry, irz
+    !$acc declare create(irx, iry, irz)
 
     type(int_bounds_info) :: is1, is2, is3
     !$acc declare create(is1, is2, is3)
 
-    type(int_bounds_info) :: ixt, iyt, izt
-    !$acc declare create(ixt, iyt, izt)
-
     !> @name Saved fluxes for testing
     !> @{
     type(scalar_field) :: alf_sum
@@ -241,37 +232,18 @@ contains
 
         integer :: i, j, k, l, id !< Generic loop iterators
 
-        ! Configuring Coordinate Direction Indexes =========================
-        ibounds_interior(1)%beg = 0; ibounds_interior(2)%beg = 0; ibounds_interior(3)%beg = 0
-        ibounds_interior(1)%end = m; ibounds_interior(2)%end = n; ibounds_interior(3)%end = p
-
-        ibounds_wbuffer(1)%beg = -buff_size
-        ibounds_wbuffer(2)%beg = 0
-        ibounds_wbuffer(3)%beg = 0
-        if (num_dims > 1) then ibounds_wbuffer(2)%beg = -buff_size; end if
-        if (num_dims > 2) then ibounds_wbuffer(3)%beg = -buff_size; end if
-
-        ibounds_wbuffer(1)%end = m + buff_size
-        ibounds_wbuffer(2)%end = n + buff_size
-        ibounds_wbuffer(3)%end = p + buff_size
-
-        ix = ibounds_interior(1); iy = ibounds_interior(2); iz = ibounds_interior(3)
-        ! ==================================================================
-
-        !$acc enter data copyin(ibounds_interior, ibounds_wbuffer, ix, iy, iz)
-        !$acc update device(ibounds_interior, ibounds_wbuffer, ix, iy, iz)
-
-        ixt = ix; iyt = iy; izt = iz
+        !$acc enter data copyin(idwbuff, idwbuff)
+        !$acc update device(idwbuff, idwbuff)
 
         @:ALLOCATE(q_cons_qp%vf(1:sys_size))
         @:ALLOCATE(q_prim_qp%vf(1:sys_size))
 
         do l = 1, sys_size
-            @:ALLOCATE(q_cons_qp%vf(l)%sf(ix%beg:ix%end, iy%beg:iy%end, iz%beg:iz%end))
+            @:ALLOCATE(q_cons_qp%vf(l)%sf(idwbuff(1)%beg:idwbuff(1)%end, idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end))
         end do
 
         do l = mom_idx%beg, E_idx
-            @:ALLOCATE(q_prim_qp%vf(l)%sf(ix%beg:ix%end, iy%beg:iy%end, iz%beg:iz%end))
+            @:ALLOCATE(q_prim_qp%vf(l)%sf(idwbuff(1)%beg:idwbuff(1)%end, idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end))
         end do
 
         if (.not. f_is_default(sigma)) then
@@ -279,11 +251,11 @@ contains
             ! the last equation. If this changes then this logic will
             ! need updated
             do l = adv_idx%end + 1, sys_size - 1
-                @:ALLOCATE(q_prim_qp%vf(l)%sf(ix%beg:ix%end, iy%beg:iy%end, iz%beg:iz%end))
+                @:ALLOCATE(q_prim_qp%vf(l)%sf(idwbuff(1)%beg:idwbuff(1)%end, idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end))
             end do
         else
             do l = adv_idx%end + 1, sys_size
-                @:ALLOCATE(q_prim_qp%vf(l)%sf(ix%beg:ix%end, iy%beg:iy%end, iz%beg:iz%end))
+                @:ALLOCATE(q_prim_qp%vf(l)%sf(idwbuff(1)%beg:idwbuff(1)%end, idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end))
             end do
 
         end if
@@ -312,14 +284,14 @@ contains
         if (any(Re_size > 0)) then
             @:ALLOCATE_GLOBAL(tau_Re_vf(1:sys_size))
             do i = 1, num_dims
-                @:ALLOCATE(tau_Re_vf(cont_idx%end + i)%sf(ix%beg:ix%end, &
-                                                    &  iy%beg:iy%end, &
-                                                    &  iz%beg:iz%end))
+                @:ALLOCATE(tau_Re_vf(cont_idx%end + i)%sf(idwbuff(1)%beg:idwbuff(1)%end, &
+                                                    &  idwbuff(2)%beg:idwbuff(2)%end, &
+                                                    &  idwbuff(3)%beg:idwbuff(3)%end))
                 @:ACC_SETUP_SFs(tau_Re_vf(cont_idx%end + i))
             end do
-            @:ALLOCATE(tau_Re_vf(E_idx)%sf(ix%beg:ix%end, &
-                                        & iy%beg:iy%end, &
-                                        & iz%beg:iz%end))
+            @:ALLOCATE(tau_Re_vf(E_idx)%sf(idwbuff(1)%beg:idwbuff(1)%end, &
+                                        & idwbuff(2)%beg:idwbuff(2)%end, &
+                                        & idwbuff(3)%beg:idwbuff(3)%end))
             @:ACC_SETUP_SFs(tau_Re_vf(E_idx))
         end if
 
@@ -332,9 +304,9 @@ contains
                 do j = 0, 2
                     do k = 1, nb
                         @:ALLOCATE(mom_3d(i, j, k)%sf( &
-                                      & ix%beg:ix%end, &
-                                      & iy%beg:iy%end, &
-                                      & iz%beg:iz%end))
+                                      & idwbuff(1)%beg:idwbuff(1)%end, &
+                                      & idwbuff(2)%beg:idwbuff(2)%end, &
+                                      & idwbuff(3)%beg:idwbuff(3)%end))
                         @:ACC_SETUP_SFs(mom_3d(i, j, k))
                     end do
                 end do
@@ -342,9 +314,9 @@ contains
 
             do i = 1, nmomsp
                 @:ALLOCATE(mom_sp(i)%sf( &
-                        & ix%beg:ix%end, &
-                        & iy%beg:iy%end, &
-                        & iz%beg:iz%end))
+                        & idwbuff(1)%beg:idwbuff(1)%end, &
+                        & idwbuff(2)%beg:idwbuff(2)%end, &
+                        & idwbuff(3)%beg:idwbuff(3)%end))
                 @:ACC_SETUP_SFs(mom_sp(i))
             end do
         end if
@@ -357,45 +329,45 @@ contains
             @:ALLOCATE(qL_prim(i)%vf(1:sys_size))
             @:ALLOCATE(qR_prim(i)%vf(1:sys_size))
             do l = mom_idx%beg, mom_idx%end
-                @:ALLOCATE(qL_prim(i)%vf(l)%sf(ix%beg:ix%end, iy%beg:iy%end, iz%beg:iz%end))
-                @:ALLOCATE(qR_prim(i)%vf(l)%sf(ix%beg:ix%end, iy%beg:iy%end, iz%beg:iz%end))
+                @:ALLOCATE(qL_prim(i)%vf(l)%sf(idwbuff(1)%beg:idwbuff(1)%end, idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end))
+                @:ALLOCATE(qR_prim(i)%vf(l)%sf(idwbuff(1)%beg:idwbuff(1)%end, idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end))
             end do
             @:ACC_SETUP_VFs(qL_prim(i), qR_prim(i))
         end do
 
         if (mpp_lim .and. bubbles) then
-            @:ALLOCATE(alf_sum%sf(ix%beg:ix%end, iy%beg:iy%end, iz%beg:iz%end))
+            @:ALLOCATE(alf_sum%sf(idwbuff(1)%beg:idwbuff(1)%end, idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end))
         end if
         ! END: Allocation/Association of qK_cons_n and qK_prim_n ======
 
-        @:ALLOCATE_GLOBAL(qL_rsx_vf(ix%beg:ix%end, &
-            iy%beg:iy%end, iz%beg:iz%end, 1:sys_size))
-        @:ALLOCATE_GLOBAL(qR_rsx_vf(ix%beg:ix%end, &
-            iy%beg:iy%end, iz%beg:iz%end, 1:sys_size))
+        @:ALLOCATE_GLOBAL(qL_rsx_vf(idwbuff(1)%beg:idwbuff(1)%end, &
+            idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end, 1:sys_size))
+        @:ALLOCATE_GLOBAL(qR_rsx_vf(idwbuff(1)%beg:idwbuff(1)%end, &
+            idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end, 1:sys_size))
 
         if (n > 0) then
 
-            @:ALLOCATE_GLOBAL(qL_rsy_vf(iy%beg:iy%end, &
-                ix%beg:ix%end, iz%beg:iz%end, 1:sys_size))
-            @:ALLOCATE_GLOBAL(qR_rsy_vf(iy%beg:iy%end, &
-                ix%beg:ix%end, iz%beg:iz%end, 1:sys_size))
+            @:ALLOCATE_GLOBAL(qL_rsy_vf(idwbuff(2)%beg:idwbuff(2)%end, &
+                idwbuff(1)%beg:idwbuff(1)%end, idwbuff(3)%beg:idwbuff(3)%end, 1:sys_size))
+            @:ALLOCATE_GLOBAL(qR_rsy_vf(idwbuff(2)%beg:idwbuff(2)%end, &
+                idwbuff(1)%beg:idwbuff(1)%end, idwbuff(3)%beg:idwbuff(3)%end, 1:sys_size))
         else
-            @:ALLOCATE_GLOBAL(qL_rsy_vf(ix%beg:ix%end, &
-                iy%beg:iy%end, iz%beg:iz%end, 1:sys_size))
-            @:ALLOCATE_GLOBAL(qR_rsy_vf(ix%beg:ix%end, &
-                iy%beg:iy%end, iz%beg:iz%end, 1:sys_size))
+            @:ALLOCATE_GLOBAL(qL_rsy_vf(idwbuff(1)%beg:idwbuff(1)%end, &
+                idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end, 1:sys_size))
+            @:ALLOCATE_GLOBAL(qR_rsy_vf(idwbuff(1)%beg:idwbuff(1)%end, &
+                idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end, 1:sys_size))
         end if
 
         if (p > 0) then
-            @:ALLOCATE_GLOBAL(qL_rsz_vf(iz%beg:iz%end, &
-                iy%beg:iy%end, ix%beg:ix%end, 1:sys_size))
-            @:ALLOCATE_GLOBAL(qR_rsz_vf(iz%beg:iz%end, &
-                iy%beg:iy%end, ix%beg:ix%end, 1:sys_size))
+            @:ALLOCATE_GLOBAL(qL_rsz_vf(idwbuff(3)%beg:idwbuff(3)%end, &
+                idwbuff(2)%beg:idwbuff(2)%end, idwbuff(1)%beg:idwbuff(1)%end, 1:sys_size))
+            @:ALLOCATE_GLOBAL(qR_rsz_vf(idwbuff(3)%beg:idwbuff(3)%end, &
+                idwbuff(2)%beg:idwbuff(2)%end, idwbuff(1)%beg:idwbuff(1)%end, 1:sys_size))
         else
-            @:ALLOCATE_GLOBAL(qL_rsz_vf(ix%beg:ix%end, &
-                iy%beg:iy%end, iz%beg:iz%end, 1:sys_size))
-            @:ALLOCATE_GLOBAL(qR_rsz_vf(ix%beg:ix%end, &
-                iy%beg:iy%end, iz%beg:iz%end, 1:sys_size))
+            @:ALLOCATE_GLOBAL(qL_rsz_vf(idwbuff(1)%beg:idwbuff(1)%end, &
+                idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end, 1:sys_size))
+            @:ALLOCATE_GLOBAL(qR_rsz_vf(idwbuff(1)%beg:idwbuff(1)%end, &
+                idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end, 1:sys_size))
 
         end if
 
@@ -413,9 +385,9 @@ contains
 
                 do l = mom_idx%beg, mom_idx%end
                     @:ALLOCATE(dq_prim_dx_qp(1)%vf(l)%sf( &
-                              & ix%beg:ix%end, &
-                              & iy%beg:iy%end, &
-                              & iz%beg:iz%end))
+                              & idwbuff(1)%beg:idwbuff(1)%end, &
+                              & idwbuff(2)%beg:idwbuff(2)%end, &
+                              & idwbuff(3)%beg:idwbuff(3)%end))
                 end do
 
                 @:ACC_SETUP_VFs(dq_prim_dx_qp(1))
@@ -424,9 +396,9 @@ contains
 
                     do l = mom_idx%beg, mom_idx%end
                         @:ALLOCATE(dq_prim_dy_qp(1)%vf(l)%sf( &
-                                 & ix%beg:ix%end, &
-                                 & iy%beg:iy%end, &
-                                 & iz%beg:iz%end))
+                                 & idwbuff(1)%beg:idwbuff(1)%end, &
+                                 & idwbuff(2)%beg:idwbuff(2)%end, &
+                                 & idwbuff(3)%beg:idwbuff(3)%end))
                     end do
 
                     @:ACC_SETUP_VFs(dq_prim_dy_qp(1))
@@ -435,9 +407,9 @@ contains
 
                         do l = mom_idx%beg, mom_idx%end
                             @:ALLOCATE(dq_prim_dz_qp(1)%vf(l)%sf( &
-                                     & ix%beg:ix%end, &
-                                     & iy%beg:iy%end, &
-                                     & iz%beg:iz%end))
+                                     & idwbuff(1)%beg:idwbuff(1)%end, &
+                                     & idwbuff(2)%beg:idwbuff(2)%end, &
+                                     & idwbuff(3)%beg:idwbuff(3)%end))
                         end do
                         @:ACC_SETUP_VFs(dq_prim_dz_qp(1))
                     end if
@@ -487,38 +459,38 @@ contains
 
                     do l = mom_idx%beg, mom_idx%end
                         @:ALLOCATE(dqL_prim_dx_n(i)%vf(l)%sf( &
-                                 & ix%beg:ix%end, &
-                                 & iy%beg:iy%end, &
-                                 & iz%beg:iz%end))
+                                 & idwbuff(1)%beg:idwbuff(1)%end, &
+                                 & idwbuff(2)%beg:idwbuff(2)%end, &
+                                 & idwbuff(3)%beg:idwbuff(3)%end))
                         @:ALLOCATE(dqR_prim_dx_n(i)%vf(l)%sf( &
-                                 & ix%beg:ix%end, &
-                                 & iy%beg:iy%end, &
-                                 & iz%beg:iz%end))
+                                 & idwbuff(1)%beg:idwbuff(1)%end, &
+                                 & idwbuff(2)%beg:idwbuff(2)%end, &
+                                 & idwbuff(3)%beg:idwbuff(3)%end))
                     end do
 
                     if (n > 0) then
                         do l = mom_idx%beg, mom_idx%end
                             @:ALLOCATE(dqL_prim_dy_n(i)%vf(l)%sf( &
-                                     & ix%beg:ix%end, &
-                                     & iy%beg:iy%end, &
-                                     & iz%beg:iz%end))
+                                     & idwbuff(1)%beg:idwbuff(1)%end, &
+                                     & idwbuff(2)%beg:idwbuff(2)%end, &
+                                     & idwbuff(3)%beg:idwbuff(3)%end))
                             @:ALLOCATE(dqR_prim_dy_n(i)%vf(l)%sf( &
-                                     & ix%beg:ix%end, &
-                                     & iy%beg:iy%end, &
-                                     & iz%beg:iz%end))
+                                     & idwbuff(1)%beg:idwbuff(1)%end, &
+                                     & idwbuff(2)%beg:idwbuff(2)%end, &
+                                     & idwbuff(3)%beg:idwbuff(3)%end))
                         end do
                     end if
 
                     if (p > 0) then
                         do l = mom_idx%beg, mom_idx%end
                             @:ALLOCATE(dqL_prim_dz_n(i)%vf(l)%sf( &
-                                     & ix%beg:ix%end, &
-                                     & iy%beg:iy%end, &
-                                     & iz%beg:iz%end))
+                                     & idwbuff(1)%beg:idwbuff(1)%end, &
+                                     & idwbuff(2)%beg:idwbuff(2)%end, &
+                                     & idwbuff(3)%beg:idwbuff(3)%end))
                             @:ALLOCATE(dqR_prim_dz_n(i)%vf(l)%sf( &
-                                     & ix%beg:ix%end, &
-                                     & iy%beg:iy%end, &
-                                     & iz%beg:iz%end))
+                                     & idwbuff(1)%beg:idwbuff(1)%end, &
+                                     & idwbuff(2)%beg:idwbuff(2)%end, &
+                                     & idwbuff(3)%beg:idwbuff(3)%end))
                         end do
                     end if
 
@@ -532,35 +504,35 @@ contains
 
         if (any(Re_size > 0)) then
             if (weno_Re_flux) then
-                @:ALLOCATE_GLOBAL(dqL_rsx_vf(ix%beg:ix%end, &
-                    iy%beg:iy%end, iz%beg:iz%end, mom_idx%beg:mom_idx%end))
-                @:ALLOCATE_GLOBAL(dqR_rsx_vf(ix%beg:ix%end, &
-                    iy%beg:iy%end, iz%beg:iz%end, mom_idx%beg:mom_idx%end))
+                @:ALLOCATE_GLOBAL(dqL_rsx_vf(idwbuff(1)%beg:idwbuff(1)%end, &
+                    idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end, mom_idx%beg:mom_idx%end))
+                @:ALLOCATE_GLOBAL(dqR_rsx_vf(idwbuff(1)%beg:idwbuff(1)%end, &
+                    idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end, mom_idx%beg:mom_idx%end))
 
                 if (n > 0) then
 
-                    @:ALLOCATE_GLOBAL(dqL_rsy_vf(iy%beg:iy%end, &
-                        ix%beg:ix%end, iz%beg:iz%end, mom_idx%beg:mom_idx%end))
-                    @:ALLOCATE_GLOBAL(dqR_rsy_vf(iy%beg:iy%end, &
-                        ix%beg:ix%end, iz%beg:iz%end, mom_idx%beg:mom_idx%end))
+                    @:ALLOCATE_GLOBAL(dqL_rsy_vf(idwbuff(2)%beg:idwbuff(2)%end, &
+                        idwbuff(1)%beg:idwbuff(1)%end, idwbuff(3)%beg:idwbuff(3)%end, mom_idx%beg:mom_idx%end))
+                    @:ALLOCATE_GLOBAL(dqR_rsy_vf(idwbuff(2)%beg:idwbuff(2)%end, &
+                        idwbuff(1)%beg:idwbuff(1)%end, idwbuff(3)%beg:idwbuff(3)%end, mom_idx%beg:mom_idx%end))
                 else
-                    @:ALLOCATE_GLOBAL(dqL_rsy_vf(ix%beg:ix%end, &
-                        iy%beg:iy%end, iz%beg:iz%end, mom_idx%beg:mom_idx%end))
-                    @:ALLOCATE_GLOBAL(dqR_rsy_vf(ix%beg:ix%end, &
-                        iy%beg:iy%end, iz%beg:iz%end, mom_idx%beg:mom_idx%end))
+                    @:ALLOCATE_GLOBAL(dqL_rsy_vf(idwbuff(1)%beg:idwbuff(1)%end, &
+                        idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end, mom_idx%beg:mom_idx%end))
+                    @:ALLOCATE_GLOBAL(dqR_rsy_vf(idwbuff(1)%beg:idwbuff(1)%end, &
+                        idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end, mom_idx%beg:mom_idx%end))
 
                 end if
 
                 if (p > 0) then
-                    @:ALLOCATE_GLOBAL(dqL_rsz_vf(iz%beg:iz%end, &
-                        iy%beg:iy%end, ix%beg:ix%end, mom_idx%beg:mom_idx%end))
-                    @:ALLOCATE_GLOBAL(dqR_rsz_vf(iz%beg:iz%end, &
-                        iy%beg:iy%end, ix%beg:ix%end, mom_idx%beg:mom_idx%end))
+                    @:ALLOCATE_GLOBAL(dqL_rsz_vf(idwbuff(3)%beg:idwbuff(3)%end, &
+                        idwbuff(2)%beg:idwbuff(2)%end, idwbuff(1)%beg:idwbuff(1)%end, mom_idx%beg:mom_idx%end))
+                    @:ALLOCATE_GLOBAL(dqR_rsz_vf(idwbuff(3)%beg:idwbuff(3)%end, &
+                        idwbuff(2)%beg:idwbuff(2)%end, idwbuff(1)%beg:idwbuff(1)%end, mom_idx%beg:mom_idx%end))
                 else
-                    @:ALLOCATE_GLOBAL(dqL_rsz_vf(ix%beg:ix%end, &
-                        iy%beg:iy%end, iz%beg:iz%end, mom_idx%beg:mom_idx%end))
-                    @:ALLOCATE_GLOBAL(dqR_rsz_vf(ix%beg:ix%end, &
-                        iy%beg:iy%end, iz%beg:iz%end, mom_idx%beg:mom_idx%end))
+                    @:ALLOCATE_GLOBAL(dqL_rsz_vf(idwbuff(1)%beg:idwbuff(1)%end, &
+                        idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end, mom_idx%beg:mom_idx%end))
+                    @:ALLOCATE_GLOBAL(dqR_rsz_vf(idwbuff(1)%beg:idwbuff(1)%end, &
+                        idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end, mom_idx%beg:mom_idx%end))
 
                 end if
             end if
@@ -587,53 +559,53 @@ contains
             if (i == 1) then
                 do l = 1, sys_size
                     @:ALLOCATE(flux_n(i)%vf(l)%sf( &
-                             & ix%beg:ix%end, &
-                             & iy%beg:iy%end, &
-                             & iz%beg:iz%end))
+                             & idwbuff(1)%beg:idwbuff(1)%end, &
+                             & idwbuff(2)%beg:idwbuff(2)%end, &
+                             & idwbuff(3)%beg:idwbuff(3)%end))
                     @:ALLOCATE(flux_gsrc_n(i)%vf(l)%sf( &
-                            & ix%beg:ix%end, &
-                            & iy%beg:iy%end, &
-                            & iz%beg:iz%end))
+                            & idwbuff(1)%beg:idwbuff(1)%end, &
+                            & idwbuff(2)%beg:idwbuff(2)%end, &
+                            & idwbuff(3)%beg:idwbuff(3)%end))
                 end do
 
                 if (any(Re_size > 0) .or. (.not. f_is_default(sigma))) then
                     do l = mom_idx%beg, E_idx
                         @:ALLOCATE(flux_src_n(i)%vf(l)%sf( &
-                                 & ix%beg:ix%end, &
-                                 & iy%beg:iy%end, &
-                                 & iz%beg:iz%end))
+                                 & idwbuff(1)%beg:idwbuff(1)%end, &
+                                 & idwbuff(2)%beg:idwbuff(2)%end, &
+                                 & idwbuff(3)%beg:idwbuff(3)%end))
                     end do
                 end if
 
                 @:ALLOCATE(flux_src_n(i)%vf(adv_idx%beg)%sf( &
-                         & ix%beg:ix%end, &
-                         & iy%beg:iy%end, &
-                         & iz%beg:iz%end))
+                         & idwbuff(1)%beg:idwbuff(1)%end, &
+                         & idwbuff(2)%beg:idwbuff(2)%end, &
+                         & idwbuff(3)%beg:idwbuff(3)%end))
 
                 if (riemann_solver == 1) then
                     do l = adv_idx%beg + 1, adv_idx%end
                         @:ALLOCATE(flux_src_n(i)%vf(l)%sf( &
-                                 & ix%beg:ix%end, &
-                                 & iy%beg:iy%end, &
-                                 & iz%beg:iz%end))
+                                 & idwbuff(1)%beg:idwbuff(1)%end, &
+                                 & idwbuff(2)%beg:idwbuff(2)%end, &
+                                 & idwbuff(3)%beg:idwbuff(3)%end))
                     end do
                 end if
 
                 if (chemistry) then
                     do l = chemxb, chemxe
                         @:ALLOCATE(flux_src_n(i)%vf(l)%sf( &
-                                 & ix%beg:ix%end, &
-                                 & iy%beg:iy%end, &
-                                 & iz%beg:iz%end))
+                                 & idwbuff(1)%beg:idwbuff(1)%end, &
+                                 & idwbuff(2)%beg:idwbuff(2)%end, &
+                                 & idwbuff(3)%beg:idwbuff(3)%end))
                     end do
                 end if
 
             else
                 do l = 1, sys_size
                     @:ALLOCATE(flux_gsrc_n(i)%vf(l)%sf( &
-                        ix%beg:ix%end, &
-                        iy%beg:iy%end, &
-                        iz%beg:iz%end))
+                        idwbuff(1)%beg:idwbuff(1)%end, &
+                        idwbuff(2)%beg:idwbuff(2)%end, &
+                        idwbuff(3)%beg:idwbuff(3)%end))
                 end do
             end if
 
@@ -761,22 +733,13 @@ contains
 
         call nvtxStartRange("Compute_RHS")
 
-        ! Configuring Coordinate Direction Indexes =========================
-        ix%beg = -buff_size; iy%beg = 0; iz%beg = 0
-
-        if (n > 0) iy%beg = -buff_size; if (p > 0) iz%beg = -buff_size
-
-        ix%end = m - ix%beg; iy%end = n - iy%beg; iz%end = p - iz%beg
-        ! ==================================================================
-
-        !$acc update device(ix, iy, iz)
         call cpu_time(t_start)
         ! Association/Population of Working Variables ======================
         !$acc parallel loop collapse(4) gang vector default(present)
         do i = 1, sys_size
-            do l = iz%beg, iz%end
-                do k = iy%beg, iy%end
-                    do j = ix%beg, ix%end
+            do l = idwbuff(3)%beg, idwbuff(3)%end
+                do k = idwbuff(2)%beg, idwbuff(2)%end
+                    do j = idwbuff(1)%beg, idwbuff(1)%end
                         q_cons_qp%vf(i)%sf(j, k, l) = q_cons_vf(i)%sf(j, k, l)
                     end do
                 end do
@@ -789,9 +752,9 @@ contains
 
         if (mpp_lim .and. bubbles) then
             !$acc parallel loop collapse(3) gang vector default(present)
-            do l = iz%beg, iz%end
-                do k = iy%beg, iy%end
-                    do j = ix%beg, ix%end
+            do l = idwbuff(3)%beg, idwbuff(3)%end
+                do k = idwbuff(2)%beg, idwbuff(2)%end
+                    do j = idwbuff(1)%beg, idwbuff(1)%end
                         alf_sum%sf(j, k, l) = 0d0
                         !$acc loop seq
                         do i = advxb, advxe - 1
@@ -811,8 +774,8 @@ contains
         call s_convert_conservative_to_primitive_variables( &
             q_cons_qp%vf, &
             q_prim_qp%vf, &
-            gm_alpha_qp%vf, &
-            ibounds_interior(1), ibounds_interior(2), ibounds_interior(3))
+            idwint, &
+            gm_alpha_qp%vf)
         call nvtxEndRange
 
         call nvtxStartRange("RHS-MPI")
@@ -827,7 +790,7 @@ contains
 
         ! ==================================================================
 
-        if (qbmm) call s_mom_inv(q_cons_qp%vf, q_prim_qp%vf, mom_sp, mom_3d, pb, rhs_pb, mv, rhs_mv, ix, iy, iz, nbub)
+        if (qbmm) call s_mom_inv(q_cons_qp%vf, q_prim_qp%vf, mom_sp, mom_3d, pb, rhs_pb, mv, rhs_mv, idwbuff(1), idwbuff(2), idwbuff(3), nbub)
 
         call nvtxStartRange("Viscous")
         if (any(Re_size > 0)) call s_get_viscous(qL_rsx_vf, qL_rsy_vf, qL_rsz_vf, &
@@ -838,24 +801,16 @@ contains
                                                  qR_prim, &
                                                  q_prim_qp, &
                                                  dq_prim_dx_qp, dq_prim_dy_qp, dq_prim_dz_qp, &
-                                                 ix, iy, iz)
+                                                 idwbuff(1), idwbuff(2), idwbuff(3))
         call nvtxEndRange
 
         call nvtxStartRange("Surface_Tension")
         if (.not. f_is_default(sigma)) call s_get_capilary(q_prim_qp%vf)
         call nvtxEndRange
-
         ! Dimensional Splitting Loop =======================================
 
         do id = 1, num_dims
 
-            ! Configuring Coordinate Direction Indexes ======================
-            ix%beg = -buff_size; iy%beg = 0; iz%beg = 0
-
-            if (n > 0) iy%beg = -buff_size; if (p > 0) iz%beg = -buff_size
-
-            ix%end = m - ix%beg; iy%end = n - iy%beg; iz%end = p - iz%beg
-            ! ===============================================================
             ! Reconstructing Primitive/Conservative Variables ===============
 
             call nvtxStartRange("RHS-WENO")
@@ -899,21 +854,21 @@ contains
                     dqL_rsx_vf, dqL_rsy_vf, dqL_rsz_vf, &
                     dqR_rsx_vf, dqR_rsy_vf, dqR_rsz_vf, &
                     id, dqL_prim_dx_n(id)%vf(iv%beg:iv%end), dqR_prim_dx_n(id)%vf(iv%beg:iv%end), &
-                    ix, iy, iz)
+                    idwbuff(1), idwbuff(2), idwbuff(3))
                 if (n > 0) then
                     call s_reconstruct_cell_boundary_values_visc_deriv( &
                         dq_prim_dy_qp(1)%vf(iv%beg:iv%end), &
                         dqL_rsx_vf, dqL_rsy_vf, dqL_rsz_vf, &
                         dqR_rsx_vf, dqR_rsy_vf, dqR_rsz_vf, &
                         id, dqL_prim_dy_n(id)%vf(iv%beg:iv%end), dqR_prim_dy_n(id)%vf(iv%beg:iv%end), &
-                        ix, iy, iz)
+                        idwbuff(1), idwbuff(2), idwbuff(3))
                     if (p > 0) then
                         call s_reconstruct_cell_boundary_values_visc_deriv( &
                             dq_prim_dz_qp(1)%vf(iv%beg:iv%end), &
                             dqL_rsx_vf, dqL_rsy_vf, dqL_rsz_vf, &
                             dqR_rsx_vf, dqR_rsy_vf, dqR_rsz_vf, &
                             id, dqL_prim_dz_n(id)%vf(iv%beg:iv%end), dqR_prim_dz_n(id)%vf(iv%beg:iv%end), &
-                            ix, iy, iz)
+                            idwbuff(1), idwbuff(2), idwbuff(3))
                     end if
                 end if
             end if
@@ -922,14 +877,15 @@ contains
 
             ! Configuring Coordinate Direction Indexes ======================
             if (id == 1) then
-                ix%beg = -1; iy%beg = 0; iz%beg = 0
+                irx%beg = -1; iry%beg = 0; irz%beg = 0
             elseif (id == 2) then
-                ix%beg = 0; iy%beg = -1; iz%beg = 0
+                irx%beg = 0; iry%beg = -1; irz%beg = 0
             else
-                ix%beg = 0; iy%beg = 0; iz%beg = -1
+                irx%beg = 0; iry%beg = 0; irz%beg = -1
             end if
-            ix%end = m; iy%end = n; iz%end = p
+            irx%end = m; iry%end = n; irz%end = p
             ! ===============================================================
+
             call nvtxStartRange("RHS_riemann_solver")
 
             ! Computing Riemann Solver Flux and Source Flux =================
@@ -948,7 +904,7 @@ contains
                                   flux_n(id)%vf, &
                                   flux_src_n(id)%vf, &
                                   flux_gsrc_n(id)%vf, &
-                                  id, ix, iy, iz)
+                                  id, irx, iry, irz)
             call nvtxEndRange
 
             ! Additional physics and source terms ==============================
@@ -978,8 +934,7 @@ contains
                                                       flux_src_n(id)%vf, &
                                                       dq_prim_dx_qp(1)%vf, &
                                                       dq_prim_dy_qp(1)%vf, &
-                                                      dq_prim_dz_qp(1)%vf, &
-                                                      ixt, iyt, izt)
+                                                      dq_prim_dz_qp(1)%vf)
             end if
             call nvtxEndRange
 
@@ -1055,25 +1010,20 @@ contains
         ! END: Additional pphysics and source terms ============================
 
         if (run_time_info .or. probe_wrt .or. ib) then
-
-            ix%beg = -buff_size; iy%beg = 0; iz%beg = 0
-            if (n > 0) iy%beg = -buff_size; 
-            if (p > 0) iz%beg = -buff_size; 
-            ix%end = m - ix%beg; iy%end = n - iy%beg; iz%end = p - iz%beg
-            !$acc update device(ix, iy, iz)
-
             !$acc parallel loop collapse(4) gang vector default(present)
             do i = 1, sys_size
-                do l = iz%beg, iz%end
-                    do k = iy%beg, iy%end
-                        do j = ix%beg, ix%end
+                do l = idwbuff(3)%beg, idwbuff(3)%end
+                    do k = idwbuff(2)%beg, idwbuff(2)%end
+                        do j = idwbuff(1)%beg, idwbuff(1)%end
                             q_prim_vf(i)%sf(j, k, l) = q_prim_qp%vf(i)%sf(j, k, l)
                         end do
                     end do
                 end do
             end do
         end if
+
         call cpu_time(t_finish)
+
         if (t_step >= 4) then
             time_avg = (abs(t_finish - t_start) + (t_step - 4)*time_avg)/(t_step - 3)
         else
@@ -1082,6 +1032,7 @@ contains
         ! ==================================================================
 
         call nvtxEndRange
+
     end subroutine s_compute_rhs
 
     subroutine s_compute_advection_source_term(idir, rhs_vf, q_cons_vf, q_prim_vf, flux_src_n_vf)
@@ -1122,12 +1073,12 @@ contains
 
             if (bc_x%beg <= -5 .and. bc_x%beg >= -13) then
                 call s_cbc(q_prim_vf%vf, flux_n(idir)%vf, &
-                           flux_src_n(idir)%vf, idir, -1, ix, iy, iz)
+                           flux_src_n(idir)%vf, idir, -1, irx, iry, irz)
             end if
 
             if (bc_x%end <= -5 .and. bc_x%end >= -13) then
                 call s_cbc(q_prim_vf%vf, flux_n(idir)%vf, &
-                           flux_src_n(idir)%vf, idir, 1, ix, iy, iz)
+                           flux_src_n(idir)%vf, idir, 1, irx, iry, irz)
             end if
 
             !$acc parallel loop collapse(4) gang vector default(present)
@@ -1231,12 +1182,12 @@ contains
 
             if (bc_y%beg <= -5 .and. bc_y%beg >= -13) then
                 call s_cbc(q_prim_vf%vf, flux_n(idir)%vf, &
-                           flux_src_n(idir)%vf, idir, -1, ix, iy, iz)
+                           flux_src_n(idir)%vf, idir, -1, irx, iry, irz)
             end if
 
             if (bc_y%end <= -5 .and. bc_y%end >= -13) then
                 call s_cbc(q_prim_vf%vf, flux_n(idir)%vf, &
-                           flux_src_n(idir)%vf, idir, 1, ix, iy, iz)
+                           flux_src_n(idir)%vf, idir, 1, irx, iry, irz)
             end if
 
             !$acc parallel loop collapse(4) gang vector default(present)
@@ -1405,12 +1356,12 @@ contains
 
             if (bc_z%beg <= -5 .and. bc_z%beg >= -13) then
                 call s_cbc(q_prim_vf%vf, flux_n(idir)%vf, &
-                           flux_src_n(idir)%vf, idir, -1, ix, iy, iz)
+                           flux_src_n(idir)%vf, idir, -1, irx, iry, irz)
             end if
 
             if (bc_z%end <= -5 .and. bc_z%end >= -13) then
                 call s_cbc(q_prim_vf%vf, flux_n(idir)%vf, &
-                           flux_src_n(idir)%vf, idir, 1, ix, iy, iz)
+                           flux_src_n(idir)%vf, idir, 1, irx, iry, irz)
             end if
 
             if (grid_geometry == 3) then ! Cylindrical Coordinates
@@ -1640,14 +1591,13 @@ contains
     end subroutine s_compute_advection_source_term
 
     subroutine s_compute_additional_physics_rhs(idir, q_prim_vf, rhs_vf, flux_src_n, &
-                                                dq_prim_dx_vf, dq_prim_dy_vf, dq_prim_dz_vf, ixt, iyt, izt)
+                                                dq_prim_dx_vf, dq_prim_dy_vf, dq_prim_dz_vf)
 
         integer, intent(in) :: idir
         type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
         type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf
         type(scalar_field), dimension(sys_size), intent(in) :: flux_src_n
         type(scalar_field), dimension(sys_size), intent(in) :: dq_prim_dx_vf, dq_prim_dy_vf, dq_prim_dz_vf
-        type(int_bounds_info), intent(in) :: ixt, iyt, izt
 
         integer :: i, j, k, l, q
 
@@ -1708,14 +1658,14 @@ contains
                                                              dq_prim_dy_vf(mom_idx%beg:mom_idx%end), &
                                                              dq_prim_dz_vf(mom_idx%beg:mom_idx%end), &
                                                              tau_Re_vf, &
-                                                             ixt, iyt, izt)
+                                                             idwbuff(1), idwbuff(2), idwbuff(3))
                     else
                         call s_compute_viscous_stress_tensor(q_prim_vf, &
                                                              dq_prim_dx_vf(mom_idx%beg:mom_idx%end), &
                                                              dq_prim_dy_vf(mom_idx%beg:mom_idx%end), &
                                                              dq_prim_dy_vf(mom_idx%beg:mom_idx%end), &
                                                              tau_Re_vf, &
-                                                             ixt, iyt, izt)
+                                                             idwbuff(1), idwbuff(2), idwbuff(3))
                     end if
 
                     !$acc parallel loop collapse(2) gang vector default(present)
@@ -2149,17 +2099,17 @@ contains
         ! Reconstruction in s1-direction ===================================
 
         if (norm_dir == 1) then
-            is1 = ix; is2 = iy; is3 = iz
+            is1 = idwbuff(1); is2 = idwbuff(2); is3 = idwbuff(3)
             weno_dir = 1; is1%beg = is1%beg + weno_polyn
             is1%end = is1%end - weno_polyn
 
         elseif (norm_dir == 2) then
-            is1 = iy; is2 = ix; is3 = iz
+            is1 = idwbuff(2); is2 = idwbuff(1); is3 = idwbuff(3)
             weno_dir = 2; is1%beg = is1%beg + weno_polyn
             is1%end = is1%end - weno_polyn
 
         else
-            is1 = iz; is2 = iy; is3 = ix
+            is1 = idwbuff(3); is2 = idwbuff(2); is3 = idwbuff(1)
             weno_dir = 3; is1%beg = is1%beg + weno_polyn
             is1%end = is1%end - weno_polyn
 
@@ -2203,17 +2153,17 @@ contains
         ! Reconstruction in s1-direction ===================================
 
         if (norm_dir == 1) then
-            is1 = ix; is2 = iy; is3 = iz
+            is1 = idwbuff(1); is2 = idwbuff(2); is3 = idwbuff(3)
             recon_dir = 1; is1%beg = is1%beg + weno_polyn
             is1%end = is1%end - weno_polyn
 
         elseif (norm_dir == 2) then
-            is1 = iy; is2 = ix; is3 = iz
+            is1 = idwbuff(2); is2 = idwbuff(1); is3 = idwbuff(3)
             recon_dir = 2; is1%beg = is1%beg + weno_polyn
             is1%end = is1%end - weno_polyn
 
         else
-            is1 = iz; is2 = iy; is3 = ix
+            is1 = idwbuff(3); is2 = idwbuff(2); is3 = idwbuff(1)
             recon_dir = 3; is1%beg = is1%beg + weno_polyn
             is1%end = is1%end - weno_polyn
 
@@ -2310,8 +2260,8 @@ contains
         end if
 
         if (mpp_lim .and. bubbles) then
-            !deallocate(alf_sum%sf(ix%beg:ix%end, iy%beg:iy%end, iz%beg:iz%end))
-            !$acc exit data delete(alf_sum%sf(ix%beg:ix%end, iy%beg:iy%end, iz%beg:iz%end))
+            !$acc exit data delete(alf_sum%sf)
+            deallocate (alf_sum%sf)
         end if
 
         if (any(Re_size > 0)) then
diff --git a/src/simulation/m_surface_tension.fpp b/src/simulation/m_surface_tension.fpp
index 2e9ca8fed..f4dfb4bb0 100644
--- a/src/simulation/m_surface_tension.fpp
+++ b/src/simulation/m_surface_tension.fpp
@@ -50,8 +50,8 @@ module m_surface_tension
     !$acc declare create(gL_x, gR_x, gL_y, gR_y, gL_z, gR_z)
 #endif
 
-    type(int_bounds_info) :: ix, iy, iz, is1, is2, is3, iv
-    !$acc declare create(ix, iy, iz, is1, is2, is3, iv)
+    type(int_bounds_info) :: is1, is2, is3, iv
+    !$acc declare create(is1, is2, is3, iv)
 
     integer :: j, k, l, i
 
@@ -59,30 +59,22 @@ contains
 
     subroutine s_initialize_surface_tension_module
 
-        ! Configuring Coordinate Direction Indexes =========================
-        ix%beg = -buff_size; iy%beg = 0; iz%beg = 0
-
-        if (n > 0) iy%beg = -buff_size; if (p > 0) iz%beg = -buff_size
-
-        ix%end = m - ix%beg; iy%end = n - iy%beg; iz%end = p - iz%beg
-        ! ==================================================================
-
         @:ALLOCATE_GLOBAL(c_divs(1:num_dims + 1))
 
         do j = 1, num_dims + 1
-            @:ALLOCATE(c_divs(j)%sf(ix%beg:ix%end, iy%beg:iy%end, iz%beg:iz%end))
+            @:ALLOCATE(c_divs(j)%sf(idwbuff(1)%beg:idwbuff(1)%end, idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end))
             @:ACC_SETUP_SFs(c_divs(j))
         end do
 
-        @:ALLOCATE_GLOBAL(gL_x(ix%beg:ix%end, iy%beg:iy%end, iz%beg:iz%end, num_dims + 1))
-        @:ALLOCATE_GLOBAL(gR_x(ix%beg:ix%end, iy%beg:iy%end, iz%beg:iz%end, num_dims + 1))
+        @:ALLOCATE_GLOBAL(gL_x(idwbuff(1)%beg:idwbuff(1)%end, idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end, num_dims + 1))
+        @:ALLOCATE_GLOBAL(gR_x(idwbuff(1)%beg:idwbuff(1)%end, idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end, num_dims + 1))
 
-        @:ALLOCATE_GLOBAL(gL_y(iy%beg:iy%end, ix%beg:ix%end, iz%beg:iz%end, num_dims + 1))
-        @:ALLOCATE_GLOBAL(gR_y(iy%beg:iy%end, ix%beg:ix%end, iz%beg:iz%end, num_dims + 1))
+        @:ALLOCATE_GLOBAL(gL_y(idwbuff(2)%beg:idwbuff(2)%end, idwbuff(1)%beg:idwbuff(1)%end, idwbuff(3)%beg:idwbuff(3)%end, num_dims + 1))
+        @:ALLOCATE_GLOBAL(gR_y(idwbuff(2)%beg:idwbuff(2)%end, idwbuff(1)%beg:idwbuff(1)%end, idwbuff(3)%beg:idwbuff(3)%end, num_dims + 1))
 
         if (p > 0) then
-            @:ALLOCATE_GLOBAL(gL_z(iz%beg:iz%end, ix%beg:ix%end, iy%beg:iy%end, num_dims + 1))
-            @:ALLOCATE_GLOBAL(gR_z(iz%beg:iz%end, ix%beg:ix%end, iy%beg:iy%end, num_dims + 1))
+            @:ALLOCATE_GLOBAL(gL_z(idwbuff(3)%beg:idwbuff(3)%end, idwbuff(1)%beg:idwbuff(1)%end, idwbuff(2)%beg:idwbuff(2)%end, num_dims + 1))
+            @:ALLOCATE_GLOBAL(gR_z(idwbuff(3)%beg:idwbuff(3)%end, idwbuff(1)%beg:idwbuff(1)%end, idwbuff(2)%beg:idwbuff(2)%end, num_dims + 1))
         end if
     end subroutine s_initialize_surface_tension_module
 
@@ -335,17 +327,17 @@ contains
         ! Reconstruction in s1-direction ===================================
 
         if (norm_dir == 1) then
-            is1 = ix; is2 = iy; is3 = iz
+            is1 = idwbuff(1); is2 = idwbuff(2); is3 = idwbuff(3)
             recon_dir = 1; is1%beg = is1%beg + weno_polyn
             is1%end = is1%end - weno_polyn
 
         elseif (norm_dir == 2) then
-            is1 = iy; is2 = ix; is3 = iz
+            is1 = idwbuff(2); is2 = idwbuff(1); is3 = idwbuff(3)
             recon_dir = 2; is1%beg = is1%beg + weno_polyn
             is1%end = is1%end - weno_polyn
 
         else
-            is1 = iz; is2 = iy; is3 = ix
+            is1 = idwbuff(3); is2 = idwbuff(2); is3 = idwbuff(1)
             recon_dir = 3; is1%beg = is1%beg + weno_polyn
             is1%end = is1%end - weno_polyn
 
diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp
index 8f0131ef0..4c8aa126d 100644
--- a/src/simulation/m_time_steppers.fpp
+++ b/src/simulation/m_time_steppers.fpp
@@ -99,9 +99,6 @@ contains
         !!      other procedures that are necessary to setup the module.
     subroutine s_initialize_time_steppers_module
 
-        type(int_bounds_info) :: ix_t, iy_t, iz_t !<
-            !! Indical bounds in the x-, y- and z-directions
-
         integer :: i, j !< Generic loop iterators
 
         ! Setting number of time-stages for selected time-stepping scheme
@@ -111,22 +108,6 @@ contains
             num_ts = 2
         end if
 
-        ! Setting the indical bounds in the x-, y- and z-directions
-        ix_t%beg = -buff_size; ix_t%end = m + buff_size
-
-        if (n > 0) then
-            iy_t%beg = -buff_size; iy_t%end = n + buff_size
-
-            if (p > 0) then
-                iz_t%beg = -buff_size; iz_t%end = p + buff_size
-            else
-                iz_t%beg = 0; iz_t%end = 0
-            end if
-        else
-            iy_t%beg = 0; iy_t%end = 0
-            iz_t%beg = 0; iz_t%end = 0
-        end if
-
         ! Allocating the cell-average conservative variables
         @:ALLOCATE_GLOBAL(q_cons_ts(1:num_ts))
 
@@ -136,9 +117,9 @@ contains
 
         do i = 1, num_ts
             do j = 1, sys_size
-                @:ALLOCATE(q_cons_ts(i)%vf(j)%sf(ix_t%beg:ix_t%end, &
-                    iy_t%beg:iy_t%end, &
-                    iz_t%beg:iz_t%end))
+                @:ALLOCATE(q_cons_ts(i)%vf(j)%sf(idwbuff(1)%beg:idwbuff(1)%end, &
+                    idwbuff(2)%beg:idwbuff(2)%end, &
+                    idwbuff(3)%beg:idwbuff(3)%end))
             end do
             @:ACC_SETUP_VFs(q_cons_ts(i))
         end do
@@ -153,9 +134,9 @@ contains
 
             do i = 0, 3
                 do j = 1, sys_size
-                    @:ALLOCATE(q_prim_ts(i)%vf(j)%sf(ix_t%beg:ix_t%end, &
-                        iy_t%beg:iy_t%end, &
-                        iz_t%beg:iz_t%end))
+                    @:ALLOCATE(q_prim_ts(i)%vf(j)%sf(idwbuff(1)%beg:idwbuff(1)%end, &
+                        idwbuff(2)%beg:idwbuff(2)%end, &
+                        idwbuff(3)%beg:idwbuff(3)%end))
                 end do
             end do
 
@@ -168,23 +149,23 @@ contains
         @:ALLOCATE_GLOBAL(q_prim_vf(1:sys_size))
 
         do i = 1, adv_idx%end
-            @:ALLOCATE(q_prim_vf(i)%sf(ix_t%beg:ix_t%end, &
-                iy_t%beg:iy_t%end, &
-                iz_t%beg:iz_t%end))
+            @:ALLOCATE(q_prim_vf(i)%sf(idwbuff(1)%beg:idwbuff(1)%end, &
+                idwbuff(2)%beg:idwbuff(2)%end, &
+                idwbuff(3)%beg:idwbuff(3)%end))
             @:ACC_SETUP_SFs(q_prim_vf(i))
         end do
 
         if (bubbles) then
             do i = bub_idx%beg, bub_idx%end
-                @:ALLOCATE(q_prim_vf(i)%sf(ix_t%beg:ix_t%end, &
-                    iy_t%beg:iy_t%end, &
-                    iz_t%beg:iz_t%end))
+                @:ALLOCATE(q_prim_vf(i)%sf(idwbuff(1)%beg:idwbuff(1)%end, &
+                    idwbuff(2)%beg:idwbuff(2)%end, &
+                    idwbuff(3)%beg:idwbuff(3)%end))
                 @:ACC_SETUP_SFs(q_prim_vf(i))
             end do
             if (adv_n) then
-                @:ALLOCATE(q_prim_vf(n_idx)%sf(ix_t%beg:ix_t%end, &
-                    iy_t%beg:iy_t%end, &
-                    iz_t%beg:iz_t%end))
+                @:ALLOCATE(q_prim_vf(n_idx)%sf(idwbuff(1)%beg:idwbuff(1)%end, &
+                    idwbuff(2)%beg:idwbuff(2)%end, &
+                    idwbuff(3)%beg:idwbuff(3)%end))
                 @:ACC_SETUP_SFs(q_prim_vf(n_idx))
             end if
         end if
@@ -192,106 +173,106 @@ contains
         if (hypoelasticity) then
 
             do i = stress_idx%beg, stress_idx%end
-                @:ALLOCATE(q_prim_vf(i)%sf(ix_t%beg:ix_t%end, &
-                    iy_t%beg:iy_t%end, &
-                    iz_t%beg:iz_t%end))
+                @:ALLOCATE(q_prim_vf(i)%sf(idwbuff(1)%beg:idwbuff(1)%end, &
+                    idwbuff(2)%beg:idwbuff(2)%end, &
+                    idwbuff(3)%beg:idwbuff(3)%end))
                 @:ACC_SETUP_SFs(q_prim_vf(i))
             end do
         end if
 
         if (model_eqns == 3) then
             do i = internalEnergies_idx%beg, internalEnergies_idx%end
-                @:ALLOCATE(q_prim_vf(i)%sf(ix_t%beg:ix_t%end, &
-                    iy_t%beg:iy_t%end, &
-                    iz_t%beg:iz_t%end))
+                @:ALLOCATE(q_prim_vf(i)%sf(idwbuff(1)%beg:idwbuff(1)%end, &
+                    idwbuff(2)%beg:idwbuff(2)%end, &
+                    idwbuff(3)%beg:idwbuff(3)%end))
                 @:ACC_SETUP_SFs(q_prim_vf(i))
             end do
         end if
 
         if (.not. f_is_default(sigma)) then
-            @:ALLOCATE(q_prim_vf(c_idx)%sf(ix_t%beg:ix_t%end, &
-                iy_t%beg:iy_t%end, &
-                iz_t%beg:iz_t%end))
+            @:ALLOCATE(q_prim_vf(c_idx)%sf(idwbuff(1)%beg:idwbuff(1)%end, &
+                idwbuff(2)%beg:idwbuff(2)%end, &
+                idwbuff(3)%beg:idwbuff(3)%end))
             @:ACC_SETUP_SFs(q_prim_vf(c_idx))
         end if
 
         if (chemistry) then
             do i = chemxb, chemxe
-                @:ALLOCATE(q_prim_vf(i)%sf(ix_t%beg:ix_t%end, &
-                    iy_t%beg:iy_t%end, &
-                    iz_t%beg:iz_t%end))
+                @:ALLOCATE(q_prim_vf(i)%sf(idwbuff(1)%beg:idwbuff(1)%end, &
+                    idwbuff(2)%beg:idwbuff(2)%end, &
+                    idwbuff(3)%beg:idwbuff(3)%end))
                 @:ACC_SETUP_SFs(q_prim_vf(i))
             end do
 
-            @:ALLOCATE(q_prim_vf(tempxb)%sf(ix_t%beg:ix_t%end, &
-                iy_t%beg:iy_t%end, &
-                iz_t%beg:iz_t%end))
+            @:ALLOCATE(q_prim_vf(tempxb)%sf(idwbuff(1)%beg:idwbuff(1)%end, &
+                idwbuff(2)%beg:idwbuff(2)%end, &
+                idwbuff(3)%beg:idwbuff(3)%end))
             @:ACC_SETUP_SFs(q_prim_vf(tempxb))
         end if
 
         @:ALLOCATE_GLOBAL(pb_ts(1:2))
         !Initialize bubble variables pb and mv at all quadrature nodes for all R0 bins
         if (qbmm .and. (.not. polytropic)) then
-            @:ALLOCATE(pb_ts(1)%sf(ix_t%beg:ix_t%end, &
-                iy_t%beg:iy_t%end, &
-                iz_t%beg:iz_t%end, 1:nnode, 1:nb))
+            @:ALLOCATE(pb_ts(1)%sf(idwbuff(1)%beg:idwbuff(1)%end, &
+                idwbuff(2)%beg:idwbuff(2)%end, &
+                idwbuff(3)%beg:idwbuff(3)%end, 1:nnode, 1:nb))
             @:ACC_SETUP_SFs(pb_ts(1))
 
-            @:ALLOCATE(pb_ts(2)%sf(ix_t%beg:ix_t%end, &
-                iy_t%beg:iy_t%end, &
-                iz_t%beg:iz_t%end, 1:nnode, 1:nb))
+            @:ALLOCATE(pb_ts(2)%sf(idwbuff(1)%beg:idwbuff(1)%end, &
+                idwbuff(2)%beg:idwbuff(2)%end, &
+                idwbuff(3)%beg:idwbuff(3)%end, 1:nnode, 1:nb))
             @:ACC_SETUP_SFs(pb_ts(2))
 
-            @:ALLOCATE_GLOBAL(rhs_pb(ix_t%beg:ix_t%end, &
-                iy_t%beg:iy_t%end, &
-                iz_t%beg:iz_t%end, 1:nnode, 1:nb))
+            @:ALLOCATE_GLOBAL(rhs_pb(idwbuff(1)%beg:idwbuff(1)%end, &
+                idwbuff(2)%beg:idwbuff(2)%end, &
+                idwbuff(3)%beg:idwbuff(3)%end, 1:nnode, 1:nb))
         else if (qbmm .and. polytropic) then
-            @:ALLOCATE(pb_ts(1)%sf(ix_t%beg:ix_t%beg + 1, &
-                iy_t%beg:iy_t%beg + 1, &
-                iz_t%beg:iz_t%beg + 1, 1:nnode, 1:nb))
+            @:ALLOCATE(pb_ts(1)%sf(idwbuff(1)%beg:idwbuff(1)%beg + 1, &
+                idwbuff(2)%beg:idwbuff(2)%beg + 1, &
+                idwbuff(3)%beg:idwbuff(3)%beg + 1, 1:nnode, 1:nb))
             @:ACC_SETUP_SFs(pb_ts(1))
 
-            @:ALLOCATE(pb_ts(2)%sf(ix_t%beg:ix_t%beg + 1, &
-                iy_t%beg:iy_t%beg + 1, &
-                iz_t%beg:iz_t%beg + 1, 1:nnode, 1:nb))
+            @:ALLOCATE(pb_ts(2)%sf(idwbuff(1)%beg:idwbuff(1)%beg + 1, &
+                idwbuff(2)%beg:idwbuff(2)%beg + 1, &
+                idwbuff(3)%beg:idwbuff(3)%beg + 1, 1:nnode, 1:nb))
             @:ACC_SETUP_SFs(pb_ts(2))
 
-            @:ALLOCATE_GLOBAL(rhs_pb(ix_t%beg:ix_t%beg + 1, &
-                iy_t%beg:iy_t%beg + 1, &
-                iz_t%beg:iz_t%beg + 1, 1:nnode, 1:nb))
+            @:ALLOCATE_GLOBAL(rhs_pb(idwbuff(1)%beg:idwbuff(1)%beg + 1, &
+                idwbuff(2)%beg:idwbuff(2)%beg + 1, &
+                idwbuff(3)%beg:idwbuff(3)%beg + 1, 1:nnode, 1:nb))
         end if
 
         @:ALLOCATE_GLOBAL(mv_ts(1:2))
 
         if (qbmm .and. (.not. polytropic)) then
-            @:ALLOCATE(mv_ts(1)%sf(ix_t%beg:ix_t%end, &
-                iy_t%beg:iy_t%end, &
-                iz_t%beg:iz_t%end, 1:nnode, 1:nb))
+            @:ALLOCATE(mv_ts(1)%sf(idwbuff(1)%beg:idwbuff(1)%end, &
+                idwbuff(2)%beg:idwbuff(2)%end, &
+                idwbuff(3)%beg:idwbuff(3)%end, 1:nnode, 1:nb))
             @:ACC_SETUP_SFs(mv_ts(1))
 
-            @:ALLOCATE(mv_ts(2)%sf(ix_t%beg:ix_t%end, &
-                iy_t%beg:iy_t%end, &
-                iz_t%beg:iz_t%end, 1:nnode, 1:nb))
+            @:ALLOCATE(mv_ts(2)%sf(idwbuff(1)%beg:idwbuff(1)%end, &
+                idwbuff(2)%beg:idwbuff(2)%end, &
+                idwbuff(3)%beg:idwbuff(3)%end, 1:nnode, 1:nb))
             @:ACC_SETUP_SFs(mv_ts(2))
 
-            @:ALLOCATE_GLOBAL(rhs_mv(ix_t%beg:ix_t%end, &
-                iy_t%beg:iy_t%end, &
-                iz_t%beg:iz_t%end, 1:nnode, 1:nb))
+            @:ALLOCATE_GLOBAL(rhs_mv(idwbuff(1)%beg:idwbuff(1)%end, &
+                idwbuff(2)%beg:idwbuff(2)%end, &
+                idwbuff(3)%beg:idwbuff(3)%end, 1:nnode, 1:nb))
 
         else if (qbmm .and. polytropic) then
-            @:ALLOCATE(mv_ts(1)%sf(ix_t%beg:ix_t%beg + 1, &
-                iy_t%beg:iy_t%beg + 1, &
-                iz_t%beg:iz_t%beg + 1, 1:nnode, 1:nb))
+            @:ALLOCATE(mv_ts(1)%sf(idwbuff(1)%beg:idwbuff(1)%beg + 1, &
+                idwbuff(2)%beg:idwbuff(2)%beg + 1, &
+                idwbuff(3)%beg:idwbuff(3)%beg + 1, 1:nnode, 1:nb))
             @:ACC_SETUP_SFs(mv_ts(1))
 
-            @:ALLOCATE(mv_ts(2)%sf(ix_t%beg:ix_t%beg + 1, &
-                iy_t%beg:iy_t%beg + 1, &
-                iz_t%beg:iz_t%beg + 1, 1:nnode, 1:nb))
+            @:ALLOCATE(mv_ts(2)%sf(idwbuff(1)%beg:idwbuff(1)%beg + 1, &
+                idwbuff(2)%beg:idwbuff(2)%beg + 1, &
+                idwbuff(3)%beg:idwbuff(3)%beg + 1, 1:nnode, 1:nb))
             @:ACC_SETUP_SFs(mv_ts(2))
 
-            @:ALLOCATE_GLOBAL(rhs_mv(ix_t%beg:ix_t%beg + 1, &
-                iy_t%beg:iy_t%beg + 1, &
-                iz_t%beg:iz_t%beg + 1, 1:nnode, 1:nb))
+            @:ALLOCATE_GLOBAL(rhs_mv(idwbuff(1)%beg:idwbuff(1)%beg + 1, &
+                idwbuff(2)%beg:idwbuff(2)%beg + 1, &
+                idwbuff(3)%beg:idwbuff(3)%beg + 1, 1:nnode, 1:nb))
         end if
 
         ! Allocating the cell-average RHS variables
@@ -913,18 +894,15 @@ contains
 
         integer, intent(in) :: t_step
 
-        type(int_bounds_info) :: ix, iy, iz
         type(vector_field) :: gm_alpha_qp
 
         integer :: i, j, k, l, q !< Generic loop iterator
 
-        ix%beg = 0; iy%beg = 0; iz%beg = 0
-        ix%end = m; iy%end = n; iz%end = p
         call s_convert_conservative_to_primitive_variables( &
             q_cons_ts(1)%vf, &
             q_prim_vf, &
-            gm_alpha_qp%vf, &
-            ix, iy, iz)
+            idwint, &
+            gm_alpha_qp%vf)
 
         call s_compute_bubble_source(q_cons_ts(1)%vf, q_prim_vf, t_step, rhs_vf)
 
@@ -946,17 +924,13 @@ contains
         real(kind(0d0)), dimension(2) :: Re         !< Cell-avg. Reynolds numbers
         type(vector_field) :: gm_alpha_qp
         real(kind(0d0)) :: dt_local
-        type(int_bounds_info) :: ix, iy, iz
         integer :: i, j, k, l, q !< Generic loop iterators
 
-        ix%beg = 0; iy%beg = 0; iz%beg = 0
-        ix%end = m; iy%end = n; iz%end = p
-
         call s_convert_conservative_to_primitive_variables( &
             q_cons_ts(1)%vf, &
             q_prim_vf, &
-            gm_alpha_qp%vf, &
-            ix, iy, iz)
+            idwint, &
+            gm_alpha_qp%vf)
 
         !$acc parallel loop collapse(3) gang vector default(present) private(vel, alpha, Re)
         do l = 0, p
diff --git a/src/simulation/m_viscous.fpp b/src/simulation/m_viscous.fpp
index c8c83f6b4..139fccc7c 100644
--- a/src/simulation/m_viscous.fpp
+++ b/src/simulation/m_viscous.fpp
@@ -41,15 +41,6 @@ contains
     subroutine s_initialize_viscous_module
 
         integer :: i, j !< generic loop iterators
-        type(int_bounds_info) :: ix, iy, iz
-
-        ! Configuring Coordinate Direction Indexes =========================
-        ix%beg = -buff_size; iy%beg = 0; iz%beg = 0
-
-        if (n > 0) iy%beg = -buff_size; if (p > 0) iz%beg = -buff_size
-
-        ix%end = m - ix%beg; iy%end = n - iy%beg; iz%end = p - iz%beg
-        ! ==================================================================
 
         @:ALLOCATE_GLOBAL(Res_viscous(1:2, 1:maxval(Re_size)))
 
@@ -553,7 +544,7 @@ contains
                              dqL_prim_dz_n, dqR_prim_dz_n
 
         type(vector_field), dimension(1), intent(inout) :: dq_prim_dx_qp, dq_prim_dy_qp, dq_prim_dz_qp
-        type(int_bounds_info), intent(inout) :: ix, iy, iz
+        type(int_bounds_info), intent(in) :: ix, iy, iz
 
         integer :: i, j, k, l
 
@@ -949,8 +940,7 @@ contains
                         call s_compute_fd_gradient(q_prim_qp%vf(i), &
                                                    dq_prim_dx_qp(1)%vf(i), &
                                                    dq_prim_dy_qp(1)%vf(i), &
-                                                   dq_prim_dz_qp(1)%vf(i), &
-                                                   ix, iy, iz, buff_size)
+                                                   dq_prim_dz_qp(1)%vf(i))
                     end do
 
                 else
@@ -959,8 +949,7 @@ contains
                         call s_compute_fd_gradient(q_prim_qp%vf(i), &
                                                    dq_prim_dx_qp(1)%vf(i), &
                                                    dq_prim_dy_qp(1)%vf(i), &
-                                                   dq_prim_dy_qp(1)%vf(i), &
-                                                   ix, iy, iz, buff_size)
+                                                   dq_prim_dy_qp(1)%vf(i))
                     end do
 
                 end if
@@ -971,8 +960,7 @@ contains
                     call s_compute_fd_gradient(q_prim_qp%vf(i), &
                                                dq_prim_dx_qp(1)%vf(i), &
                                                dq_prim_dx_qp(1)%vf(i), &
-                                               dq_prim_dx_qp(1)%vf(i), &
-                                               ix, iy, iz, buff_size)
+                                               dq_prim_dx_qp(1)%vf(i))
                 end do
 
             end if
@@ -1312,29 +1300,27 @@ contains
         !!  @param grad_y Second coordinate direction component of the derivative
         !!  @param grad_z Third coordinate direction component of the derivative
         !!  @param norm Norm of the gradient vector
-    subroutine s_compute_fd_gradient(var, grad_x, grad_y, grad_z, &
-                                     ix, iy, iz, buff_size_in)
+    subroutine s_compute_fd_gradient(var, grad_x, grad_y, grad_z)
 
         type(scalar_field), intent(in) :: var
         type(scalar_field), intent(inout) :: grad_x
         type(scalar_field), intent(inout) :: grad_y
         type(scalar_field), intent(inout) :: grad_z
-        type(int_bounds_info), intent(inout) :: ix, iy, iz
-        integer, intent(in) :: buff_size_in
+        type(int_bounds_info) :: ix, iy, iz
 
         integer :: j, k, l !< Generic loop iterators
 
-        ix%beg = -buff_size_in; ix%end = m + buff_size_in; 
+        ix%beg = 1 - buff_size; ix%end = m + buff_size - 1
         if (n > 0) then
-            iy%beg = -buff_size_in; iy%end = n + buff_size_in
+            iy%beg = 1 - buff_size; iy%end = n + buff_size - 1
         else
-            iy%beg = -1; iy%end = 1
+            iy%beg = 0; iy%end = 0
         end if
 
         if (p > 0) then
-            iz%beg = -buff_size_in; iz%end = p + buff_size_in
+            iz%beg = 1 - buff_size; iz%end = p + buff_size - 1
         else
-            iz%beg = -1; iz%end = 1
+            iz%beg = 0; iz%end = 0
         end if
 
         is1_viscous = ix; is2_viscous = iy; is3_viscous = iz
@@ -1342,9 +1328,9 @@ contains
         !$acc update device(is1_viscous, is2_viscous, is3_viscous)
 
         !$acc parallel loop collapse(3) gang vector default(present)
-        do l = is3_viscous%beg + 1, is3_viscous%end - 1
-            do k = is2_viscous%beg + 1, is2_viscous%end - 1
-                do j = is1_viscous%beg + 1, is1_viscous%end - 1
+        do l = is3_viscous%beg, is3_viscous%end
+            do k = is2_viscous%beg, is2_viscous%end
+                do j = is1_viscous%beg, is1_viscous%end
                     grad_x%sf(j, k, l) = &
                         (var%sf(j + 1, k, l) - var%sf(j - 1, k, l))/ &
                         (x_cc(j + 1) - x_cc(j - 1))
@@ -1354,9 +1340,9 @@ contains
 
         if (n > 0) then
             !$acc parallel loop collapse(3) gang vector
-            do l = is3_viscous%beg + 1, is3_viscous%end - 1
-                do k = is2_viscous%beg + 1, is2_viscous%end - 1
-                    do j = is1_viscous%beg + 1, is1_viscous%end - 1
+            do l = is3_viscous%beg, is3_viscous%end
+                do k = is2_viscous%beg, is2_viscous%end
+                    do j = is1_viscous%beg, is1_viscous%end
                         grad_y%sf(j, k, l) = &
                             (var%sf(j, k + 1, l) - var%sf(j, k - 1, l))/ &
                             (y_cc(k + 1) - y_cc(k - 1))
@@ -1367,9 +1353,9 @@ contains
 
         if (p > 0) then
             !$acc parallel loop collapse(3) gang vector
-            do l = is3_viscous%beg + 1, is3_viscous%end - 1
-                do k = is2_viscous%beg + 1, is2_viscous%end - 1
-                    do j = is1_viscous%beg + 1, is1_viscous%end - 1
+            do l = is3_viscous%beg, is3_viscous%end
+                do k = is2_viscous%beg, is2_viscous%end
+                    do j = is1_viscous%beg, is1_viscous%end
                         grad_z%sf(j, k, l) = &
                             (var%sf(j, k, l + 1) - var%sf(j, k, l - 1))/ &
                             (z_cc(l + 1) - z_cc(l - 1))
@@ -1378,54 +1364,39 @@ contains
             end do
         end if
 
-        is1_viscous%beg = -buff_size_in; is1_viscous%end = m + buff_size_in; 
-        if (n > 0) then
-            is2_viscous%beg = -buff_size_in; is2_viscous%end = n + buff_size_in
-        else
-            is2_viscous%beg = 0; is2_viscous%end = 0
-        end if
-
-        if (p > 0) then
-            is3_viscous%beg = -buff_size_in; is3_viscous%end = p + buff_size_in
-        else
-            is3_viscous%beg = 0; is3_viscous%end = 0
-        end if
-
-        !$acc update device(is1_viscous, is2_viscous, is3_viscous)
-
         !$acc parallel loop collapse(2) gang vector default(present)
-        do l = is3_viscous%beg, is3_viscous%end
-            do k = is2_viscous%beg, is2_viscous%end
-                grad_x%sf(is1_viscous%beg, k, l) = &
-                    (-3d0*var%sf(is1_viscous%beg, k, l) + 4d0*var%sf(is1_viscous%beg + 1, k, l) - var%sf(is1_viscous%beg + 2, k, l))/ &
-                    (x_cc(is1_viscous%beg + 2) - x_cc(is1_viscous%beg))
-                grad_x%sf(is1_viscous%end, k, l) = &
-                    (3d0*var%sf(is1_viscous%end, k, l) - 4d0*var%sf(is1_viscous%end - 1, k, l) + var%sf(is1_viscous%end - 2, k, l))/ &
-                    (x_cc(is1_viscous%end) - x_cc(is1_viscous%end - 2))
+        do l = idwbuff(3)%beg, idwbuff(3)%end
+            do k = idwbuff(2)%beg, idwbuff(2)%end
+                grad_x%sf(idwbuff(1)%beg, k, l) = &
+                    (-3d0*var%sf(idwbuff(1)%beg, k, l) + 4d0*var%sf(idwbuff(1)%beg + 1, k, l) - var%sf(idwbuff(1)%beg + 2, k, l))/ &
+                    (x_cc(idwbuff(1)%beg + 2) - x_cc(idwbuff(1)%beg))
+                grad_x%sf(idwbuff(1)%end, k, l) = &
+                    (+3d0*var%sf(idwbuff(1)%end, k, l) - 4d0*var%sf(idwbuff(1)%end - 1, k, l) + var%sf(idwbuff(1)%end - 2, k, l))/ &
+                    (x_cc(idwbuff(1)%end) - x_cc(idwbuff(1)%end - 2))
             end do
         end do
         if (n > 0) then
             !$acc parallel loop collapse(2) gang vector default(present)
-            do l = is3_viscous%beg, is3_viscous%end
-                do j = is1_viscous%beg, is1_viscous%end
-                    grad_y%sf(j, is2_viscous%beg, l) = &
-                        (-3d0*var%sf(j, is2_viscous%beg, l) + 4d0*var%sf(j, is2_viscous%beg + 1, l) - var%sf(j, is2_viscous%beg + 2, l))/ &
-                        (y_cc(is2_viscous%beg + 2) - y_cc(is2_viscous%beg))
-                    grad_y%sf(j, is2_viscous%end, l) = &
-                        (3d0*var%sf(j, is2_viscous%end, l) - 4d0*var%sf(j, is2_viscous%end - 1, l) + var%sf(j, is2_viscous%end - 2, l))/ &
-                        (y_cc(is2_viscous%end) - y_cc(is2_viscous%end - 2))
+            do l = idwbuff(3)%beg, idwbuff(3)%end
+                do j = idwbuff(1)%beg, idwbuff(1)%end
+                    grad_y%sf(j, idwbuff(2)%beg, l) = &
+                        (-3d0*var%sf(j, idwbuff(2)%beg, l) + 4d0*var%sf(j, idwbuff(2)%beg + 1, l) - var%sf(j, idwbuff(2)%beg + 2, l))/ &
+                        (y_cc(idwbuff(2)%beg + 2) - y_cc(idwbuff(2)%beg))
+                    grad_y%sf(j, idwbuff(2)%end, l) = &
+                        (+3d0*var%sf(j, idwbuff(2)%end, l) - 4d0*var%sf(j, idwbuff(2)%end - 1, l) + var%sf(j, idwbuff(2)%end - 2, l))/ &
+                        (y_cc(idwbuff(2)%end) - y_cc(idwbuff(2)%end - 2))
                 end do
             end do
             if (p > 0) then
                 !$acc parallel loop collapse(2) gang vector default(present)
-                do k = is2_viscous%beg, is2_viscous%end
-                    do j = is1_viscous%beg, is1_viscous%end
-                        grad_z%sf(j, k, is3_viscous%beg) = &
-                            (-3d0*var%sf(j, k, is3_viscous%beg) + 4d0*var%sf(j, k, is3_viscous%beg + 1) - var%sf(j, k, is3_viscous%beg + 2))/ &
-                            (z_cc(is3_viscous%beg + 2) - z_cc(is3_viscous%beg))
-                        grad_z%sf(j, k, is3_viscous%end) = &
-                            (3d0*var%sf(j, k, is3_viscous%end) - 4d0*var%sf(j, k, is3_viscous%end - 1) + var%sf(j, k, is3_viscous%end - 2))/ &
-                            (z_cc(is3_viscous%end) - z_cc(is3_viscous%end - 2))
+                do k = idwbuff(2)%beg, idwbuff(2)%end
+                    do j = idwbuff(1)%beg, idwbuff(1)%end
+                        grad_z%sf(j, k, idwbuff(3)%beg) = &
+                            (-3d0*var%sf(j, k, idwbuff(3)%beg) + 4d0*var%sf(j, k, idwbuff(3)%beg + 1) - var%sf(j, k, idwbuff(3)%beg + 2))/ &
+                            (z_cc(idwbuff(3)%beg + 2) - z_cc(is3_viscous%beg))
+                        grad_z%sf(j, k, idwbuff(3)%end) = &
+                            (+3d0*var%sf(j, k, idwbuff(3)%end) - 4d0*var%sf(j, k, idwbuff(3)%end - 1) + var%sf(j, k, idwbuff(3)%end - 2))/ &
+                            (z_cc(idwbuff(3)%end) - z_cc(idwbuff(3)%end - 2))
                     end do
                 end do
             end if
@@ -1433,8 +1404,8 @@ contains
 
         if (bc_x%beg <= -3) then
             !$acc parallel loop collapse(2) gang vector default(present)
-            do l = is3_viscous%beg, is3_viscous%end
-                do k = is2_viscous%beg, is2_viscous%end
+            do l = idwbuff(3)%beg, idwbuff(3)%end
+                do k = idwbuff(2)%beg, idwbuff(2)%end
                     grad_x%sf(0, k, l) = (-3d0*var%sf(0, k, l) + 4d0*var%sf(1, k, l) - var%sf(2, k, l))/ &
                                          (x_cc(2) - x_cc(0))
                 end do
@@ -1442,8 +1413,8 @@ contains
         end if
         if (bc_x%end <= -3) then
             !$acc parallel loop collapse(2) gang vector default(present)
-            do l = is3_viscous%beg, is3_viscous%end
-                do k = is2_viscous%beg, is2_viscous%end
+            do l = idwbuff(3)%beg, idwbuff(3)%end
+                do k = idwbuff(2)%beg, idwbuff(2)%end
                     grad_x%sf(m, k, l) = (3d0*var%sf(m, k, l) - 4d0*var%sf(m - 1, k, l) + var%sf(m - 2, k, l))/ &
                                          (x_cc(m) - x_cc(m - 2))
                 end do
@@ -1452,8 +1423,8 @@ contains
         if (n > 0) then
             if (bc_y%beg <= -3 .and. bc_y%beg /= -13) then
                 !$acc parallel loop collapse(2) gang vector default(present)
-                do l = is3_viscous%beg, is3_viscous%end
-                    do j = is1_viscous%beg, is1_viscous%end
+                do l = idwbuff(3)%beg, idwbuff(3)%end
+                    do j = idwbuff(1)%beg, idwbuff(1)%end
                         grad_y%sf(j, 0, l) = (-3d0*var%sf(j, 0, l) + 4d0*var%sf(j, 1, l) - var%sf(j, 2, l))/ &
                                              (y_cc(2) - y_cc(0))
                     end do
@@ -1461,8 +1432,8 @@ contains
             end if
             if (bc_y%end <= -3) then
                 !$acc parallel loop collapse(2) gang vector default(present)
-                do l = is3_viscous%beg, is3_viscous%end
-                    do j = is1_viscous%beg, is1_viscous%end
+                do l = idwbuff(3)%beg, idwbuff(3)%end
+                    do j = idwbuff(1)%beg, idwbuff(1)%end
                         grad_y%sf(j, n, l) = (3d0*var%sf(j, n, l) - 4d0*var%sf(j, n - 1, l) + var%sf(j, n - 2, l))/ &
                                              (y_cc(n) - y_cc(n - 2))
                     end do
@@ -1471,8 +1442,8 @@ contains
             if (p > 0) then
                 if (bc_z%beg <= -3) then
                     !$acc parallel loop collapse(2) gang vector default(present)
-                    do k = is2_viscous%beg, is2_viscous%end
-                        do j = is1_viscous%beg, is1_viscous%end
+                    do k = idwbuff(2)%beg, idwbuff(2)%end
+                        do j = idwbuff(1)%beg, idwbuff(1)%end
                             grad_z%sf(j, k, 0) = &
                                 (-3d0*var%sf(j, k, 0) + 4d0*var%sf(j, k, 1) - var%sf(j, k, 2))/ &
                                 (z_cc(2) - z_cc(0))
@@ -1481,8 +1452,8 @@ contains
                 end if
                 if (bc_z%end <= -3) then
                     !$acc parallel loop collapse(2) gang vector default(present)
-                    do k = is2_viscous%beg, is2_viscous%end
-                        do j = is1_viscous%beg, is1_viscous%end
+                    do k = idwbuff(2)%beg, idwbuff(2)%end
+                        do j = idwbuff(1)%beg, idwbuff(1)%end
                             grad_z%sf(j, k, p) = &
                                 (3d0*var%sf(j, k, p) - 4d0*var%sf(j, k, p - 1) + var%sf(j, k, p - 2))/ &
                                 (z_cc(p) - z_cc(p - 2))