diff --git a/Source/radiation/MGFLD.cpp b/Source/radiation/MGFLD.cpp
index 691e0df8da..d20ea28715 100644
--- a/Source/radiation/MGFLD.cpp
+++ b/Source/radiation/MGFLD.cpp
@@ -1097,97 +1097,71 @@ void Radiation::compute_limiter(int level, const BoxArray& grids,
             int reg_klo = bx.loVect3d()[2];
             int reg_khi = bx.hiVect3d()[2];
 
-            if (filter_lambda_T == 1) {
+            if (filter_lambda_T > 0) {
+                // initialize
 
-                for (int k = lam_klo; k <= lam_khi; ++k) {
-                    for (int j = lam_jlo; j <= lam_jhi; ++j) {
-#if AMREX_SPACEDIM >= 2
-                        if (Er(reg_ilo,j,k,g) == -1.e0_rt) {
-                            for (int i = reg_ilo; i <= reg_ihi; ++i) {
-                                lamfil(i,j,k) = -1.e-50_rt;
-                            }
-                            continue;
-                        }
-#endif
+                amrex::Loop(bx, [=] (int i, int j, int k) noexcept
+                {
+                    lamfil(i,j,k) = -1.e-50_rt;
+                });
+            }
 
-                        for (int i = reg_ilo; i <= reg_ihi; ++i) {
-                            lamfil(i,j,k) = ff1(0) * lam(i,j,k,g) + ff1(1) * (lam(i-1,j,k,g) + lam(i+1,j,k,g));
-                            lamfil(i,j,k) = std::min(1.e0_rt/3.e0_rt, std::max(1.e-25_rt, lamfil(i,j,k)));
-                        }
+            if (filter_lambda_T == 1) {
 
-                        if (Er(reg_ilo-1,j,k,g) == -1.e0_rt) {
-                            int i = reg_ilo;
-                            lamfil(i,j,k) = ff1b(0) * lam(i,j,k,g) + ff1b(1) * lam(i+1,j,k,g);
-                            lamfil(i,j,k) = std::min(1.e0_rt/3.e0_rt, std::max(1.e-25_rt, lamfil(i,j,k)));
-                        }
+                // filter in the x direction
 
-                        if (Er(reg_ihi+1,j,k,g) == -1.e0_rt) {
-                            int i = reg_ihi;
-                            lamfil(i,j,k) = ff1b(1) * lam(i-1,j,k,g) + ff1b(0) * lam(i,j,k,g);
-                            lamfil(i,j,k) = std::min(1.e0_rt/3.e0_rt, std::max(1.e-25_rt, lamfil(i,j,k)));
-                        }
+                amrex::Loop(bx, [=] (int i, int j, int k) noexcept
+                {
+                    if (Er(i-1,j,k,g) == -1.e0_rt) {
+                        lamfil(i,j,k) = ff1b(0) * lam(i,j,k,g) + ff1b(1) * lam(i+1,j,k,g);
                     }
-                }
-
-#if AMREX_SPACEDIM >= 2
-                for (int k = lam_klo; k <= lam_khi; ++k) {
-                    if (Er(reg_ilo,reg_jlo,k,g) == -1.e0_rt) {
-                        continue;
+                    else if (Er(i+1,j,k,g) == -1.e0_rt) {
+                        lamfil(i,j,k) = ff1b(1) * lam(i-1,j,k,g) + ff1b(0) * lam(i,j,k,g);
                     }
-
-                    for (int j = reg_jlo; j <= reg_jhi; ++j) {
-                        for (int i = reg_ilo; i <= reg_ihi; ++i) {
-                            lam(i,j,k,g) = ff1(0) * lamfil(i,j,k) + ff1(1) * (lamfil(i,j-1,k) + lamfil(i,j+1,k));
-                        }
+                    else {
+                        lamfil(i,j,k) = ff1(0) * lam(i,j,k,g) + ff1(1) * (lam(i-1,j,k,g) + lam(i+1,j,k,g));
                     }
+                    lamfil(i,j,k) = std::min(1.e0_rt/3.e0_rt, std::max(1.e-25_rt, lamfil(i,j,k)));
+                });
 
-                    if (Er(reg_ilo,reg_jlo-1,k,g) == -1.e0_rt) {
-                        int j = reg_jlo;
-                        for (int i = reg_ilo; i <= reg_ihi; ++i) {
-                            lam(i,j,k,g) = ff1b(0) * lamfil(i,j,k) + ff1b(1) * lamfil(i,j+1,k);
-                        }
-                    }
+                // filter in the y direction
 
-                    if (Er(reg_ilo,reg_jhi+1,k,g) == -1.e0_rt) {
-                        int j = reg_jhi;
-                        for (int i = reg_ilo; i <= reg_ihi; ++i) {
-                            lam(i,j,k,g) = ff1b(1) * lamfil(i,j-1,k) + ff1b(0) * lamfil(i,j,k);
-                        }
+#if AMREX_SPACEDIM >= 2
+                amrex::Loop(bx, [=] (int i, int j, int k) noexcept
+                {
+                    if (Er(i,j-1,k,g) == -1.e0_rt) {
+                        lamfil(i,j,k,g) = ff1b(0) * lamfil(i,j,k) + ff1b(1) * lamfil(i,j+1,k);
                     }
-                }
+                    else if (Er(i,j+1,k,g) == -1.e0_rt) {
+                        lamfil(i,j,k,g) = ff1b(1) * lamfil(i,j-1,k) + ff1b(0) * lamfil(i,j,k);
+                    }
+                    else {
+                        lamfil(i,j,k,g) = ff1(0) * lamfil(i,j,k) + ff1(1) * (lamfil(i,j-1,k) + lamfil(i,j+1,k));
+                    }
+                    lamfil(i,j,k) = std::min(1.e0_rt/3.e0_rt, std::max(1.e-25_rt, lamfil(i,j,k)));
+                });
 #endif
 
+                // filter in the z direction
+
 #if AMREX_SPACEDIM == 3
-                for (int k = reg_klo; k <= reg_khi; ++k) {
-                    for (int j = reg_jlo; j <= reg_jhi; ++j) {
-                        for (int i = reg_ilo; i <= reg_ihi; ++i) {
-                            lamfil(i,j,k) = ff1(0) * lam(i,j,k,g) + ff1(1) * (lam(i,j,k-1,g) + lam(i,j,k+1,g));
-                            lamfil(i,j,k) = std::min(1.e0_rt/3.e0_rt, std::max(1.e-25_rt, lamfil(i,j,k)));
-                        }
+                amrex::Loop(bx, [=] (int i, int j, int k) noexcept
+                {
+                    if (Er(i,j,k-1,g) == -1.e0_rt) {
+                        lamfil(i,j,k) = ff1b(0) * lam(i,j,k,g) + ff1b(1) * lam(i,j,k+1,g);
                     }
-                }
-
-                if (Er(reg_ilo,reg_jlo,reg_klo-1,g) == -1.e0_rt) {
-                    int k = reg_klo;
-                    for (int j = reg_jlo; j <= reg_jhi; ++j) {
-                        for (int i = reg_ilo; i <= reg_ihi; ++i) {
-                            lamfil(i,j,k) = ff1b(0) * lam(i,j,k,g) + ff1b(1) * lam(i,j,k+1,g);
-                            lamfil(i,j,k) = std::min(1.e0_rt/3.e0_rt, std::max(1.e-25_rt, lamfil(i,j,k)));
-                        }
+                    else if (Er(i,j,k+1,g) == -1.e0_rt) {
+                        lamfil(i,j,k) = ff1b(1) * lam(i,j,k-1,g) + ff1b(0) * lam(i,j,k,g);
                     }
-                }
-
-                if (Er(reg_ilo,reg_jlo,reg_khi+1,g) == -1.e0_rt) {
-                    int k = reg_khi;
-                    for (int j = reg_jlo; j <= reg_jhi; ++j) {
-                        for (int i = reg_ilo; i <= reg_ihi; ++i) {
-                            lamfil(i,j,k) = ff1b(1) * lam(i,j,k-1,g) + ff1b(0) * lam(i,j,k,g);
-                            lamfil(i,j,k) = std::min(1.e0_rt/3.e0_rt, std::max(1.e-25_rt, lamfil(i,j,k)));
-                        }
+                    else {
+                        lamfil(i,j,k) = ff1(0) * lam(i,j,k,g) + ff1(1) * (lam(i,j,k-1,g) + lam(i,j,k+1,g));
                     }
-                }
+                    lamfil(i,j,k) = std::min(1.e0_rt/3.e0_rt, std::max(1.e-25_rt, lamfil(i,j,k)));
+                });
 #endif
 
+                // store the final results
+
                 amrex::Loop(bx, [=] (int i, int j, int k) noexcept
                 {
                     lam(i,j,k,g) = lamfil(i,j,k);