diff --git a/Source/ERF.cpp b/Source/ERF.cpp
index 1275e3d28..72f325e88 100644
--- a/Source/ERF.cpp
+++ b/Source/ERF.cpp
@@ -1525,33 +1525,33 @@ ERF::AverageDownTo (int crse_lev) // NOLINT
 void
 ERF::define_grids_to_evolve (int lev) // NOLINT
 {
-   int width = wrfbdy_width - 1;
-   Box domain(geom[lev].Domain());
-   if (lev == 0 && ( init_type == "real" || init_type == "metgrid" ) )
-   {
-      Box shrunk_domain(domain);
-      shrunk_domain.grow(0,-width);
-      shrunk_domain.grow(1,-width);
-      grids_to_evolve[lev] = amrex::intersect(grids[lev],shrunk_domain);
-   } else if (lev == 1) {
-      Box shrunk_domain(boxes_at_level[lev][0]);
-      shrunk_domain.grow(0,-width);
-      shrunk_domain.grow(1,-width);
-      grids_to_evolve[lev] = amrex::intersect(grids[lev],shrunk_domain);
+    int width = wrfbdy_set_width;
+    Box domain(geom[lev].Domain());
+    if (lev == 0 && ( init_type == "real" || init_type == "metgrid" ) )
+    {
+        Box shrunk_domain(domain);
+        shrunk_domain.grow(0,-width);
+        shrunk_domain.grow(1,-width);
+        grids_to_evolve[lev] = amrex::intersect(grids[lev],shrunk_domain);
+    } else if (lev == 1) {
+        Box shrunk_domain(boxes_at_level[lev][0]);
+        shrunk_domain.grow(0,-width);
+        shrunk_domain.grow(1,-width);
+        grids_to_evolve[lev] = amrex::intersect(grids[lev],shrunk_domain);
 #if 0
-      if (num_boxes_at_level[lev] > 1) {
-          for (int i = 1; i < num_boxes_at_level[lev]; i++) {
-              Box shrunk_domain(boxes_at_level[lev][i]);
-              shrunk_domain.grow(0,-width);
-              shrunk_domain.grow(1,-width);
-              grids_to_evolve[lev] = amrex::intersect(grids_to_evolve[lev],shrunk_domain);
-          }
-      }
+        if (num_boxes_at_level[lev] > 1) {
+            for (int i = 1; i < num_boxes_at_level[lev]; i++) {
+                Box shrunk_domain(boxes_at_level[lev][i]);
+                shrunk_domain.grow(0,-width);
+                shrunk_domain.grow(1,-width);
+                grids_to_evolve[lev] = amrex::intersect(grids_to_evolve[lev],shrunk_domain);
+            }
+        }
 #endif
-   } else {
-      // Just copy grids...
-      grids_to_evolve[lev] = grids[lev];
-   }
+    } else {
+        // Just copy grids...
+        grids_to_evolve[lev] = grids[lev];
+    }
 }
 
 #ifdef ERF_USE_MULTIBLOCK
diff --git a/Source/TimeIntegration/TI_slow_rhs_fun.H b/Source/TimeIntegration/TI_slow_rhs_fun.H
index 355f5c147..58666f797 100644
--- a/Source/TimeIntegration/TI_slow_rhs_fun.H
+++ b/Source/TimeIntegration/TI_slow_rhs_fun.H
@@ -194,7 +194,7 @@
 
 #ifdef ERF_USE_NETCDF
         // Populate RHS for relaxation zones
