diff --git a/Source/radiation/RAD_1D.F90 b/Source/radiation/RAD_1D.F90
index bef1b8620e..b960b07211 100644
--- a/Source/radiation/RAD_1D.F90
+++ b/Source/radiation/RAD_1D.F90
@@ -26,19 +26,4 @@ subroutine rfface(fine, &
   fine(fbox_l1) = crse(cbox_l1)
 end subroutine rfface
 
-subroutine bextrp(f, fbox_l1, fbox_h1, reg_l1, reg_h1) bind(C, name="bextrp")
-  use amrex_fort_module, only : rt => amrex_real
-  implicit none
-  integer :: fbox_l1, fbox_h1
-  integer ::  reg_l1,  reg_h1
-  real(rt)         :: f(fbox_l1:fbox_h1)
-  integer :: i
-
-  !     i direction first:
-  i = reg_l1
-  f(i-1) = 2.e0_rt * f(i) - f(i+1)
-  i = reg_h1
-  f(i+1) = 2.e0_rt * f(i) - f(i-1)
-end subroutine bextrp
-
 end module rad_module
diff --git a/Source/radiation/RAD_2D.F90 b/Source/radiation/RAD_2D.F90
index f32f5a8ea5..936db98070 100644
--- a/Source/radiation/RAD_2D.F90
+++ b/Source/radiation/RAD_2D.F90
@@ -35,33 +35,4 @@ subroutine rfface(fine, &
   endif
 end subroutine rfface
 
-
-subroutine bextrp(f, fbox_l1, fbox_l2, fbox_h1, fbox_h2, &
-                  reg_l1, reg_l2, reg_h1, reg_h2) bind(C, name="bextrp")
-
-  use amrex_fort_module, only : rt => amrex_real
-  integer :: fbox_l1, fbox_l2, fbox_h1, fbox_h2
-  integer ::  reg_l1,  reg_l2,  reg_h1,  reg_h2
-  real(rt)         :: f(fbox_l1:fbox_h1,fbox_l2:fbox_h2)
-  integer :: i, j
-
-  !     i direction first:
-  do j = reg_l2, reg_h2
-     i = reg_l1
-     f(i-1,j) = 2.e0_rt * f(i,j) - f(i+1,j)
-     i = reg_h1
-     f(i+1,j) = 2.e0_rt * f(i,j) - f(i-1,j)
-  enddo
-
-  !     j direction second, including corners:
-  do i = reg_l1 - 1, reg_h1 + 1
-     j = reg_l2
-     f(i,j-1) = 2.e0_rt * f(i,j) - f(i,j+1)
-     j = reg_h2
-     f(i,j+1) = 2.e0_rt * f(i,j) - f(i,j-1)
-  enddo
-
-  !  corner results are the same whichever direction we extrapolate first
-end subroutine bextrp
-
 end module rad_module
diff --git a/Source/radiation/RAD_3D.F90 b/Source/radiation/RAD_3D.F90
index ebdbdee21a..5b0f3dd71a 100644
--- a/Source/radiation/RAD_3D.F90
+++ b/Source/radiation/RAD_3D.F90
@@ -52,49 +52,5 @@ subroutine rfface(fine, &
   endif
 end subroutine rfface
 
-
-subroutine bextrp( &
-                  f, fbox_l1, fbox_l2, fbox_l3, fbox_h1, fbox_h2, fbox_h3, &
-                  reg_l1, reg_l2, reg_l3, reg_h1, reg_h2, reg_h3) bind(C, name="bextrp")
-
-  use amrex_fort_module, only : rt => amrex_real
-  integer :: fbox_l1, fbox_l2, fbox_l3, fbox_h1, fbox_h2, fbox_h3
-  integer ::  reg_l1,  reg_l2,  reg_l3,  reg_h1,  reg_h2,  reg_h3
-  real(rt)         :: f(fbox_l1:fbox_h1,fbox_l2:fbox_h2,fbox_l3:fbox_h3)
-  integer :: i, j, k
-
-  !     i direction first:
-  do k = reg_l3, reg_h3
-     do j = reg_l2, reg_h2
-        i = reg_l1
-        f(i-1,j,k) = 2.e0_rt * f(i,j,k) - f(i+1,j,k)
-        i = reg_h1
-        f(i+1,j,k) = 2.e0_rt * f(i,j,k) - f(i-1,j,k)
-     enddo
-  enddo
-
-  !     j direction second, including edges:
-  do k = reg_l3, reg_h3
-     do i = reg_l1 - 1, reg_h1 + 1
-        j = reg_l2
-        f(i,j-1,k) = 2.e0_rt * f(i,j,k) - f(i,j+1,k)
-        j = reg_h2
-        f(i,j+1,k) = 2.e0_rt * f(i,j,k) - f(i,j-1,k)
-     enddo
-  enddo
-
-  !     k direction third, including corners:
-  do j = reg_l2 - 1, reg_h2 + 1
-     do i = reg_l1 - 1, reg_h1 + 1
-        k = reg_l3
-        f(i,j,k-1) = 2.e0_rt * f(i,j,k) - f(i,j,k+1)
-        k = reg_h3
-        f(i,j,k+1) = 2.e0_rt * f(i,j,k) - f(i,j,k-1)
-     enddo
-  enddo
-
-  !   corner results are the same whichever direction we extrapolate first
-end subroutine bextrp
-
 end module rad_module
 
diff --git a/Source/radiation/RAD_F.H b/Source/radiation/RAD_F.H
index f77512daa0..2306d83fa6 100644
--- a/Source/radiation/RAD_F.H
+++ b/Source/radiation/RAD_F.H
@@ -153,9 +153,6 @@ extern "C" {
                       int* tfab, ARLIM_P(dlo), ARLIM_P(dhi),
                       const amrex::Real* dx, const amrex::Real* xlo, const amrex::Real& time);
 
-  void bextrp(BL_FORT_FAB_ARG(f), 
-              ARLIM_P(reglo), ARLIM_P(reghi));
-
 #ifdef __cplusplus
 }
 #endif
diff --git a/Source/radiation/Radiation.cpp b/Source/radiation/Radiation.cpp
index 975f07202b..ace2a6969b 100644
--- a/Source/radiation/Radiation.cpp
+++ b/Source/radiation/Radiation.cpp
@@ -1211,21 +1211,54 @@ void Radiation::state_update(MultiFab& state, MultiFab& frhoes)
 
 void Radiation::extrapolateBorders(MultiFab& f, int indx)
 {
-  BL_PROFILE("Radiation::extrapolateBorders");
+    BL_PROFILE("Radiation::extrapolateBorders");
 
-  BL_ASSERT(f.nGrow() >= 1);
+    BL_ASSERT(f.nGrow() >= 1);
 
 #ifdef _OPENMP
 #pragma omp parallel
 #endif
-  for(MFIter mfi(f); mfi.isValid(); ++mfi) {
-    int i = mfi.index();
+    for (MFIter mfi(f); mfi.isValid(); ++mfi) {
+        // Note no tiling in the current implementation.
+        const Box& bx = mfi.validbox();
+        const Box& grownbx = amrex::grow(bx, 1);
 
-    const Box& reg  = f.box(i);
+        Array4<Real> const f_arr = f[mfi].array(indx);
 
-    bextrp(BL_TO_FORTRAN_N(f[mfi],indx),
-           ARLIM(reg.loVect()), ARLIM(reg.hiVect()));
-  }
+        amrex::LoopOnCpu(grownbx, [=] (int i, int j, int k) noexcept
+        {
+            // Note that the results on the corners will be the same
+            // regardless of which order we do the loop in.
+
+            if (i == bx.smallEnd(0)) {
+                f_arr(i-1,j,k) = 2.0_rt * f_arr(i,j,k) - f_arr(i+1,j,k);
+            }
+
+            if (i == bx.bigEnd(0)) {
+                f_arr(i+1,j,k) = 2.0_rt * f_arr(i,j,k) - f_arr(i-1,j,k);
+            }
+
+#if AMREX_SPACEDIM >= 2
+            if (j == bx.smallEnd(1)) {
+                f_arr(i,j-1,k) = 2.0_rt * f_arr(i,j,k) - f_arr(i,j+1,k);
+            }
+
+            if (j == bx.bigEnd(1)) {
+                f_arr(i,j+1,k) = 2.0_rt * f_arr(i,j,k) - f_arr(i,j-1,k);
+            }
+#endif
+
+#if AMREX_SPACEDIM == 3
+            if (k == bx.smallEnd(2)) {
+                f_arr(i,j,k-1) = 2.0_rt * f_arr(i,j,k) - f_arr(i,j,k+1);
+            }
+
+            if (k == bx.bigEnd(2)) {
+                f_arr(i,j,k+1) = 2.0_rt * f_arr(i,j,k) - f_arr(i,j,k-1);
+            }
+#endif
+        });
+    }
 }