diff --git a/amr-wind/diffusion/incflo_diffusion.cpp b/amr-wind/diffusion/incflo_diffusion.cpp
index da951dfa4c..ada217cc52 100644
--- a/amr-wind/diffusion/incflo_diffusion.cpp
+++ b/amr-wind/diffusion/incflo_diffusion.cpp
@@ -223,6 +223,7 @@ void viscosity_to_uniform_space(
         repo.get_mesh_mapping_det_j(amr_wind::FieldLoc::ZFACE);
 
     // beta accounted for mesh mapping (x-face) = J/fac^2 * mu
+    // needs openmp pragma?
     for (amrex::MFIter mfi(b[0]); mfi.isValid(); ++mfi) {
         amrex::Array4<amrex::Real> const& mu = b[0].array(mfi);
         amrex::Array4<amrex::Real const> const& fac =
@@ -237,6 +238,7 @@ void viscosity_to_uniform_space(
             });
     }
     // beta accounted for mesh mapping (y-face) = J/fac^2 * mu
+    // needs openmp pragma?
     for (amrex::MFIter mfi(b[1]); mfi.isValid(); ++mfi) {
         amrex::Array4<amrex::Real> const& mu = b[1].array(mfi);
         amrex::Array4<amrex::Real const> const& fac =
@@ -251,6 +253,7 @@ void viscosity_to_uniform_space(
             });
     }
     // beta accounted for mesh mapping (z-face) = J/fac^2 * mu