-        if (init_type == "real" && level == 0)
+        if (init_type=="real" && level==0 && nrk==0)
             compute_interior_ghost_RHS(bdy_time_interval, old_stage_time, slow_dt,
                                        wrfbdy_width-1, wrfbdy_set_width, fine_geom,
                                        S_rhs, S_data, bdy_data_xlo, bdy_data_xhi,
diff --git a/Source/Utils/InteriorGhostCells.cpp b/Source/Utils/InteriorGhostCells.cpp
index 9bd72fdb3..e5ca25986 100644
--- a/Source/Utils/InteriorGhostCells.cpp
+++ b/Source/Utils/InteriorGhostCells.cpp
@@ -420,7 +420,7 @@ compute_interior_ghost_RHS(const Real& bdy_time_interval,
             // 4 halo boxes w/ no ghost cells for CC vars
             Box domain = geom.Domain();
             Box bx_xlo, bx_xhi, bx_ylo, bx_yhi;
-            compute_interior_ghost_bxs_xy(tbx, domain, width, 0,
+            compute_interior_ghost_bxs_xy(tbx, domain, set_width, 0,
                                           bx_xlo, bx_xhi,
                                           bx_ylo, bx_yhi);
             S_rhs[IntVar::cons][mfi].template setVal<RunOn::Device>(0.,bx_xlo,Rho_comp,NVAR);
@@ -437,7 +437,7 @@ compute_interior_ghost_RHS(const Real& bdy_time_interval,
             Box domain = geom.Domain();
             domain.convert(S_rhs[IntVar::zmom].boxArray().ixType());
             Box bx_xlo, bx_xhi, bx_ylo, bx_yhi;
-            compute_interior_ghost_bxs_xy(tbx, domain, width, 0,
+            compute_interior_ghost_bxs_xy(tbx, domain, set_width, 0,
                                           bx_xlo, bx_xhi,
                                           bx_ylo, bx_yhi);
             S_rhs[IntVar::zmom][mfi].template setVal<RunOn::Device>(0.,bx_xlo);
@@ -457,6 +457,7 @@ compute_interior_ghost_RHS(const Real& bdy_time_interval,
         Box domain = geom.Domain();
         domain.convert(S_data[var_idx].boxArray().ixType());
         const auto& dom_hi = ubound(domain);
+        const auto& dom_lo = lbound(domain);
 
 #ifdef _OPENMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
@@ -517,27 +518,32 @@ compute_interior_ghost_RHS(const Real& bdy_time_interval,
             [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
                 // Corners with x boxes
-                int j_lo = std::min(j,width-1);
+                int j_lo = std::min(j-dom_lo.y,width-1);
                 int j_hi = std::min(dom_hi.y-j,width-1);
                 int jj   = std::min(j_lo,j_hi);
-                int n    = std::min(i,jj) + 1;
+                int n    = std::min(i-dom_lo.x,jj) + 1;
                 if (n <= Spec_z) {
                     rhs_arr(i,j,k,icomp) = 0.0;
                 } else {
                     Real Factor   = (num - Real(n))/denom;
-                    Real delta    = arr_xlo(i  ,j  ,k) - data_arr(i  ,j  ,k,icomp);
-                    Real delta_xp = arr_xlo(i+1,j  ,k) - data_arr(i+1,j  ,k,icomp);
-                    Real delta_xm = arr_xlo(i-1,j  ,k) - data_arr(i-1,j  ,k,icomp);
-                    Real delta_yp = arr_xlo(i  ,j+1,k) - data_arr(i  ,j+1,k,icomp);
-                    Real delta_ym = arr_xlo(i  ,j-1,k) - data_arr(i  ,j-1,k,icomp);
+                    Real d        = data_arr(i  ,j  ,k,icomp) + delta_t*rhs_arr(i  , j  , k,icomp);
+                    Real d_ip1    = data_arr(i+1,j  ,k,icomp) + delta_t*rhs_arr(i+1, j  , k,icomp);
+                    Real d_im1    = data_arr(i-1,j  ,k,icomp) + delta_t*rhs_arr(i-1, j  , k,icomp);
+                    Real d_jp1    = data_arr(i  ,j+1,k,icomp) + delta_t*rhs_arr(i  , j+1, k,icomp);
+                    Real d_jm1    = data_arr(i  ,j-1,k,icomp) + delta_t*rhs_arr(i  , j-1, k,icomp);
+                    Real delta    = arr_xlo(i  ,j  ,k) - d;
+                    Real delta_xp = arr_xlo(i+1,j  ,k) - d_ip1;
+                    Real delta_xm = arr_xlo(i-1,j  ,k) - d_im1;
+                    Real delta_yp = arr_xlo(i  ,j+1,k) - d_jp1;
+                    Real delta_ym = arr_xlo(i  ,j-1,k) - d_jm1;
                     Real Laplacian = delta_xp + delta_xm + delta_yp + delta_ym - 4.0*delta;
-                    rhs_arr(i,j,k,icomp) = (F1*delta - F2*Laplacian) * Factor;
+                    rhs_arr(i,j,k,icomp) += (F1*delta - F2*Laplacian) * Factor;
                 }
             },
             [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
                 // Corners with x boxes
-                int j_lo = std::min(j,width-1);
+                int j_lo = std::min(j-dom_lo.y,width-1);
                 int j_hi = std::min(dom_hi.y-j,width-1);
                 int jj   = std::min(j_lo,j_hi);
                 int n    = std::min(dom_hi.x-i,jj) + 1;
@@ -545,13 +551,18 @@ compute_interior_ghost_RHS(const Real& bdy_time_interval,
                     rhs_arr(i,j,k,icomp) = 0.0;
                 } else {
                     Real Factor   = (num - Real(n))/denom;
-                    Real delta    = arr_xhi(i  ,j  ,k) - data_arr(i  ,j  ,k,icomp);
-                    Real delta_xp = arr_xhi(i+1,j  ,k) - data_arr(i+1,j  ,k,icomp);
-                    Real delta_xm = arr_xhi(i-1,j  ,k) - data_arr(i-1,j  ,k,icomp);
-                    Real delta_yp = arr_xhi(i  ,j+1,k) - data_arr(i  ,j+1,k,icomp);
-                    Real delta_ym = arr_xhi(i  ,j-1,k) - data_arr(i  ,j-1,k,icomp);
+                    Real d        = data_arr(i  ,j  ,k,icomp) + delta_t*rhs_arr(i  , j  , k,icomp);
+                    Real d_ip1    = data_arr(i+1,j  ,k,icomp) + delta_t*rhs_arr(i+1, j  , k,icomp);
+                    Real d_im1    = data_arr(i-1,j  ,k,icomp) + delta_t*rhs_arr(i-1, j  , k,icomp);
+                    Real d_jp1    = data_arr(i  ,j+1,k,icomp) + delta_t*rhs_arr(i  , j+1, k,icomp);
+                    Real d_jm1    = data_arr(i  ,j-1,k,icomp) + delta_t*rhs_arr(i  , j-1, k,icomp);
+                    Real delta    = arr_xlo(i  ,j  ,k) - d;
+                    Real delta_xp = arr_xlo(i+1,j  ,k) - d_ip1;
+                    Real delta_xm = arr_xlo(i-1,j  ,k) - d_im1;
+                    Real delta_yp = arr_xlo(i  ,j+1,k) - d_jp1;
+                    Real delta_ym = arr_xlo(i  ,j-1,k) - d_jm1;
                     Real Laplacian = delta_xp + delta_xm + delta_yp + delta_ym - 4.0*delta;
-                    rhs_arr(i,j,k,icomp) = (F1*delta - F2*Laplacian) * Factor;
+                    rhs_arr(i,j,k,icomp) += (F1*delta - F2*Laplacian) * Factor;
                 }
             });
 
@@ -559,18 +570,23 @@ compute_interior_ghost_RHS(const Real& bdy_time_interval,
             [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
                 // No corners for y boxes
-                int n = j + 1;
+                int n = j -dom_lo.y + 1;
                 if (n <= Spec_z) {
                     rhs_arr(i,j,k,icomp) = 0.0;
                 } else {
                     Real Factor   = (num - Real(n))/denom;
-                    Real delta    = arr_ylo(i  ,j  ,k) - data_arr(i  ,j  ,k,icomp);
-                    Real delta_xp = arr_ylo(i+1,j  ,k) - data_arr(i+1,j  ,k,icomp);
-                    Real delta_xm = arr_ylo(i-1,j  ,k) - data_arr(i-1,j  ,k,icomp);
-                    Real delta_yp = arr_ylo(i  ,j+1,k) - data_arr(i  ,j+1,k,icomp);
-                    Real delta_ym = arr_ylo(i  ,j-1,k) - data_arr(i  ,j-1,k,icomp);
+                    Real d        = data_arr(i  ,j  ,k,icomp) + delta_t*rhs_arr(i  , j  , k,icomp);
+                    Real d_ip1    = data_arr(i+1,j  ,k,icomp) + delta_t*rhs_arr(i+1, j  , k,icomp);
+                    Real d_im1    = data_arr(i-1,j  ,k,icomp) + delta_t*rhs_arr(i-1, j  , k,icomp);
+                    Real d_jp1    = data_arr(i  ,j+1,k,icomp) + delta_t*rhs_arr(i  , j+1, k,icomp);
+                    Real d_jm1    = data_arr(i  ,j-1,k,icomp) + delta_t*rhs_arr(i  , j-1, k,icomp);
+                    Real delta    = arr_xlo(i  ,j  ,k) - d;
+                    Real delta_xp = arr_xlo(i+1,j  ,k) - d_ip1;
+                    Real delta_xm = arr_xlo(i-1,j  ,k) - d_im1;
+                    Real delta_yp = arr_xlo(i  ,j+1,k) - d_jp1;
+                    Real delta_ym = arr_xlo(i  ,j-1,k) - d_jm1;
                     Real Laplacian = delta_xp + delta_xm + delta_yp + delta_ym - 4.0*delta;
-                    rhs_arr(i,j,k,icomp) = (F1*delta - F2*Laplacian) * Factor;
+                    rhs_arr(i,j,k,icomp) += (F1*delta - F2*Laplacian) * Factor;
                 }
             },
             [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
@@ -581,13 +597,18 @@ compute_interior_ghost_RHS(const Real& bdy_time_interval,
                     rhs_arr(i,j,k,icomp) = 0.0;
                 } else {
                     Real Factor   = (num - Real(n))/denom;
-                    Real delta    = arr_yhi(i  ,j  ,k) - data_arr(i  ,j  ,k,icomp);
-                    Real delta_xp = arr_yhi(i+1,j  ,k) - data_arr(i+1,j  ,k,icomp);
-                    Real delta_xm = arr_yhi(i-1,j  ,k) - data_arr(i-1,j  ,k,icomp);
-                    Real delta_yp = arr_yhi(i  ,j+1,k) - data_arr(i  ,j+1,k,icomp);
-                    Real delta_ym = arr_yhi(i  ,j-1,k) - data_arr(i  ,j-1,k,icomp);
+                    Real d        = data_arr(i  ,j  ,k,icomp) + delta_t*rhs_arr(i  , j  , k,icomp);
+                    Real d_ip1    = data_arr(i+1,j  ,k,icomp) + delta_t*rhs_arr(i+1, j  , k,icomp);
+                    Real d_im1    = data_arr(i-1,j  ,k,icomp) + delta_t*rhs_arr(i-1, j  , k,icomp);
+                    Real d_jp1    = data_arr(i  ,j+1,k,icomp) + delta_t*rhs_arr(i  , j+1, k,icomp);
+                    Real d_jm1    = data_arr(i  ,j-1,k,icomp) + delta_t*rhs_arr(i  , j-1, k,icomp);
+                    Real delta    = arr_xlo(i  ,j  ,k) - d;
+                    Real delta_xp = arr_xlo(i+1,j  ,k) - d_ip1;
+                    Real delta_xm = arr_xlo(i-1,j  ,k) - d_im1;
+                    Real delta_yp = arr_xlo(i  ,j+1,k) - d_jp1;
+                    Real delta_ym = arr_xlo(i  ,j-1,k) - d_jm1;
                     Real Laplacian = delta_xp + delta_xm + delta_yp + delta_ym - 4.0*delta;
-                    rhs_arr(i,j,k,icomp) = (F1*delta - F2*Laplacian) * Factor;
+                    rhs_arr(i,j,k,icomp) += (F1*delta - F2*Laplacian) * Factor;
                 }
             });
         } // mfi