+    // needs openmp pragma?
     for (amrex::MFIter mfi(b[2]); mfi.isValid(); ++mfi) {
         amrex::Array4<amrex::Real> const& mu = b[2].array(mfi);
         amrex::Array4<amrex::Real const> const& fac =
diff --git a/amr-wind/equation_systems/CompRHSOps.H b/amr-wind/equation_systems/CompRHSOps.H
index b0a83415c6..80fd1face0 100644
--- a/amr-wind/equation_systems/CompRHSOps.H
+++ b/amr-wind/equation_systems/CompRHSOps.H
@@ -79,6 +79,7 @@ struct ComputeRHSOp
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
+            // this one
             for (amrex::MFIter mfi(field(lev)); mfi.isValid(); ++mfi) {
                 const auto& bx = mfi.tilebox();
                 auto fld = field(lev).array(mfi);
@@ -199,6 +200,7 @@ struct ComputeRHSOp
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
+            // this one
             for (amrex::MFIter mfi(field(lev)); mfi.isValid(); ++mfi) {
                 const auto& bx = mfi.tilebox();
                 auto fld = field(lev).array(mfi);
diff --git a/amr-wind/equation_systems/DiffusionOps.cpp b/amr-wind/equation_systems/DiffusionOps.cpp
index 73a9cb5c1e..e8091fb0c3 100644
--- a/amr-wind/equation_systems/DiffusionOps.cpp
+++ b/amr-wind/equation_systems/DiffusionOps.cpp
@@ -128,6 +128,7 @@ void DiffSolverIface<LinOp>::linsys_solve_impl()
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
+        // this one
         for (amrex::MFIter mfi(rhs, amrex::TilingIfNotGPU()); mfi.isValid();
              ++mfi) {
             const auto& bx = mfi.tilebox();
diff --git a/amr-wind/equation_systems/PDEOps.H b/amr-wind/equation_systems/PDEOps.H
index 6dcdd7e955..2e64d0f34d 100644
--- a/amr-wind/equation_systems/PDEOps.H
+++ b/amr-wind/equation_systems/PDEOps.H
@@ -100,6 +100,7 @@ struct SrcTermOpBase
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
+            // this one
             for (amrex::MFIter mfi(src_term, amrex::TilingIfNotGPU());
                  mfi.isValid(); ++mfi) {
                 const auto& bx = mfi.tilebox();
diff --git a/amr-wind/equation_systems/density/density_ops.H b/amr-wind/equation_systems/density/density_ops.H
index d1e15eb7fd..fdd1b62bf9 100644
--- a/amr-wind/equation_systems/density/density_ops.H
+++ b/amr-wind/equation_systems/density/density_ops.H
@@ -30,6 +30,7 @@ struct ComputeRHSOp<Density, Scheme>
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
+            // this one
             for (amrex::MFIter mfi(field(lev)); mfi.isValid(); ++mfi) {
                 const auto& bx = mfi.tilebox();
                 auto rho = field(lev).array(mfi);
@@ -62,6 +63,7 @@ struct ComputeRHSOp<Density, Scheme>
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
+            // this one
             for (amrex::MFIter mfi(field(lev)); mfi.isValid(); ++mfi) {
                 const auto& bx = mfi.tilebox();
                 auto rho = field(lev).array(mfi);
diff --git a/amr-wind/equation_systems/icns/icns_advection.cpp b/amr-wind/equation_systems/icns/icns_advection.cpp
index 92415663f2..ea4cfb1e7a 100644
--- a/amr-wind/equation_systems/icns/icns_advection.cpp
+++ b/amr-wind/equation_systems/icns/icns_advection.cpp
@@ -337,6 +337,8 @@ void MacProjOp::mac_proj_to_uniform_space(
     // scale U^mac to accommodate for mesh mapping -> U^bar = J/fac *
     // U^mac beta accounted for mesh mapping = J/fac^2 * 1/rho construct
     // rho and mesh map u_mac on x-face
+    // this one
+    // needs openmp pragma?
     for (amrex::MFIter mfi(*(rho_face[0])); mfi.isValid(); ++mfi) {
         amrex::Array4<amrex::Real> const& u = u_mac(lev).array(mfi);
         amrex::Array4<amrex::Real> const& rho = rho_face[0]->array(mfi);
@@ -353,6 +355,8 @@ void MacProjOp::mac_proj_to_uniform_space(
             });
     }
     // construct rho on y-face
+    // this one
+    // needs openmp pragma?
     for (amrex::MFIter mfi(*(rho_face[1])); mfi.isValid(); ++mfi) {
         amrex::Array4<amrex::Real> const& v = v_mac(lev).array(mfi);
         amrex::Array4<amrex::Real> const& rho = rho_face[1]->array(mfi);
@@ -369,6 +373,8 @@ void MacProjOp::mac_proj_to_uniform_space(
             });
     }
     // construct rho on z-face
+    // this one
+    // needs openmp pragma?
     for (amrex::MFIter mfi(*(rho_face[2])); mfi.isValid(); ++mfi) {
         amrex::Array4<amrex::Real> const& w = w_mac(lev).array(mfi);
         amrex::Array4<amrex::Real> const& rho = rho_face[2]->array(mfi);
diff --git a/amr-wind/equation_systems/icns/icns_diffusion.H b/amr-wind/equation_systems/icns/icns_diffusion.H
index 2792e27c1b..0234df82df 100644
--- a/amr-wind/equation_systems/icns/icns_diffusion.H
+++ b/amr-wind/equation_systems/icns/icns_diffusion.H
@@ -169,6 +169,7 @@ public:
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
+                // this one
                 for (amrex::MFIter mfi(divtau(lev), amrex::TilingIfNotGPU());
                      mfi.isValid(); ++mfi) {
                     const auto& bx = mfi.tilebox();
@@ -240,6 +241,7 @@ public:
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
+            // this one
             for (amrex::MFIter mfi(rhs, amrex::TilingIfNotGPU()); mfi.isValid();
                  ++mfi) {
                 const auto& bx = mfi.tilebox();
@@ -401,6 +403,7 @@ public:
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
+                // this one
                 for (amrex::MFIter mfi(divtau(lev), amrex::TilingIfNotGPU());
                      mfi.isValid(); ++mfi) {
                     const auto& bx = mfi.tilebox();
@@ -483,6 +486,7 @@ public:
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
+            // this one
             for (amrex::MFIter mfi(rhs, amrex::TilingIfNotGPU()); mfi.isValid();
                  ++mfi) {
                 const auto& bx = mfi.tilebox();
diff --git a/amr-wind/equation_systems/icns/icns_ops.H b/amr-wind/equation_systems/icns/icns_ops.H
index 7acb631144..44ff96e630 100644
--- a/amr-wind/equation_systems/icns/icns_ops.H
+++ b/amr-wind/equation_systems/icns/icns_ops.H
@@ -93,6 +93,7 @@ struct SrcTermOp<ICNS> : SrcTermOpBase<ICNS>
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
+            // this one
             for (amrex::MFIter mfi(src_term, amrex::TilingIfNotGPU());
                  mfi.isValid(); ++mfi) {
                 const auto& bx = mfi.tilebox();
diff --git a/amr-wind/equation_systems/levelset/levelset_ops.H b/amr-wind/equation_systems/levelset/levelset_ops.H
index 3b29b9299e..3742157acf 100644
--- a/amr-wind/equation_systems/levelset/levelset_ops.H
+++ b/amr-wind/equation_systems/levelset/levelset_ops.H
@@ -69,6 +69,7 @@ struct ComputeRHSOp<Levelset, Scheme>
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
+            // this one
             for (amrex::MFIter mfi(field(lev)); mfi.isValid(); ++mfi) {
                 const auto& bx = mfi.tilebox();
                 auto phi = field(lev).array(mfi);
@@ -101,6 +102,7 @@ struct ComputeRHSOp<Levelset, Scheme>
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
+            // this one
             for (amrex::MFIter mfi(field(lev)); mfi.isValid(); ++mfi) {
                 const auto& bx = mfi.tilebox();
                 auto phi = field(lev).array(mfi);
diff --git a/amr-wind/equation_systems/tke/tke_ops.H b/amr-wind/equation_systems/tke/tke_ops.H
index c395100b1f..93f96bc8cc 100644
--- a/amr-wind/equation_systems/tke/tke_ops.H
+++ b/amr-wind/equation_systems/tke/tke_ops.H
@@ -42,6 +42,7 @@ struct PostSolveOp<TKE>
         const int nlevels = repo.num_active_levels();
         const auto clip_value = m_clip_value;
         for (int lev = 0; lev < nlevels; ++lev) {
+            // this one
             for (amrex::MFIter mfi(field(lev)); mfi.isValid(); ++mfi) {
                 const auto& bx = mfi.tilebox();
                 const auto& field_arr = field(lev).array(mfi);
diff --git a/amr-wind/equation_systems/vof/vof_hybridsolver_ops.H b/amr-wind/equation_systems/vof/vof_hybridsolver_ops.H
index f22fa4ffc3..cb0a0ffd1f 100644
--- a/amr-wind/equation_systems/vof/vof_hybridsolver_ops.H
+++ b/amr-wind/equation_systems/vof/vof_hybridsolver_ops.H
@@ -19,6 +19,7 @@ static void replace_masked_vof(
         auto& vof = f_vof(lev);
         const auto& vof_new = f_vof_new(lev);
 
+        // this one
         for (amrex::MFIter mfi(iblank); mfi.isValid(); ++mfi) {
             const auto& gbx = mfi.growntilebox();
             const amrex::Array4<const int>& native_flag =
diff --git a/amr-wind/immersed_boundary/bluff_body/bluff_body_ops.cpp b/amr-wind/immersed_boundary/bluff_body/bluff_body_ops.cpp
index 516de2f489..4e49f3438a 100644
--- a/amr-wind/immersed_boundary/bluff_body/bluff_body_ops.cpp
+++ b/amr-wind/immersed_boundary/bluff_body/bluff_body_ops.cpp
@@ -47,6 +47,7 @@ void apply_mms_vel(CFDSim& sim)
         const auto& dx = geom[lev].CellSizeArray();
         const auto& problo = geom[lev].ProbLoArray();
 
+        // this one
         for (amrex::MFIter mfi(levelset(lev)); mfi.isValid(); ++mfi) {
             const auto& bx = mfi.growntilebox();
             const auto& phi = levelset(lev).const_array(mfi);
@@ -90,6 +91,7 @@ void apply_dirichlet_vel(CFDSim& sim, const amrex::Vector<amrex::Real>& vel_bc)
         // Defining the "ghost-cell" band distance
         amrex::Real phi_b = std::cbrt(dx[0] * dx[1] * dx[2]);
 
+        // this one
         for (amrex::MFIter mfi(levelset(lev)); mfi.isValid(); ++mfi) {
             const auto& bx = mfi.tilebox();
             const auto& varr = velocity(lev).array(mfi);
diff --git a/amr-wind/immersed_boundary/bluff_body/box_ops.H b/amr-wind/immersed_boundary/bluff_body/box_ops.H
index 11567c6366..f7864616da 100644
--- a/amr-wind/immersed_boundary/bluff_body/box_ops.H
+++ b/amr-wind/immersed_boundary/bluff_body/box_ops.H
@@ -59,6 +59,7 @@ struct InitDataOp<Box>
         for (int lev = 0; lev < nlevels; ++lev) {
             const auto& problo = geom[lev].ProbLoArray();
             const auto& dx = geom[lev].CellSizeArray();
+            // this one
             for (amrex::MFIter mfi(levelset(lev)); mfi.isValid(); ++mfi) {
                 const auto& bx = mfi.growntilebox();
                 const auto& epsilon_node = mask_node(lev).array(mfi);
diff --git a/amr-wind/immersed_boundary/bluff_body/cylinder_ops.H b/amr-wind/immersed_boundary/bluff_body/cylinder_ops.H
index 2002f661b6..46440c0c33 100644
--- a/amr-wind/immersed_boundary/bluff_body/cylinder_ops.H
+++ b/amr-wind/immersed_boundary/bluff_body/cylinder_ops.H
@@ -58,6 +58,7 @@ struct InitDataOp<Cylinder>
             const auto& problo = geom[lev].ProbLoArray();
             const auto& dx = geom[lev].CellSizeArray();
 
+            // this one
             for (amrex::MFIter mfi(levelset(lev)); mfi.isValid(); ++mfi) {
                 const auto& phi = levelset(lev).array(mfi);
 
diff --git a/amr-wind/immersed_boundary/bluff_body/sphere_ops.H b/amr-wind/immersed_boundary/bluff_body/sphere_ops.H
index ec8fc94709..cda68c77ff 100644
--- a/amr-wind/immersed_boundary/bluff_body/sphere_ops.H
+++ b/amr-wind/immersed_boundary/bluff_body/sphere_ops.H
@@ -56,6 +56,7 @@ struct InitDataOp<Sphere>
             const auto& problo = geom[lev].ProbLoArray();
             const auto& dx = geom[lev].CellSizeArray();
 
+            // this one
             for (amrex::MFIter mfi(levelset(lev)); mfi.isValid(); ++mfi) {
                 const auto& bx = mfi.growntilebox();
                 const auto& phi = levelset(lev).array(mfi);
diff --git a/amr-wind/mesh_mapping_models/ChannelFlowMap.cpp b/amr-wind/mesh_mapping_models/ChannelFlowMap.cpp
index f56f85651c..843f3a1fc0 100644
--- a/amr-wind/mesh_mapping_models/ChannelFlowMap.cpp
+++ b/amr-wind/mesh_mapping_models/ChannelFlowMap.cpp
@@ -70,6 +70,7 @@ void ChannelFlowMap::create_cell_node_map(int lev, const amrex::Geometry& geom)
         prob_hi[0] - prob_lo[0], prob_hi[1] - prob_lo[1],
         prob_hi[2] - prob_lo[2]};
 
+    // this one
     for (amrex::MFIter mfi((*m_mesh_scale_fac_cc)(lev)); mfi.isValid(); ++mfi) {
 
         const auto& bx = mfi.growntilebox();
@@ -149,6 +150,7 @@ void ChannelFlowMap::create_face_map(int lev, const amrex::Geometry& geom)
         prob_hi[0] - prob_lo[0], prob_hi[1] - prob_lo[1],
         prob_hi[2] - prob_lo[2]};
 
+    // this one
     for (amrex::MFIter mfi((*m_mesh_scale_fac_xf)(lev)); mfi.isValid(); ++mfi) {
 
         const auto& bx = mfi.growntilebox();
@@ -181,6 +183,7 @@ void ChannelFlowMap::create_face_map(int lev, const amrex::Geometry& geom)
             });
     }
 
+    // this one
     for (amrex::MFIter mfi((*m_mesh_scale_fac_yf)(lev)); mfi.isValid(); ++mfi) {
 
         const auto& bx = mfi.growntilebox();
@@ -213,6 +216,7 @@ void ChannelFlowMap::create_face_map(int lev, const amrex::Geometry& geom)
             });
     }
 
+    // this one
     for (amrex::MFIter mfi((*m_mesh_scale_fac_zf)(lev)); mfi.isValid(); ++mfi) {
 
         const auto& bx = mfi.growntilebox();
@@ -277,6 +281,7 @@ void ChannelFlowMap::create_non_uniform_mesh(
         probhi_physical[0] - prob_lo[0], probhi_physical[1] - prob_lo[1],
         probhi_physical[2] - prob_lo[2]};
 
+    // this one
     for (amrex::MFIter mfi((*m_non_uniform_coord_cc)(lev)); mfi.isValid();
          ++mfi) {
 
diff --git a/amr-wind/mesh_mapping_models/ConstantMap.cpp b/amr-wind/mesh_mapping_models/ConstantMap.cpp
index 56292e8d42..0184ccaa36 100644
--- a/amr-wind/mesh_mapping_models/ConstantMap.cpp
+++ b/amr-wind/mesh_mapping_models/ConstantMap.cpp
@@ -27,6 +27,7 @@ void ConstantMap::create_cell_node_map(int lev)
     amrex::Real fac_y = m_fac[1];
     amrex::Real fac_z = m_fac[2];
 
+    // this one
     for (amrex::MFIter mfi((*m_mesh_scale_fac_cc)(lev)); mfi.isValid(); ++mfi) {
 
         const auto& bx = mfi.growntilebox();
@@ -69,6 +70,7 @@ void ConstantMap::create_face_map(int lev)
     amrex::Real fac_y = m_fac[1];
     amrex::Real fac_z = m_fac[2];
 
+    // this one
     for (amrex::MFIter mfi((*m_mesh_scale_fac_xf)(lev)); mfi.isValid(); ++mfi) {
         const auto& bx = mfi.growntilebox();
         amrex::Array4<amrex::Real> const& scale_fac_xf =
@@ -86,6 +88,7 @@ void ConstantMap::create_face_map(int lev)
             });
     }
 
+    // this one
     for (amrex::MFIter mfi((*m_mesh_scale_fac_yf)(lev)); mfi.isValid(); ++mfi) {
         const auto& bx = mfi.growntilebox();
         amrex::Array4<amrex::Real> const& scale_fac_yf =
@@ -103,6 +106,7 @@ void ConstantMap::create_face_map(int lev)
             });
     }
 
+    // this one
     for (amrex::MFIter mfi((*m_mesh_scale_fac_zf)(lev)); mfi.isValid(); ++mfi) {
         const auto& bx = mfi.growntilebox();
         amrex::Array4<amrex::Real> const& scale_fac_zf =
@@ -129,6 +133,7 @@ void ConstantMap::create_non_uniform_mesh(int lev, const amrex::Geometry& geom)
     const auto& problo = geom.ProbLoArray();
     const auto& dx = geom.CellSizeArray();
 
+    // this one
     for (amrex::MFIter mfi((*m_non_uniform_coord_cc)(lev)); mfi.isValid();
          ++mfi) {
 
diff --git a/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H b/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H
index d03a918386..1c9900b842 100644
--- a/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H
+++ b/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H
@@ -44,6 +44,7 @@ struct InitDataOp<LinearWaves>
 
         const auto& problo = geom.ProbLoArray();
         const auto& dx = geom.CellSizeArray();
+        // this one
         for (amrex::MFIter mfi(m_levelset(level)); mfi.isValid(); ++mfi) {
 
             const auto& phi = m_levelset(level).array(mfi);
@@ -132,6 +133,7 @@ struct UpdateRelaxZonesOp<LinearWaves>
             const auto& problo = geom[lev].ProbLoArray();
             const auto& dx = geom[lev].CellSizeArray();
 
+            // this one
             for (amrex::MFIter mfi(m_ow_levelset(lev)); mfi.isValid(); ++mfi) {
                 const auto& phi = m_ow_levelset(lev).array(mfi);
                 const auto& vel = m_ow_velocity(lev).array(mfi);
diff --git a/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.cpp b/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.cpp
index 1b43164083..3a529b7bef 100644
--- a/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.cpp
+++ b/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.cpp
@@ -66,6 +66,7 @@ void apply_relaxation_zones(CFDSim& sim, const RelaxZonesBaseData& wdata)
         auto& target_vof = m_ow_vof(lev);
         const auto& dx = geom[lev].CellSizeArray();
 
+        // this one
         for (amrex::MFIter mfi(ls); mfi.isValid(); ++mfi) {
             const auto& gbx = mfi.growntilebox(2);
             const amrex::Array4<amrex::Real>& phi = ls.array(mfi);
@@ -89,6 +90,7 @@ void apply_relaxation_zones(CFDSim& sim, const RelaxZonesBaseData& wdata)
     auto& density = sim.repo().get_field("density");
 
     for (int lev = 0; lev < nlevels; ++lev) {
+        // this one
         for (amrex::MFIter mfi(vof(lev)); mfi.isValid(); ++mfi) {
             const auto& gbx = mfi.growntilebox(2);
             const auto& dx = geom[lev].CellSizeArray();
diff --git a/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H b/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H
index e0d3ab39e7..9ce7a4d8d0 100644
--- a/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H
+++ b/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H
@@ -83,6 +83,7 @@ struct InitDataOp<StokesWaves>
         auto& velocity = sim.repo().get_field("velocity");
         const auto& problo = geom.ProbLoArray();
         const auto& dx = geom.CellSizeArray();
+        // this one
         for (amrex::MFIter mfi(m_levelset(level)); mfi.isValid(); ++mfi) {
 
             const auto& phi = m_levelset(level).array(mfi);
@@ -158,6 +159,7 @@ struct UpdateRelaxZonesOp<StokesWaves>
             const auto& problo = geom[lev].ProbLoArray();
             const auto& dx = geom[lev].CellSizeArray();
 
+            // this one
             for (amrex::MFIter mfi(m_ow_levelset(lev)); mfi.isValid(); ++mfi) {
                 const auto& phi = m_ow_levelset(lev).array(mfi);
                 const auto& vel = m_ow_velocity(lev).array(mfi);
diff --git a/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H b/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H
index 0c11fca116..1cb0b6465c 100644
--- a/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H
+++ b/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H
@@ -82,6 +82,7 @@ void postprocess_velocity_mfab_liquid(
     amrex::MultiFab& lvs_mfab,
     const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> dx)
 {
+    // this one
     for (amrex::MFIter mfi(vel_mfab); mfi.isValid(); ++mfi) {
         auto vel = vel_mfab.array(mfi);
         auto phi = lvs_mfab.const_array(mfi);
@@ -427,6 +428,7 @@ struct InitDataOp<W2AWaves>
         // Blank initialization if asked for
         if (!init_waves) {
             // Loop to populate field data
+            // this one
             for (amrex::MFIter mfi(levelset(level)); mfi.isValid(); ++mfi) {
                 auto phi = levelset(level).array(mfi);
                 auto vel = velocity(level).array(mfi);
@@ -653,6 +655,7 @@ struct UpdateRelaxZonesOp<W2AWaves>
 
         // Temporally interpolate at every timestep to get target solution
         for (int lev = 0; lev < nlevels; ++lev) {
+            // this one
             for (amrex::MFIter mfi(m_ow_levelset(lev)); mfi.isValid(); ++mfi) {
                 auto phi = m_ow_levelset(lev).array(mfi);
                 auto vel = m_ow_velocity(lev).array(mfi);
diff --git a/amr-wind/overset/OversetOps.cpp b/amr-wind/overset/OversetOps.cpp
index 0ef0671e24..e13857de4a 100644
--- a/amr-wind/overset/OversetOps.cpp
+++ b/amr-wind/overset/OversetOps.cpp
@@ -137,6 +137,7 @@ void OversetOps::update_gradp()
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
+        // this one
         for (MFIter mfi(grad_p(lev), TilingIfNotGPU()); mfi.isValid(); ++mfi) {
             Box const& tbx = mfi.tilebox();
             Array4<Real> const& gp_lev = grad_p(lev).array(mfi);
diff --git a/amr-wind/physics/BoussinesqBubble.cpp b/amr-wind/physics/BoussinesqBubble.cpp
index 444e42134d..c8dc9d67f8 100644
--- a/amr-wind/physics/BoussinesqBubble.cpp
+++ b/amr-wind/physics/BoussinesqBubble.cpp
@@ -32,6 +32,7 @@ void BoussinesqBubble::initialize_fields(int level, const amrex::Geometry& geom)
     auto& density = m_density(level);
     auto& scalars = (*m_temperature)(level);
 
+    // needs openmp pragma?
     for (amrex::MFIter mfi(density); mfi.isValid(); ++mfi) {
         const auto& vbx = mfi.validbox();
 
diff --git a/amr-wind/physics/BurggrafFlow.cpp b/amr-wind/physics/BurggrafFlow.cpp
index 5399048308..c279370207 100644
--- a/amr-wind/physics/BurggrafFlow.cpp
+++ b/amr-wind/physics/BurggrafFlow.cpp
@@ -63,6 +63,7 @@ void BurggrafFlow::initialize_fields(int level, const amrex::Geometry& geom)
 
     UExact u_exact;
     VExact v_exact;
+    // this one
     for (amrex::MFIter mfi(velocity); mfi.isValid(); ++mfi) {
         const auto& vbx = mfi.validbox();
 
@@ -181,4 +182,4 @@ void BurggrafFlow::post_init_actions() { output_error(); }
 
 void BurggrafFlow::post_advance_work() { output_error(); }
 
-} // namespace amr_wind::burggraf
\ No newline at end of file
+} // namespace amr_wind::burggraf
diff --git a/amr-wind/physics/ChannelFlow.cpp b/amr-wind/physics/ChannelFlow.cpp
index f49d99e57e..ce6477b72d 100644
--- a/amr-wind/physics/ChannelFlow.cpp
+++ b/amr-wind/physics/ChannelFlow.cpp
@@ -148,6 +148,7 @@ void ChannelFlow::initialize_fields(
             const auto C0 = m_C0;
             const auto C1 = m_C1;
             const auto dpdx = m_dpdx;
+            // this one
             for (amrex::MFIter mfi(velocity); mfi.isValid(); ++mfi) {
                 const auto& vbx = mfi.validbox();
                 auto vel = velocity.array(mfi);
@@ -181,6 +182,7 @@ void ChannelFlow::initialize_fields(
             }
             auto& walldist = m_repo.get_field("wall_dist")(level);
 
+            // this one
             for (amrex::MFIter mfi(velocity); mfi.isValid(); ++mfi) {
                 const auto& vbx = mfi.validbox();
                 auto vel = velocity.array(mfi);
diff --git a/amr-wind/physics/ConvectingTaylorVortex.cpp b/amr-wind/physics/ConvectingTaylorVortex.cpp
index 08f1822caa..e30499c7c9 100644
--- a/amr-wind/physics/ConvectingTaylorVortex.cpp
+++ b/amr-wind/physics/ConvectingTaylorVortex.cpp
@@ -156,6 +156,7 @@ void ConvectingTaylorVortex::initialize_fields(
     GpyExact gpy_exact;
     GpzExact gpz_exact;
 
+    // this one
     for (amrex::MFIter mfi(velocity); mfi.isValid(); ++mfi) {
         const auto& vbx = mfi.validbox();
 
@@ -240,6 +241,7 @@ amrex::Real ConvectingTaylorVortex::compute_error(const Field& field)
         }
 
         if (m_sim.has_overset()) {
+            // this one
             for (amrex::MFIter mfi(field(lev)); mfi.isValid(); ++mfi) {
                 const auto& vbx = mfi.validbox();
 
diff --git a/amr-wind/physics/EkmanSpiral.cpp b/amr-wind/physics/EkmanSpiral.cpp
index 25dda0ba95..c3bf5491c3 100644
--- a/amr-wind/physics/EkmanSpiral.cpp
+++ b/amr-wind/physics/EkmanSpiral.cpp
@@ -111,6 +111,7 @@ void EkmanSpiral::initialize_fields(int level, const amrex::Geometry& geom)
     UExact u_exact;
     VExact v_exact;
 
+    // this one
     for (amrex::MFIter mfi(velocity); mfi.isValid(); ++mfi) {
         const auto& vbx = mfi.validbox();
 
diff --git a/amr-wind/physics/HybridRANSLESABL.cpp b/amr-wind/physics/HybridRANSLESABL.cpp
index 94a20d2e50..a447365aaa 100644
--- a/amr-wind/physics/HybridRANSLESABL.cpp
+++ b/amr-wind/physics/HybridRANSLESABL.cpp
@@ -20,6 +20,7 @@ void HybridRANSLESABL::initialize_fields(int level, const amrex::Geometry& geom)
     const amrex::Real dz = geom.CellSize()[2];
     const amrex::Real ds = std::cbrt(dx * dy * dz);
 
+    // this one
     for (amrex::MFIter mfi((*m_tke)(level)); mfi.isValid(); ++mfi) {
         const auto& bx = mfi.growntilebox();
         const auto& tke_arr = (*m_tke)(level).array(mfi);
@@ -88,6 +89,7 @@ void HybridRANSLESABL::compute_sdr_impl()
         const amrex::Real dz = geom.CellSize()[2];
         const amrex::Real ds = std::cbrt(dx * dy * dz);
 
+        // this one
         for (amrex::MFIter mfi((*tke)(lev)); mfi.isValid(); ++mfi) {
             const auto& bx = mfi.growntilebox();
             const auto& tke_arr = (*tke)(lev).array(mfi);
diff --git a/amr-wind/physics/RayleighTaylor.cpp b/amr-wind/physics/RayleighTaylor.cpp
index e0756f4fd8..1aa2d14e24 100644
--- a/amr-wind/physics/RayleighTaylor.cpp
+++ b/amr-wind/physics/RayleighTaylor.cpp
@@ -24,6 +24,7 @@ void RayleighTaylor::initialize_fields(int level, const amrex::Geometry& geom)
 
     velocity.setVal(0.0);
 
+    // needs openmp pragma?
     for (amrex::MFIter mfi(density); mfi.isValid(); ++mfi) {
         const auto& vbx = mfi.validbox();
 
diff --git a/amr-wind/physics/ScalarAdvection.cpp b/amr-wind/physics/ScalarAdvection.cpp
index 93a1bcd71b..61206d295c 100644
--- a/amr-wind/physics/ScalarAdvection.cpp
+++ b/amr-wind/physics/ScalarAdvection.cpp
@@ -162,6 +162,7 @@ void ScalarAdvection::initialize_fields(
     auto& density = m_density(level);
     density.setVal(m_rho);
 
+    // this one
     for (amrex::MFIter mfi(velocity); mfi.isValid(); ++mfi) {
         const auto& vbx = mfi.validbox();
         auto vel = velocity.array(mfi);
@@ -215,6 +216,7 @@ void ScalarAdvection::initialize_scalar(const Shape& scalar_function)
         const auto& problo = m_repo.mesh().Geom(level).ProbLoArray();
         auto& scalar = (*m_scalar)(level);
 
+        // this one
         for (amrex::MFIter mfi(scalar); mfi.isValid(); ++mfi) {
             const auto& dx = m_repo.mesh().Geom(level).CellSizeArray();
             const auto& nbx = mfi.nodaltilebox();
diff --git a/amr-wind/physics/SyntheticTurbulence.cpp b/amr-wind/physics/SyntheticTurbulence.cpp
index 2f71d8edac..0bb0bbaafe 100644
--- a/amr-wind/physics/SyntheticTurbulence.cpp
+++ b/amr-wind/physics/SyntheticTurbulence.cpp
@@ -574,6 +574,7 @@ void SyntheticTurbulence::update_impl(
         const auto& gauss_scaling = m_gauss_scaling;
         const auto& epsilon = m_epsilon;
 
+        // this one
         for (amrex::MFIter mfi(m_turb_force(lev)); mfi.isValid(); ++mfi) {
             const auto& bx = mfi.tilebox();
             const auto& turb_force_arr = m_turb_force(lev).array(mfi);
diff --git a/amr-wind/physics/TaylorGreenVortex.cpp b/amr-wind/physics/TaylorGreenVortex.cpp
index 2de400cdda..ac17a6b399 100644
--- a/amr-wind/physics/TaylorGreenVortex.cpp
+++ b/amr-wind/physics/TaylorGreenVortex.cpp
@@ -32,6 +32,7 @@ void TaylorGreenVortex::initialize_fields(
     const amrex::Real Ly = probhi[1] - problo[1];
     const amrex::Real Lz = probhi[2] - problo[2];
 
+    // this one
     for (amrex::MFIter mfi(velocity); mfi.isValid(); ++mfi) {
         const auto& vbx = mfi.validbox();
         const auto& dx = geom.CellSizeArray();
diff --git a/amr-wind/physics/VortexDipole.cpp b/amr-wind/physics/VortexDipole.cpp
index 259b70697c..92f50ce4d0 100644
--- a/amr-wind/physics/VortexDipole.cpp
+++ b/amr-wind/physics/VortexDipole.cpp
@@ -36,6 +36,7 @@ void VortexDipole::initialize_fields(int level, const amrex::Geometry& geom)
 
     const auto& problo = geom.ProbLoArray();
 
+    // this one
     for (amrex::MFIter mfi(velocity); mfi.isValid(); ++mfi) {
         const auto& vbx = mfi.validbox();
         const auto& dx = geom.CellSizeArray();
@@ -74,4 +75,4 @@ void VortexDipole::initialize_fields(int level, const amrex::Geometry& geom)
     }
 }
 
-} // namespace amr_wind
\ No newline at end of file
+} // namespace amr_wind
diff --git a/amr-wind/physics/VortexRing.cpp b/amr-wind/physics/VortexRing.cpp
index 49368852b2..962aefc130 100644
--- a/amr-wind/physics/VortexRing.cpp
+++ b/amr-wind/physics/VortexRing.cpp
@@ -164,6 +164,7 @@ void VortexRing::initialize_velocity(const VortexRingType& vorticity_theta)
 
         const auto& problo = m_repo.mesh().Geom(level).ProbLoArray();
 
+        // this one
         for (amrex::MFIter mfi(m_velocity(level)); mfi.isValid(); ++mfi) {
             const auto& dx = m_repo.mesh().Geom(level).CellSizeArray();
             const auto& nbx = mfi.nodaltilebox();
@@ -250,6 +251,7 @@ void VortexRing::initialize_velocity(const VortexRingType& vorticity_theta)
     for (int level = 0; level <= m_repo.mesh().finestLevel(); ++level) {
         auto& velocity = m_velocity(level);
 
+        // this one
         for (amrex::MFIter mfi(velocity); mfi.isValid(); ++mfi) {
             const auto& vbx = mfi.validbox();
 
diff --git a/amr-wind/physics/multiphase/BreakingWaves.cpp b/amr-wind/physics/multiphase/BreakingWaves.cpp
index de203b33f0..a635135486 100644
--- a/amr-wind/physics/multiphase/BreakingWaves.cpp
+++ b/amr-wind/physics/multiphase/BreakingWaves.cpp
@@ -69,6 +69,7 @@ void BreakingWaves::initialize_fields(int level, const amrex::Geometry& geom)
         sdr.setVal(m_sdr0);
     }
 
+    // this one
     for (amrex::MFIter mfi(levelset); mfi.isValid(); ++mfi) {
         const auto& vbx = mfi.growntilebox();
         auto vel = velocity.array(mfi);
diff --git a/amr-wind/physics/multiphase/DamBreak.cpp b/amr-wind/physics/multiphase/DamBreak.cpp
index 3ffa69ca7e..6bbd3071a7 100644
--- a/amr-wind/physics/multiphase/DamBreak.cpp
+++ b/amr-wind/physics/multiphase/DamBreak.cpp
@@ -42,6 +42,7 @@ void DamBreak::initialize_fields(int level, const amrex::Geometry& geom)
     const amrex::Real rho1 = mphase.rho1();
     const amrex::Real rho2 = mphase.rho2();
 
+    // this one
     for (amrex::MFIter mfi(levelset); mfi.isValid(); ++mfi) {
         const auto& vbx = mfi.validbox();
         auto phi = levelset.array(mfi);
diff --git a/amr-wind/physics/multiphase/MultiPhase.cpp b/amr-wind/physics/multiphase/MultiPhase.cpp
index 4195cc6bb5..17652b381a 100644
--- a/amr-wind/physics/multiphase/MultiPhase.cpp
+++ b/amr-wind/physics/multiphase/MultiPhase.cpp
@@ -304,6 +304,7 @@ void MultiPhase::set_density_via_levelset()
         auto& density = m_density(lev);
         auto& levelset = (*m_levelset)(lev);
 
+        // this one
         for (amrex::MFIter mfi(density); mfi.isValid(); ++mfi) {
             const auto& vbx = mfi.validbox();
             const auto& dx = geom[lev].CellSizeArray();
@@ -342,6 +343,7 @@ void MultiPhase::set_density_via_vof(amr_wind::FieldState fstate)
         auto& density = m_density.state(fstate)(lev);
         auto& vof = (*m_vof).state(fstate)(lev);
 
+        // this one
         for (amrex::MFIter mfi(density); mfi.isValid(); ++mfi) {
             const auto& vbx = mfi.validbox();
             const amrex::Array4<amrex::Real>& F = vof.array(mfi);
@@ -381,6 +383,7 @@ void MultiPhase::calculate_advected_facedensity()
 
     for (int lev = 0; lev < nlevels; ++lev) {
 
+        // this one
         for (amrex::MFIter mfi((*m_vof)(lev)); mfi.isValid(); ++mfi) {
             const auto& bx = mfi.tilebox();
             const auto& bxg1 = amrex::grow(bx, 1);
@@ -429,6 +432,7 @@ void MultiPhase::favre_filtering()
         auto& velocity = m_velocity(lev);
         auto& density = m_density(lev);
 
+        // this one
         for (amrex::MFIter mfi(velocity); mfi.isValid(); ++mfi) {
             const auto& bx = mfi.growntilebox(1);
             const amrex::Array4<amrex::Real>& vel = velocity.array(mfi);
@@ -450,6 +454,7 @@ void MultiPhase::favre_filtering()
         auto& vof = (*m_vof)(lev);
         auto& mom_fil = (*momentum_filter)(lev);
         auto& rho_fil = (*density_filter)(lev);
+        // this one
         for (amrex::MFIter mfi(velocity); mfi.isValid(); ++mfi) {
             const auto& vbx = mfi.validbox();
             const amrex::Array4<amrex::Real>& vel = velocity.array(mfi);
@@ -480,6 +485,7 @@ void MultiPhase::levelset2vof()
         auto& vof = (*m_vof)(lev);
         const auto& dx = geom[lev].CellSizeArray();
 
+        // this one
         for (amrex::MFIter mfi(levelset); mfi.isValid(); ++mfi) {
             const auto& vbx = mfi.validbox();
             const amrex::Array4<amrex::Real>& phi = levelset.array(mfi);
diff --git a/amr-wind/physics/multiphase/RainDrop.cpp b/amr-wind/physics/multiphase/RainDrop.cpp
index fb1f36625b..b4895c5541 100644
--- a/amr-wind/physics/multiphase/RainDrop.cpp
+++ b/amr-wind/physics/multiphase/RainDrop.cpp
@@ -32,6 +32,7 @@ void RainDrop::initialize_fields(int level, const amrex::Geometry& geom)
     const amrex::Real yc = m_loc[1];
     const amrex::Real zc = m_loc[2];
     const amrex::Real radius = m_radius;
+    // this one
     for (amrex::MFIter mfi(velocity); mfi.isValid(); ++mfi) {
         const auto& vbx = mfi.validbox();
         auto phi = levelset.array(mfi);
@@ -69,6 +70,7 @@ void RainDrop::initialize_fields(int level, const amrex::Geometry& geom)
     amrex::Gpu::copy(
         amrex::Gpu::hostToDevice, m_vel.begin(), m_vel.end(), dvel.begin());
     const auto* vptr = dvel.data();
+    // this one
     for (amrex::MFIter mfi(velocity); mfi.isValid(); ++mfi) {
         const auto& vbx = mfi.validbox();
         auto vel = velocity.array(mfi);
diff --git a/amr-wind/physics/multiphase/SloshingTank.cpp b/amr-wind/physics/multiphase/SloshingTank.cpp
index 7a78a31830..45570285fe 100644
--- a/amr-wind/physics/multiphase/SloshingTank.cpp
+++ b/amr-wind/physics/multiphase/SloshingTank.cpp
@@ -34,6 +34,7 @@ void SloshingTank::initialize_fields(int level, const amrex::Geometry& geom)
     const amrex::Real Lx = probhi[0] - problo[0];
     const amrex::Real Ly = probhi[1] - problo[1];
 
+    // this one
     for (amrex::MFIter mfi(levelset); mfi.isValid(); ++mfi) {
         const auto& vbx = mfi.validbox();
         auto phi = levelset.array(mfi);
diff --git a/amr-wind/physics/multiphase/VortexPatch.cpp b/amr-wind/physics/multiphase/VortexPatch.cpp
index 364308bd6f..9527e91654 100644
--- a/amr-wind/physics/multiphase/VortexPatch.cpp
+++ b/amr-wind/physics/multiphase/VortexPatch.cpp
@@ -45,6 +45,7 @@ void VortexPatch::initialize_fields(int level, const amrex::Geometry& geom)
     auto& v_mac = m_sim.repo().get_field("v_mac")(level);
     auto& w_mac = m_sim.repo().get_field("w_mac")(level);
 
+    // this one
     for (amrex::MFIter mfi(velocity); mfi.isValid(); ++mfi) {
         const auto& vbx = mfi.validbox();
         auto uf = u_mac.array(mfi);
@@ -117,6 +118,7 @@ void VortexPatch::pre_advance_work()
         auto& u_mac = m_sim.repo().get_field("u_mac")(lev);
         auto& v_mac = m_sim.repo().get_field("v_mac")(lev);
         auto& w_mac = m_sim.repo().get_field("w_mac")(lev);
+        // this one
         for (amrex::MFIter mfi(m_velocity(lev)); mfi.isValid(); ++mfi) {
             const auto& vbx = mfi.growntilebox(1);
             const auto& dx = geom[lev].CellSizeArray();
@@ -162,6 +164,7 @@ void VortexPatch::post_advance_work()
 
     // Overriding the velocity field
     for (int lev = 0; lev < nlevels; ++lev) {
+        // this one
         for (amrex::MFIter mfi(m_velocity(lev)); mfi.isValid(); ++mfi) {
             const auto& vbx = mfi.growntilebox();
             const auto& dx = geom[lev].CellSizeArray();
diff --git a/amr-wind/physics/multiphase/VortexPatchScalarVel.cpp b/amr-wind/physics/multiphase/VortexPatchScalarVel.cpp
index 24b3184d8c..dabc08c629 100644
--- a/amr-wind/physics/multiphase/VortexPatchScalarVel.cpp
+++ b/amr-wind/physics/multiphase/VortexPatchScalarVel.cpp
@@ -47,6 +47,7 @@ void VortexPatchScalarVel::initialize_fields(
     auto& v_mac = m_sim.repo().get_field("v_mac")(level);
     auto& w_mac = m_sim.repo().get_field("w_mac")(level);
 
+    // this one
     for (amrex::MFIter mfi(velocity); mfi.isValid(); ++mfi) {
         const auto& vbx = mfi.validbox();
         auto uf = u_mac.array(mfi);
@@ -128,6 +129,7 @@ void VortexPatchScalarVel::pre_advance_work()
         auto& u_mac = m_sim.repo().get_field("u_mac")(lev);
         auto& v_mac = m_sim.repo().get_field("v_mac")(lev);
         auto& w_mac = m_sim.repo().get_field("w_mac")(lev);
+        // this one
         for (amrex::MFIter mfi(m_velocity(lev)); mfi.isValid(); ++mfi) {
             const auto& vbx = mfi.growntilebox(1);
             const auto& dx = geom[lev].CellSizeArray();
diff --git a/amr-wind/physics/multiphase/ZalesakDisk.cpp b/amr-wind/physics/multiphase/ZalesakDisk.cpp
index d0139f1478..ac5aa8717a 100644
--- a/amr-wind/physics/multiphase/ZalesakDisk.cpp
+++ b/amr-wind/physics/multiphase/ZalesakDisk.cpp
@@ -50,6 +50,7 @@ void ZalesakDisk::initialize_fields(int level, const amrex::Geometry& geom)
     const amrex::Real hwidth = m_halfwidth;
     const amrex::Real depth = m_depth;
 
+    // this one
     for (amrex::MFIter mfi(levelset); mfi.isValid(); ++mfi) {
         const auto& gbx = mfi.growntilebox(1);
         auto uf = u_mac.array(mfi);
@@ -143,6 +144,7 @@ void ZalesakDisk::pre_advance_work()
         auto& u_mac = m_sim.repo().get_field("u_mac")(lev);
         auto& v_mac = m_sim.repo().get_field("v_mac")(lev);
         auto& w_mac = m_sim.repo().get_field("w_mac")(lev);
+        // this one
         for (amrex::MFIter mfi(m_velocity(lev)); mfi.isValid(); ++mfi) {
             const auto& gbx = mfi.growntilebox(1);
             const auto& dx = geom[lev].CellSizeArray();
@@ -175,6 +177,7 @@ void ZalesakDisk::post_advance_work()
 
     // Overriding the velocity field
     for (int lev = 0; lev < nlevels; ++lev) {
+        // this one
         for (amrex::MFIter mfi(m_velocity(lev)); mfi.isValid(); ++mfi) {
             const auto& vbx = mfi.validbox();
             const auto& dx = geom[lev].CellSizeArray();
diff --git a/amr-wind/physics/multiphase/ZalesakDiskScalarVel.cpp b/amr-wind/physics/multiphase/ZalesakDiskScalarVel.cpp
index fe1cecfa24..94f9f47bc4 100644
--- a/amr-wind/physics/multiphase/ZalesakDiskScalarVel.cpp
+++ b/amr-wind/physics/multiphase/ZalesakDiskScalarVel.cpp
@@ -83,6 +83,7 @@ void ZalesakDiskScalarVel::initialize_fields(
     const amrex::Real hwidth = m_halfwidth;
     const amrex::Real depth = m_depth;
 
+    // this one
     for (amrex::MFIter mfi(levelset); mfi.isValid(); ++mfi) {
         const auto& gbx = mfi.growntilebox(1);
         auto uf = u_mac.array(mfi);
@@ -185,6 +186,7 @@ void ZalesakDiskScalarVel::pre_advance_work()
         auto& u_mac = m_sim.repo().get_field("u_mac")(lev);
         auto& v_mac = m_sim.repo().get_field("v_mac")(lev);
         auto& w_mac = m_sim.repo().get_field("w_mac")(lev);
+        // this one
         for (amrex::MFIter mfi(m_velocity(lev)); mfi.isValid(); ++mfi) {
             const auto& gbx = mfi.growntilebox(1);
             const auto& dx = geom[lev].CellSizeArray();
@@ -250,6 +252,7 @@ amrex::Real ZalesakDiskScalarVel::compute_error(const Field& field)
         }
 
         if (m_sim.has_overset()) {
+            // this one
             for (amrex::MFIter mfi(field(lev)); mfi.isValid(); ++mfi) {
                 const auto& vbx = mfi.validbox();
 
diff --git a/amr-wind/physics/multiphase/hydrostatic_ops.H b/amr-wind/physics/multiphase/hydrostatic_ops.H
index 0cea6e1af8..f965702457 100644
--- a/amr-wind/physics/multiphase/hydrostatic_ops.H
+++ b/amr-wind/physics/multiphase/hydrostatic_ops.H
@@ -15,6 +15,7 @@ static void define_rho0(
     for (int lev = 0; lev < rho0.repo().num_active_levels(); ++lev) {
         const auto& dx = geom[lev].CellSizeArray();
         const auto& problo = geom[lev].ProbLoArray();
+        // this one
         for (amrex::MFIter mfi(rho0(lev)); mfi.isValid(); ++mfi) {
             amrex::Box const& bx = mfi.validbox();
             auto rho0_arr = rho0(lev).array(mfi);
@@ -41,6 +42,7 @@ static void define_p0(
         const auto& dx = geom[lev].CellSizeArray();
         const auto& problo = geom[lev].ProbLoArray();
         const auto& probhi = geom[lev].ProbHiArray();
+        // this one
         for (amrex::MFIter mfi(p0(lev)); mfi.isValid(); ++mfi) {
             amrex::Box const& nbx = mfi.grownnodaltilebox();
             auto p0_arr = p0(lev).array(mfi);
@@ -64,4 +66,4 @@ static void define_p0(
 
 } // namespace amr_wind::hydrostatic
 
-#endif
\ No newline at end of file
+#endif
diff --git a/amr-wind/physics/udfs/UDF.cpp b/amr-wind/physics/udfs/UDF.cpp
index 7a3427de26..52dc3ba1b4 100644
--- a/amr-wind/physics/udfs/UDF.cpp
+++ b/amr-wind/physics/udfs/UDF.cpp
@@ -46,6 +46,7 @@ void UDFImpl<T>::operator()(int level, const amrex::Geometry& geom)
     const amrex::Real time = 0.0;
     const auto ncomp = m_field.num_comp();
     const auto& dop = m_op.device_instance();
+    // this one
     for (amrex::MFIter mfi(mfab); mfi.isValid(); ++mfi) {
         const auto& bx = mfi.tilebox();
         const auto& marr = mfab.array(mfi);
diff --git a/amr-wind/projection/incflo_apply_nodal_projection.cpp b/amr-wind/projection/incflo_apply_nodal_projection.cpp
index be1b8e46d6..b9371c044b 100644
--- a/amr-wind/projection/incflo_apply_nodal_projection.cpp
+++ b/amr-wind/projection/incflo_apply_nodal_projection.cpp
@@ -210,6 +210,7 @@ void incflo::ApplyProjection(
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
+            // this one
             for (MFIter mfi(velocity(lev), TilingIfNotGPU()); mfi.isValid();
                  ++mfi) {
                 Box const& bx = mfi.tilebox();
@@ -252,6 +253,7 @@ void incflo::ApplyProjection(
                     Factory(lev));
                 // At the moment set it to zero
                 surf_tens_force.setVal(0.0);
+                // this one
                 for (MFIter mfi(velocity(lev), TilingIfNotGPU()); mfi.isValid();
                      ++mfi) {
                     Box const& bx = mfi.tilebox();
@@ -301,6 +303,7 @@ void incflo::ApplyProjection(
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
+            // this one
             for (MFIter mfi(sigma[lev], TilingIfNotGPU()); mfi.isValid();
                  ++mfi) {
                 Box const& bx = mfi.tilebox();
@@ -499,6 +502,7 @@ void incflo::ApplyProjection(
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
+        // this one
         for (MFIter mfi(grad_p(lev), TilingIfNotGPU()); mfi.isValid(); ++mfi) {
             Box const& tbx = mfi.tilebox();
             Box const& nbx = mfi.nodaltilebox();
diff --git a/amr-wind/transport_models/TwoPhaseTransport.H b/amr-wind/transport_models/TwoPhaseTransport.H
index 28ee096e9d..2a26efb6c4 100644
--- a/amr-wind/transport_models/TwoPhaseTransport.H
+++ b/amr-wind/transport_models/TwoPhaseTransport.H
@@ -82,6 +82,7 @@ public:
             const auto& vof = m_repo.get_field("vof");
 
             for (int lev = 0; lev < m_repo.num_active_levels(); ++lev) {
+                // this one
                 for (amrex::MFIter mfi((*mu)(lev)); mfi.isValid(); ++mfi) {
                     const auto& vbx = mfi.growntilebox();
                     const auto& volfrac = vof(lev).const_array(mfi);
@@ -103,6 +104,7 @@ public:
             const auto& geom = m_repo.mesh().Geom();
 
             for (int lev = 0; lev < m_repo.num_active_levels(); ++lev) {
+                // this one
                 for (amrex::MFIter mfi((*mu)(lev)); mfi.isValid(); ++mfi) {
                     const auto& vbx = mfi.growntilebox();
                     const auto& dx = geom[lev].CellSizeArray();
@@ -147,6 +149,7 @@ public:
             const auto& vof = m_repo.get_field("vof");
 
             for (int lev = 0; lev < m_repo.num_active_levels(); ++lev) {
+                // this one
                 for (amrex::MFIter mfi((*alpha)(lev)); mfi.isValid(); ++mfi) {
                     const auto& vbx = mfi.growntilebox();
                     const auto& volfrac = vof(lev).const_array(mfi);
@@ -171,6 +174,7 @@ public:
             const auto& geom = m_repo.mesh().Geom();
 
             for (int lev = 0; lev < m_repo.num_active_levels(); ++lev) {
+                // this one
                 for (amrex::MFIter mfi((*alpha)(lev)); mfi.isValid(); ++mfi) {
                     const auto& vbx = mfi.growntilebox();
                     const auto& dx = geom[lev].CellSizeArray();
diff --git a/amr-wind/turbulence/LES/AMD.cpp b/amr-wind/turbulence/LES/AMD.cpp
index 968471ed81..f48ffa038e 100644
--- a/amr-wind/turbulence/LES/AMD.cpp
+++ b/amr-wind/turbulence/LES/AMD.cpp
@@ -140,6 +140,7 @@ void AMD<Transport>::update_alphaeff(Field& alphaeff)
         const auto& geom = geom_vec[lev];
 
         const auto& dx = geom.CellSizeArray();
+        // this one
         for (amrex::MFIter mfi(alphaeff(lev)); mfi.isValid(); ++mfi) {
             const auto& bx = mfi.tilebox();
             const auto& gradVel_arr = (*gradVel)(lev).const_array(mfi);
diff --git a/amr-wind/turbulence/LES/AMDNoTherm.cpp b/amr-wind/turbulence/LES/AMDNoTherm.cpp
index 310aac57f3..2da4c6c748 100644
--- a/amr-wind/turbulence/LES/AMDNoTherm.cpp
+++ b/amr-wind/turbulence/LES/AMDNoTherm.cpp
@@ -50,6 +50,7 @@ void AMDNoTherm<Transport>::update_turbulent_viscosity(
         const auto& geom = geom_vec[lev];
 
         const auto& dx = geom.CellSizeArray();
+        // this one
         for (amrex::MFIter mfi(mu_turb(lev)); mfi.isValid(); ++mfi) {
             const auto& bx = mfi.tilebox();
             const auto& gradVel_arr = (*gradVel)(lev).array(mfi);
diff --git a/amr-wind/turbulence/LES/Kosovic.cpp b/amr-wind/turbulence/LES/Kosovic.cpp
index dfae6d5f0a..914460cab7 100644
--- a/amr-wind/turbulence/LES/Kosovic.cpp
+++ b/amr-wind/turbulence/LES/Kosovic.cpp
@@ -81,6 +81,7 @@ void Kosovic<Transport>::update_turbulent_viscosity(
         const amrex::Real locSurfaceRANSExp = m_surfaceRANSExp;
         const amrex::Real locSurfaceFactor = m_surfaceFactor;
         const amrex::Real locC1 = m_C1;
+        // this one
         for (amrex::MFIter mfi(mu_turb(lev)); mfi.isValid(); ++mfi) {
             const auto& bx = mfi.tilebox();
             const auto& mu_arr = mu_turb(lev).array(mfi);
@@ -140,6 +141,7 @@ void Kosovic<Transport>::update_alphaeff(Field& alphaeff)
     const amrex::Real muCoeff = (m_refMOL < 0) ? 3 : 1;
     const int nlevels = repo.num_active_levels();
     for (int lev = 0; lev < nlevels; ++lev) {
+        // this one
         for (amrex::MFIter mfi(mu_turb(lev)); mfi.isValid(); ++mfi) {
             const auto& bx = mfi.tilebox();
             const auto& muturb_arr = mu_turb(lev).array(mfi);
diff --git a/amr-wind/turbulence/LES/OneEqKsgs.cpp b/amr-wind/turbulence/LES/OneEqKsgs.cpp
index da1c6653e4..83d25ff0eb 100644
--- a/amr-wind/turbulence/LES/OneEqKsgs.cpp
+++ b/amr-wind/turbulence/LES/OneEqKsgs.cpp
@@ -118,6 +118,7 @@ void OneEqKsgsM84<Transport>::update_turbulent_viscosity(
         const amrex::Real dz = geom.CellSize()[2];
         const amrex::Real ds = std::cbrt(dx * dy * dz);
 
+        // this one
         for (amrex::MFIter mfi(mu_turb(lev)); mfi.isValid(); ++mfi) {
             const auto& bx = mfi.tilebox();
             const auto& mu_arr = mu_turb(lev).array(mfi);
@@ -181,6 +182,7 @@ void OneEqKsgsM84<Transport>::update_alphaeff(Field& alphaeff)
         const amrex::Real dz = geom.CellSize()[2];
         const amrex::Real ds = std::cbrt(dx * dy * dz);
 
+        // this one
         for (amrex::MFIter mfi(mu_turb(lev)); mfi.isValid(); ++mfi) {
             const auto& bx = mfi.tilebox();
             const auto& muturb_arr = mu_turb(lev).array(mfi);
@@ -247,6 +249,7 @@ void OneEqKsgsM84<Transport>::post_advance_work()
         const amrex::Real dz = geom.CellSize()[2];
         const amrex::Real ds = std::cbrt(dx * dy * dz);
 
+        // this one
         for (amrex::MFIter mfi(tke(lev)); mfi.isValid(); ++mfi) {
             const auto& bx = mfi.growntilebox();
             const auto& tke_arr = tke(lev).array(mfi);
diff --git a/amr-wind/turbulence/LES/Smagorinsky.cpp b/amr-wind/turbulence/LES/Smagorinsky.cpp
index ec5c75718d..6ea54d3c3d 100644
--- a/amr-wind/turbulence/LES/Smagorinsky.cpp
+++ b/amr-wind/turbulence/LES/Smagorinsky.cpp
@@ -47,6 +47,7 @@ void Smagorinsky<Transport>::update_turbulent_viscosity(
         const amrex::Real ds_sqr = ds * ds;
         const amrex::Real smag_factor = Cs_sqr * ds_sqr;
 
+        // this one
         for (amrex::MFIter mfi(mu_turb(lev)); mfi.isValid(); ++mfi) {
             const auto& bx = mfi.tilebox();
             const auto& mu_arr = mu_turb(lev).array(mfi);
diff --git a/amr-wind/turbulence/RANS/KOmegaSST.cpp b/amr-wind/turbulence/RANS/KOmegaSST.cpp
index 8aab69899a..42e7a5c3e8 100644
--- a/amr-wind/turbulence/RANS/KOmegaSST.cpp
+++ b/amr-wind/turbulence/RANS/KOmegaSST.cpp
@@ -113,6 +113,7 @@ void KOmegaSST<Transport>::update_turbulent_viscosity(
 
     const int nlevels = repo.num_active_levels();
     for (int lev = 0; lev < nlevels; ++lev) {
+        // this one
         for (amrex::MFIter mfi(mu_turb(lev)); mfi.isValid(); ++mfi) {
             const auto& bx = mfi.tilebox();
             const auto& lam_mu_arr = (*lam_mu)(lev).array(mfi);
@@ -281,6 +282,7 @@ void KOmegaSST<Transport>::update_scalar_diff(
         const auto& repo = deff.repo();
         const int nlevels = repo.num_active_levels();
         for (int lev = 0; lev < nlevels; ++lev) {
+            // this one
             for (amrex::MFIter mfi(deff(lev)); mfi.isValid(); ++mfi) {
                 const auto& bx = mfi.tilebox();
                 const auto& lam_mu_arr = (*lam_mu)(lev).array(mfi);
@@ -304,6 +306,7 @@ void KOmegaSST<Transport>::update_scalar_diff(
         const auto& repo = deff.repo();
         const int nlevels = repo.num_active_levels();
         for (int lev = 0; lev < nlevels; ++lev) {
+            // this one
             for (amrex::MFIter mfi(deff(lev)); mfi.isValid(); ++mfi) {
                 const auto& bx = mfi.tilebox();
                 const auto& lam_mu_arr = (*lam_mu)(lev).array(mfi);
diff --git a/amr-wind/turbulence/RANS/KOmegaSSTIDDES.cpp b/amr-wind/turbulence/RANS/KOmegaSSTIDDES.cpp
index cd3dd0ca47..2904f19bb5 100644
--- a/amr-wind/turbulence/RANS/KOmegaSSTIDDES.cpp
+++ b/amr-wind/turbulence/RANS/KOmegaSSTIDDES.cpp
@@ -130,6 +130,7 @@ void KOmegaSSTIDDES<Transport>::update_turbulent_viscosity(
         const amrex::Real hmax =
             amrex::max<amrex::Real>(amrex::max<amrex::Real>(dx, dy), dz);
 
+        // this one
         for (amrex::MFIter mfi(mu_turb(lev)); mfi.isValid(); ++mfi) {
             const auto& bx = mfi.tilebox();
             const auto& lam_mu_arr = (*lam_mu)(lev).array(mfi);
diff --git a/amr-wind/utilities/averaging/ReAveraging.cpp b/amr-wind/utilities/averaging/ReAveraging.cpp
index 10d90edfc3..364dc5cb61 100644
--- a/amr-wind/utilities/averaging/ReAveraging.cpp
+++ b/amr-wind/utilities/averaging/ReAveraging.cpp
@@ -49,6 +49,7 @@ void ReAveraging::operator()(
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
+        // this one
         for (amrex::MFIter mfi(ffab, amrex::TilingIfNotGPU()); mfi.isValid();
              ++mfi) {
             const auto& bx = mfi.tilebox();
diff --git a/amr-wind/utilities/averaging/ReynoldsStress.cpp b/amr-wind/utilities/averaging/ReynoldsStress.cpp
index a503257537..e6d59bc411 100644
--- a/amr-wind/utilities/averaging/ReynoldsStress.cpp
+++ b/amr-wind/utilities/averaging/ReynoldsStress.cpp
@@ -66,6 +66,7 @@ void ReynoldsStress::operator()(
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
+        // this one
         for (amrex::MFIter mfi(ffab, amrex::TilingIfNotGPU()); mfi.isValid();
              ++mfi) {
             const auto& bx = mfi.tilebox();
diff --git a/amr-wind/utilities/sampling/FreeSurfaceSampler.cpp b/amr-wind/utilities/sampling/FreeSurfaceSampler.cpp
index a560eec23f..d7d2250917 100644
--- a/amr-wind/utilities/sampling/FreeSurfaceSampler.cpp
+++ b/amr-wind/utilities/sampling/FreeSurfaceSampler.cpp
@@ -235,6 +235,7 @@ void FreeSurfaceSampler::initialize(const std::string& key)
             geom.ProbLoArray();
         const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> phi =
             geom.ProbHiArray();
+        // this one
         for (amrex::MFIter mfi(floc(lev)); mfi.isValid(); ++mfi) {
             auto loc_arr = floc(lev).array(mfi);
             auto idx_arr = fidx(lev).array(mfi);
@@ -431,6 +432,7 @@ bool FreeSurfaceSampler::update_sampling_locations()
                 geom.InvCellSizeArray();
             const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> plo =
                 geom.ProbLoArray();
+            // this one
             for (amrex::MFIter mfi(floc(lev)); mfi.isValid(); ++mfi) {
                 auto loc_arr = floc(lev).const_array(mfi);
                 auto idx_arr = fidx(lev).const_array(mfi);
@@ -606,6 +608,7 @@ void FreeSurfaceSampler::post_regrid_actions()
             geom.ProbLoArray();
         const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> phi =
             geom.ProbHiArray();
+        // this one
         for (amrex::MFIter mfi(floc(lev)); mfi.isValid(); ++mfi) {
             auto loc_arr = floc(lev).array(mfi);
             auto idx_arr = fidx(lev).array(mfi);
diff --git a/amr-wind/utilities/sampling/SamplingContainer.cpp b/amr-wind/utilities/sampling/SamplingContainer.cpp
index 7fc7e127c4..7930fc57ae 100644
--- a/amr-wind/utilities/sampling/SamplingContainer.cpp
+++ b/amr-wind/utilities/sampling/SamplingContainer.cpp
@@ -19,6 +19,7 @@ void SamplingContainer::setup_container(
 
     const int nlevels = m_mesh.finestLevel() + 1;
     for (int lev = 0; lev < nlevels; ++lev) {
+        // needs openmp pragma?
         for (amrex::MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) {
             DefineAndReturnParticleTile(lev, mfi.index(), mfi.LocalTileIndex());
         }
diff --git a/amr-wind/utilities/tagging/CurvatureRefinement.cpp b/amr-wind/utilities/tagging/CurvatureRefinement.cpp
index 268ffa99ab..f77894b8d2 100644
--- a/amr-wind/utilities/tagging/CurvatureRefinement.cpp
+++ b/amr-wind/utilities/tagging/CurvatureRefinement.cpp
@@ -56,6 +56,7 @@ void CurvatureRefinement::operator()(
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
+    // this one
     for (amrex::MFIter mfi(mfab, amrex::TilingIfNotGPU()); mfi.isValid();
          ++mfi) {
         const auto& bx = mfi.tilebox();
diff --git a/amr-wind/utilities/tagging/FieldRefinement.cpp b/amr-wind/utilities/tagging/FieldRefinement.cpp
index 126eec0b13..9cf92181aa 100644
--- a/amr-wind/utilities/tagging/FieldRefinement.cpp
+++ b/amr-wind/utilities/tagging/FieldRefinement.cpp
@@ -70,6 +70,7 @@ void FieldRefinement::operator()(
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
+    // this one
     for (amrex::MFIter mfi(mfab, amrex::TilingIfNotGPU()); mfi.isValid();
          ++mfi) {
         const auto& bx = mfi.tilebox();
diff --git a/amr-wind/utilities/tagging/GradientMagRefinement.cpp b/amr-wind/utilities/tagging/GradientMagRefinement.cpp
index b1af1ca1d4..a367729753 100644
--- a/amr-wind/utilities/tagging/GradientMagRefinement.cpp
+++ b/amr-wind/utilities/tagging/GradientMagRefinement.cpp
@@ -56,6 +56,7 @@ void GradientMagRefinement::operator()(
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
+    // this one
     for (amrex::MFIter mfi(mfab, amrex::TilingIfNotGPU()); mfi.isValid();
          ++mfi) {
         const auto& bx = mfi.tilebox();
diff --git a/amr-wind/utilities/tagging/OversetRefinement.cpp b/amr-wind/utilities/tagging/OversetRefinement.cpp
index 475603ef72..cf5b7be8ac 100644
--- a/amr-wind/utilities/tagging/OversetRefinement.cpp
+++ b/amr-wind/utilities/tagging/OversetRefinement.cpp
@@ -34,6 +34,7 @@ void OversetRefinement::operator()(
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
+    // this one
     for (amrex::MFIter mfi(ibfab, amrex::TilingIfNotGPU()); mfi.isValid();
          ++mfi) {
         const auto& bx = mfi.tilebox();
diff --git a/amr-wind/utilities/tagging/QCriterionRefinement.cpp b/amr-wind/utilities/tagging/QCriterionRefinement.cpp
index f7bbd4b0d8..20ed7c7001 100644
--- a/amr-wind/utilities/tagging/QCriterionRefinement.cpp
+++ b/amr-wind/utilities/tagging/QCriterionRefinement.cpp
@@ -64,6 +64,7 @@ void QCriterionRefinement::operator()(
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
+    // this one
     for (amrex::MFIter mfi(mfab, amrex::TilingIfNotGPU()); mfi.isValid();
          ++mfi) {
         const auto& bx = mfi.tilebox();
diff --git a/amr-wind/utilities/tagging/VorticityMagRefinement.cpp b/amr-wind/utilities/tagging/VorticityMagRefinement.cpp
index 5544587330..fe1797cf10 100644
--- a/amr-wind/utilities/tagging/VorticityMagRefinement.cpp
+++ b/amr-wind/utilities/tagging/VorticityMagRefinement.cpp
@@ -60,6 +60,7 @@ void VorticityMagRefinement::operator()(
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
+    // this one
     for (amrex::MFIter mfi(mfab, amrex::TilingIfNotGPU()); mfi.isValid();
          ++mfi) {
         const auto& bx = mfi.tilebox();
diff --git a/amr-wind/wind_energy/ABL.cpp b/amr-wind/wind_energy/ABL.cpp
index 7c50d78d4b..c1f703ee4e 100644
--- a/amr-wind/wind_energy/ABL.cpp
+++ b/amr-wind/wind_energy/ABL.cpp
@@ -106,6 +106,7 @@ void ABL::initialize_fields(int level, const amrex::Geometry& geom)
     }
 
     bool interp_fine_levels = false;
+    // needs openmp pragma?
     for (amrex::MFIter mfi(density); mfi.isValid(); ++mfi) {
         const auto& vbx = mfi.validbox();
 
diff --git a/amr-wind/wind_energy/ABLBoundaryPlane.cpp b/amr-wind/wind_energy/ABLBoundaryPlane.cpp
index 17c9db79a4..32ec4742fb 100644
--- a/amr-wind/wind_energy/ABLBoundaryPlane.cpp
+++ b/amr-wind/wind_energy/ABLBoundaryPlane.cpp
@@ -170,6 +170,7 @@ void InletData::read_data_native(
         bndry_n[ori].boxArray(), bndry_n[ori].DistributionMap(),
         bndry_n[ori].nComp(), 0, amrex::MFInfo());
 
+    // needs openmp pragma?
     for (amrex::MFIter mfi(bndry); mfi.isValid(); ++mfi) {
 
         const auto& vbx = mfi.validbox();
@@ -193,6 +194,7 @@ void InletData::read_data_native(
 
     bndry.copyTo((*m_data_n[ori])[lev], 0, nstart, static_cast<int>(nc));
 
+    // needs openmp pragma?
     for (amrex::MFIter mfi(bndry); mfi.isValid(); ++mfi) {
         const auto& vbx = mfi.validbox();
         const auto& bndry_np1_arr = bndry_np1[ori].array(mfi);
diff --git a/amr-wind/wind_energy/ABLFieldInit.cpp b/amr-wind/wind_energy/ABLFieldInit.cpp
index bbddad29b0..b62b4f5fce 100644
--- a/amr-wind/wind_energy/ABLFieldInit.cpp
+++ b/amr-wind/wind_energy/ABLFieldInit.cpp
@@ -299,6 +299,7 @@ void ABLFieldInit::perturb_temperature(
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
+    // this one
     for (amrex::MFIter mfi(theta_fab, amrex::TilingIfNotGPU()); mfi.isValid();
          ++mfi) {
         const auto& bx = mfi.tilebox();
@@ -336,6 +337,7 @@ void ABLFieldInit::init_tke(
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
+    // this one
     for (amrex::MFIter mfi(tke_fab, amrex::TilingIfNotGPU()); mfi.isValid();
          ++mfi) {
         const auto& bx = mfi.tilebox();
diff --git a/amr-wind/wind_energy/ABLModulatedPowerLaw.cpp b/amr-wind/wind_energy/ABLModulatedPowerLaw.cpp
index fc31e62e6e..0c99069426 100644
--- a/amr-wind/wind_energy/ABLModulatedPowerLaw.cpp
+++ b/amr-wind/wind_energy/ABLModulatedPowerLaw.cpp
@@ -173,6 +173,7 @@ void ABLModulatedPowerLaw::set_velocity(
         const auto& dbx = ori.isLow() ? amrex::adjCellLo(domain, idir, nghost)
                                       : amrex::adjCellHi(domain, idir, nghost);
 
+        // needs openmp pragma?
         for (amrex::MFIter mfi(mfab); mfi.isValid(); ++mfi) {
             auto gbx = amrex::grow(mfi.validbox(), nghost);
             const auto& bx =
@@ -248,6 +249,7 @@ void ABLModulatedPowerLaw::set_temperature(
         const auto& dbx = ori.isLow() ? amrex::adjCellLo(domain, idir, nghost)
                                       : amrex::adjCellHi(domain, idir, nghost);
 
+        // needs openmp pragma?
         for (amrex::MFIter mfi(mfab); mfi.isValid(); ++mfi) {
             auto gbx = amrex::grow(mfi.validbox(), nghost);
             const auto& bx =
diff --git a/amr-wind/wind_energy/ABLStats.cpp b/amr-wind/wind_energy/ABLStats.cpp
index 95517285cc..4612859f11 100644
--- a/amr-wind/wind_energy/ABLStats.cpp
+++ b/amr-wind/wind_energy/ABLStats.cpp
@@ -118,6 +118,7 @@ void ABLStats::calc_sfs_stress_avgs(
 
     const int nlevels = repo.num_active_levels();
     for (int lev = 0; lev < nlevels; ++lev) {
+        // this one
         for (amrex::MFIter mfi(m_mueff(lev)); mfi.isValid(); ++mfi) {
             const auto& bx = mfi.tilebox();
             const auto& mueff_arr = m_mueff(lev).array(mfi);
@@ -185,6 +186,7 @@ void ABLStats::calc_tke_diffusion(
 
     const int nlevels = m_sim.repo().num_active_levels();
     for (int lev = 0; lev < nlevels; ++lev) {
+        // this one
         for (amrex::MFIter mfi(diffusion(lev)); mfi.isValid(); ++mfi) {
             const auto& bx = mfi.tilebox();
             const auto& diffusion_arr = diffusion(lev).array(mfi);
diff --git a/amr-wind/wind_energy/actuator/ActuatorContainer.cpp b/amr-wind/wind_energy/actuator/ActuatorContainer.cpp
index 6399d27b08..1b4af4b1f9 100644
--- a/amr-wind/wind_energy/actuator/ActuatorContainer.cpp
+++ b/amr-wind/wind_energy/actuator/ActuatorContainer.cpp
@@ -96,6 +96,7 @@ void ActuatorContainer::initialize_particles(const int total_pts)
     bool assigned = false;
     const int nlevels = m_mesh.finestLevel() + 1;
     for (int lev = 0; ((lev < nlevels) && !assigned); ++lev) {
+        // needs openmp pragma?
         for (auto mfi = MakeMFIter(lev); (mfi.isValid() && !assigned); ++mfi) {
             auto& ptile = GetParticles(
                 lev)[std::make_pair(mfi.index(), mfi.LocalTileIndex())];
diff --git a/unit_tests/aw_test_utils/iter_tools.H b/unit_tests/aw_test_utils/iter_tools.H
index e91f0e4b0b..cd34b1dabd 100644
--- a/unit_tests/aw_test_utils/iter_tools.H
+++ b/unit_tests/aw_test_utils/iter_tools.H
@@ -25,6 +25,7 @@ template <typename Functor>
 void run_algorithm_single_lev(
     const int lev, amrex::MultiFab& lfab, const Functor& func)
 {
+    // needs openmp pragma?
     for (amrex::MFIter mfi(lfab); mfi.isValid(); ++mfi) {
         func(lev, mfi);
     }
@@ -34,6 +35,7 @@ template <typename Functor>
 void run_algorithm_single_lev(
     const int lev, amrex::iMultiFab& lfab, const Functor& func)
 {
+    // needs openmp pragma?
     for (amrex::MFIter mfi(lfab); mfi.isValid(); ++mfi) {
         func(lev, mfi);
     }
diff --git a/unit_tests/core/test_field_ops.cpp b/unit_tests/core/test_field_ops.cpp
index dff707d527..eab71a3951 100644
--- a/unit_tests/core/test_field_ops.cpp
+++ b/unit_tests/core/test_field_ops.cpp
@@ -21,6 +21,7 @@ class FieldOpsTest : public MeshTest
         for (int lev = 0; lev < nlevels; ++lev) {
             const auto& dx = geom[lev].CellSizeArray();
             const auto& problo = geom[lev].ProbLoArray();
+            // this one
             for (amrex::MFIter mfi(field(lev)); mfi.isValid(); ++mfi) {
                 const auto& bx = mfi.tilebox();
                 const auto& field_arr = field(lev).array(mfi);
diff --git a/unit_tests/equation_systems/test_icns_gravityforcing.cpp b/unit_tests/equation_systems/test_icns_gravityforcing.cpp
index 91aaae2e57..9362e62cfb 100644
--- a/unit_tests/equation_systems/test_icns_gravityforcing.cpp
+++ b/unit_tests/equation_systems/test_icns_gravityforcing.cpp
@@ -11,6 +11,7 @@ void init_density(amr_wind::Field& density, const int k_thresh = -3)
 
     for (int lev = 0; lev < nlevels; ++lev) {
 
+        // this one
         for (amrex::MFIter mfi(density(lev)); mfi.isValid(); ++mfi) {
             auto gbx = mfi.growntilebox();
             const auto& darr = density(lev).array(mfi);
@@ -214,4 +215,4 @@ TEST_F(GravityForcingTest, perturb_field)
         Fg, m_nx * m_ny * m_nz, m_nz, amr_wind::FieldState::New, true);
 }
 
-} // namespace amr_wind_tests
\ No newline at end of file
+} // namespace amr_wind_tests
diff --git a/unit_tests/multiphase/test_mflux_schemes.cpp b/unit_tests/multiphase/test_mflux_schemes.cpp
index fa2977d7b5..32735641f6 100644
--- a/unit_tests/multiphase/test_mflux_schemes.cpp
+++ b/unit_tests/multiphase/test_mflux_schemes.cpp
@@ -11,6 +11,7 @@ void init_scalar_increasing(amr_wind::Field& fld, int dir)
     const int nlevels = fld.repo().num_active_levels();
     for (int lev = 0; lev < nlevels; ++lev) {
 
+        // this one
         for (amrex::MFIter mfi(fld(lev)); mfi.isValid(); ++mfi) {
             auto bx = mfi.growntilebox(1);
             const auto& farr = fld(lev).array(mfi);
@@ -28,6 +29,7 @@ void init_scalar_slopechange(amr_wind::Field& fld, int dir, int center)
     const int nlevels = fld.repo().num_active_levels();
     for (int lev = 0; lev < nlevels; ++lev) {
 
+        // this one
         for (amrex::MFIter mfi(fld(lev)); mfi.isValid(); ++mfi) {
             auto bx = mfi.growntilebox(1);
             const auto& farr = fld(lev).array(mfi);
@@ -47,6 +49,7 @@ void init_scalar_uniform(amr_wind::Field& fld, amrex::Real cst)
     const int nlevels = fld.repo().num_active_levels();
     for (int lev = 0; lev < nlevels; ++lev) {
 
+        // this one
         for (amrex::MFIter mfi(fld(lev)); mfi.isValid(); ++mfi) {
             auto bx = mfi.growntilebox(1);
             const auto& farr = fld(lev).array(mfi);
@@ -81,6 +84,7 @@ void get_output_upwind(
 
         auto limiter = PPM::upwind();
 
+        // needs openmp pragma?
         for (amrex::MFIter mfi(fld(lev)); mfi.isValid(); ++mfi) {
             auto bx = mfi.validbox();
             const auto& farr = fld(lev).const_array(mfi);
@@ -127,6 +131,7 @@ void get_output_minmod(
         const auto dlo = amrex::lbound(domain);
         const auto dhi = amrex::ubound(domain);
 
+        // needs openmp pragma?
         for (amrex::MFIter mfi(fld(lev)); mfi.isValid(); ++mfi) {
             auto bx = mfi.validbox();
             const auto& farr = fld(lev).const_array(mfi);
diff --git a/unit_tests/multiphase/test_momflux.cpp b/unit_tests/multiphase/test_momflux.cpp
index bf37babbcf..6f4b8d979e 100644
--- a/unit_tests/multiphase/test_momflux.cpp
+++ b/unit_tests/multiphase/test_momflux.cpp
@@ -18,6 +18,7 @@ void init_field3(
 
     for (int lev = 0; lev < nlevels; ++lev) {
 
+        // this one
         for (amrex::MFIter mfi(fld(lev)); mfi.isValid(); ++mfi) {
             auto bx = mfi.growntilebox();
             const auto& farr = fld(lev).array(mfi);
diff --git a/unit_tests/multiphase/test_vof_BCs.cpp b/unit_tests/multiphase/test_vof_BCs.cpp
index 8311eefbe3..a3378dadfc 100644
--- a/unit_tests/multiphase/test_vof_BCs.cpp
+++ b/unit_tests/multiphase/test_vof_BCs.cpp
@@ -24,6 +24,7 @@ void get_accuracy_vofsol(
     const amr_wind::Field& vof)
 {
     /* -- Check VOF boundary values from fillpatch -- */
+    // needs openmp pragma?
     for (amrex::MFIter mfi(vof(lev)); mfi.isValid(); ++mfi) {
         auto err_arr = err_fld(lev).array(mfi);
         const auto& vof_arr = vof(lev).array(mfi);
@@ -81,6 +82,7 @@ void get_accuracy_advalpha(
     const amr_wind::Field& advalpha_f)
 {
     /* -- Check VOF boundary fluxes -- */
+    // needs openmp pragma?
     for (amrex::MFIter mfi(vof(lev)); mfi.isValid(); ++mfi) {
 
         auto err_arr = err_fld(lev).array(mfi);
diff --git a/unit_tests/ocean_waves/test_relaxation_zones.cpp b/unit_tests/ocean_waves/test_relaxation_zones.cpp
index b9f669b1ce..2cd887e189 100644
--- a/unit_tests/ocean_waves/test_relaxation_zones.cpp
+++ b/unit_tests/ocean_waves/test_relaxation_zones.cpp
@@ -95,6 +95,7 @@ void apply_relaxation_zone_field(
     const auto& geom = comp.repo().mesh().Geom();
 
     for (int lev = 0; lev < comp.repo().num_active_levels(); ++lev) {
+        // this one
         for (amrex::MFIter mfi(comp(lev)); mfi.isValid(); ++mfi) {
             const auto& gbx = mfi.growntilebox(2);
             const auto& dx = geom[lev].CellSizeArray();
diff --git a/unit_tests/turbulence/test_turbulence_LES.cpp b/unit_tests/turbulence/test_turbulence_LES.cpp
index 088b9c97fa..26a8d939ed 100644
--- a/unit_tests/turbulence/test_turbulence_LES.cpp
+++ b/unit_tests/turbulence/test_turbulence_LES.cpp
@@ -21,6 +21,7 @@ void init_field3(amr_wind::Field& fld, amrex::Real srate)
         const auto& dx = mesh.Geom(lev).CellSizeArray();
         const auto& problo = mesh.Geom(lev).ProbLoArray();
 
+        // this one
         for (amrex::MFIter mfi(fld(lev)); mfi.isValid(); ++mfi) {
             auto bx = mfi.growntilebox();
             const auto& farr = fld(lev).array(mfi);
@@ -52,6 +53,7 @@ void init_field_amd(amr_wind::Field& fld, amrex::Real scale)
         const auto& dx = mesh.Geom(lev).CellSizeArray();
         const auto& problo = mesh.Geom(lev).ProbLoArray();
 
+        // this one
         for (amrex::MFIter mfi(fld(lev)); mfi.isValid(); ++mfi) {
             auto bx = mfi.growntilebox();
             const auto& farr = fld(lev).array(mfi);
@@ -83,6 +85,7 @@ void init_field_incomp(amr_wind::Field& fld, amrex::Real scale)
         const auto& dx = mesh.Geom(lev).CellSizeArray();
         const auto& problo = mesh.Geom(lev).ProbLoArray();
 
+        // this one
         for (amrex::MFIter mfi(fld(lev)); mfi.isValid(); ++mfi) {
             auto bx = mfi.growntilebox();
             const auto& farr = fld(lev).array(mfi);
@@ -114,6 +117,7 @@ void init_field1(amr_wind::Field& fld, amrex::Real tgrad)
         const auto& dx = mesh.Geom(lev).CellSizeArray();
         const auto& problo = mesh.Geom(lev).ProbLoArray();
 
+        // this one
         for (amrex::MFIter mfi(fld(lev)); mfi.isValid(); ++mfi) {
             auto bx = mfi.growntilebox();
             const auto& farr = fld(lev).array(mfi);
diff --git a/unit_tests/turbulence/test_turbulence_LES_bc.cpp b/unit_tests/turbulence/test_turbulence_LES_bc.cpp
index 7ec3ece9b2..cc5dbc0a8d 100644
--- a/unit_tests/turbulence/test_turbulence_LES_bc.cpp
+++ b/unit_tests/turbulence/test_turbulence_LES_bc.cpp
@@ -51,6 +51,7 @@ void init_field3(amr_wind::Field& fld, amrex::Real srate)
         const auto& dx = mesh.Geom(lev).CellSizeArray();
         const auto& problo = mesh.Geom(lev).ProbLoArray();
 
+        // this one
         for (amrex::MFIter mfi(fld(lev)); mfi.isValid(); ++mfi) {
             auto bx = mfi.growntilebox();
             const auto& farr = fld(lev).array(mfi);
@@ -80,6 +81,7 @@ void init_field1(amr_wind::Field& fld, amrex::Real tgrad)
         const auto& dx = mesh.Geom(lev).CellSizeArray();
         const auto& problo = mesh.Geom(lev).ProbLoArray();
 
+        // this one
         for (amrex::MFIter mfi(fld(lev)); mfi.isValid(); ++mfi) {
             auto bx = mfi.growntilebox();
             const auto& farr = fld(lev).array(mfi);
diff --git a/unit_tests/utilities/test_diagnostics.cpp b/unit_tests/utilities/test_diagnostics.cpp
index 958d5f0b6f..cff4fb7ba4 100644
--- a/unit_tests/utilities/test_diagnostics.cpp
+++ b/unit_tests/utilities/test_diagnostics.cpp
@@ -14,6 +14,7 @@ void init_velocity(amr_wind::Field& velocity)
         const auto& dx = mesh.Geom(lev).CellSizeArray();
         const auto& problo = mesh.Geom(lev).ProbLoArray();
 
+        // this one
         for (amrex::MFIter mfi(velocity(lev)); mfi.isValid(); ++mfi) {
             auto bx = mfi.validbox();
             const auto& farr = velocity(lev).array(mfi);
@@ -51,6 +52,7 @@ void init_mac_velocity(
         const auto& dx = mesh.Geom(lev).CellSizeArray();
         const auto& problo = mesh.Geom(lev).ProbLoArray();
 
+        // this one
         for (amrex::MFIter mfi(cc(lev)); mfi.isValid(); ++mfi) {
             auto bx = mfi.growntilebox(1);
             const auto& uarr = umac(lev).array(mfi);
diff --git a/unit_tests/utilities/test_field_norms.cpp b/unit_tests/utilities/test_field_norms.cpp
index e2932d2bc3..e768966146 100644
--- a/unit_tests/utilities/test_field_norms.cpp
+++ b/unit_tests/utilities/test_field_norms.cpp
@@ -21,6 +21,7 @@ void init_velocity(
         // Apply a factor to the base level values
         const amrex::Real fac = (lev == 0) ? lev0_fac : 1.0;
 
+        // this one
         for (amrex::MFIter mfi(vel_fld(lev)); mfi.isValid(); ++mfi) {
             auto bx = mfi.growntilebox();
             const auto& farr = vel_fld(lev).array(mfi);
diff --git a/unit_tests/utilities/test_free_surface.cpp b/unit_tests/utilities/test_free_surface.cpp
index d8b7cbc680..eef1ee59b3 100644
--- a/unit_tests/utilities/test_free_surface.cpp
+++ b/unit_tests/utilities/test_free_surface.cpp
@@ -20,6 +20,7 @@ void init_vof(amr_wind::Field& vof_fld, amrex::Real water_level)
         const auto& dx = mesh.Geom(lev).CellSizeArray();
         const auto& problo = mesh.Geom(lev).ProbLoArray();
 
+        // this one
         for (amrex::MFIter mfi(vof_fld(lev)); mfi.isValid(); ++mfi) {
             auto bx = mfi.growntilebox();
             const auto& farr = vof_fld(lev).array(mfi);
@@ -56,6 +57,7 @@ void init_vof_multival(
         const auto& dx = mesh.Geom(lev).CellSizeArray();
         const auto& problo = mesh.Geom(lev).ProbLoArray();
 
+        // this one
         for (amrex::MFIter mfi(vof_fld(lev)); mfi.isValid(); ++mfi) {
             auto bx = mfi.growntilebox();
             const auto& farr = vof_fld(lev).array(mfi);
@@ -108,6 +110,7 @@ void init_vof_slope(
         const auto& dx = mesh.Geom(lev).CellSizeArray();
         const auto& problo = mesh.Geom(lev).ProbLoArray();
 
+        // this one
         for (amrex::MFIter mfi(vof_fld(lev)); mfi.isValid(); ++mfi) {
             auto bx = mfi.growntilebox();
             const auto& farr = vof_fld(lev).array(mfi);
diff --git a/unit_tests/wind_energy/test_abl_init_ncf.cpp b/unit_tests/wind_energy/test_abl_init_ncf.cpp
index a8914fe2e1..4bf51b4357 100644
--- a/unit_tests/wind_energy/test_abl_init_ncf.cpp
+++ b/unit_tests/wind_energy/test_abl_init_ncf.cpp
@@ -130,6 +130,7 @@ TEST_F(ABLMeshTest, abl_init_netcdf_multilevel)
     for (int lev = 0; lev < nlevels; ++lev) {
 
         // Fill base level using input file
+        // needs openmp pragma?
         for (amrex::MFIter mfi(velocityf(lev)); mfi.isValid(); ++mfi) {
             auto vel = velocity[lev]->array(mfi);
             const auto& bx = mfi.validbox();