diff --git a/src/common/include/2dHardcodedIC.fpp b/src/common/include/2dHardcodedIC.fpp
index aec33edb3..5dad945c1 100644
--- a/src/common/include/2dHardcodedIC.fpp
+++ b/src/common/include/2dHardcodedIC.fpp
@@ -252,9 +252,6 @@
         alpha = r/r0
         if (alpha < 1) then
             q_prim_vf(eqn_idx%B%beg)%sf(i, j, 0) = 1._wp/sqrt(4._wp*pi)*(alpha**8 - 2._wp*alpha**4 + 1._wp)
-            ! q_prim_vf(eqn_idx%B%beg)%sf(i,j,0) = 1._wp/sqrt(4000._wp*pi) * (4096._wp*r2**4 - 128._wp*r2**2 + 1._wp)
-            ! q_prim_vf(eqn_idx%B%beg)%sf(i,j,0) = 1._wp/(4._wp*pi) * (alpha**8 - 2._wp*alpha**4 + 1._wp)
-            ! q_prim_vf(eqn_idx%E)%sf(i,j,0) = 6._wp - q_prim_vf(eqn_idx%B%beg)%sf(i,j,0)**2/2._wp
         end if
     case (262)  ! Tilted 2D MHD shock-tube at \alpha = arctan2 (\approx63.4 deg)
         ! rotate by \alpha = atan(2)
diff --git a/src/common/include/macros.fpp b/src/common/include/macros.fpp
index a221fcf8c..e0ede3154 100644
--- a/src/common/include/macros.fpp
+++ b/src/common/include/macros.fpp
@@ -29,23 +29,20 @@
 
         if (nv_uvm_pref_gpu) then
             #:for arg in args
-                ! print*, "Moving ${arg}$ to GPU => ", SHAPE(${arg}$) set preferred location GPU
+                ! set preferred location GPU
                 istat = cudaMemAdvise(c_devloc(${arg}$), SIZEOF(${arg}$), cudaMemAdviseSetPreferredLocation, 0)
                 if (istat /= cudaSuccess) then
                     write (*, "('Error code: ',I0, ': ')") istat
-                    ! write(*,*) cudaGetErrorString(istat)
                 end if
                 ! set accessed by CPU
                 istat = cudaMemAdvise(c_devloc(${arg}$), SIZEOF(${arg}$), cudaMemAdviseSetAccessedBy, cudaCpuDeviceId)
                 if (istat /= cudaSuccess) then
                     write (*, "('Error code: ',I0, ': ')") istat
-                    ! write(*,*) cudaGetErrorString(istat)
                 end if
                 ! prefetch to GPU - physically populate memory pages
                 istat = cudaMemPrefetchAsync(c_devloc(${arg}$), SIZEOF(${arg}$), 0, 0)
                 if (istat /= cudaSuccess) then
                     write (*, "('Error code: ',I0, ': ')") istat
-                    ! write(*,*) cudaGetErrorString(istat)
                 end if
             #:endfor
         end if
diff --git a/src/common/m_phase_change.fpp b/src/common/m_phase_change.fpp
index b56c9a6ec..49bcb963e 100644
--- a/src/common/m_phase_change.fpp
+++ b/src/common/m_phase_change.fpp
@@ -58,14 +58,12 @@ contains
         real(wp) :: rhoe, dynE, rhos      !< total internal energy, kinetic energy, and total entropy
         real(wp) :: rho, rM, m1, m2, MCT  !< total density, total reacting mass, individual reacting masses
         real(wp) :: TvF                   !< total volume fraction
-        ! $:GPU_DECLARE(create='[pS,TS,rhoe,dynE,rhos,rho,rM,m1,m2,MCT,TvF]')
 
         #:if not MFC_CASE_OPTIMIZATION and USING_AMD
             real(wp), dimension(3) :: p_infpT, sk, hk, gk, ek, rhok
         #:else
             real(wp), dimension(num_fluids) :: p_infpT, sk, hk, gk, ek, rhok
         #:endif
-        ! $:GPU_DECLARE(create='[p_infpT,sk,hk,gk,ek,rhok]')
 
         !> Generic loop iterators
         integer :: i, j, k, l
diff --git a/src/post_process/m_data_input.f90 b/src/post_process/m_data_input.f90
index 64b3fa9ff..9617766a2 100644
--- a/src/post_process/m_data_input.f90
+++ b/src/post_process/m_data_input.f90
@@ -40,8 +40,7 @@ end subroutine s_read_abstract_data_files
     type(scalar_field), allocatable, dimension(:), public    :: q_prim_vf  !< Primitive variables
     type(integer_field), allocatable, dimension(:,:), public :: bc_type    !< Boundary condition identifiers
     type(scalar_field), public                               :: q_T_sf     !< Temperature field
-    ! type(scalar_field), public :: ib_markers !<
-    type(integer_field), public :: ib_markers
+    type(integer_field), public                              :: ib_markers
 
     procedure(s_read_abstract_data_files), pointer :: s_read_data_files => null()
 
diff --git a/src/simulation/m_body_forces.fpp b/src/simulation/m_body_forces.fpp
index 3c9519527..9810d3d10 100644
--- a/src/simulation/m_body_forces.fpp
+++ b/src/simulation/m_body_forces.fpp
@@ -14,8 +14,6 @@ module m_body_forces
     use m_helper, only: s_prng, f_unit_vector, f_cross
     use m_nvtx
 
-    ! $:USE_GPU_MODULE()
-
     implicit none
 
     private
diff --git a/src/simulation/m_bubbles_EL.fpp b/src/simulation/m_bubbles_EL.fpp
index 42191ee62..f8166ed9a 100644
--- a/src/simulation/m_bubbles_EL.fpp
+++ b/src/simulation/m_bubbles_EL.fpp
@@ -1758,9 +1758,6 @@ contains
             call my_inquire(trim(file_loc), file_exist)
             if (.not. file_exist) then
                 open (LAG_VOID_ID, FILE=trim(file_loc), form='formatted', position='rewind')
-                ! write (12, *) 'currentTime, averageVoidFraction, ', & 'maximumVoidFraction, totalParticlesVolume' write (12, *)
-                ! 'The averageVoidFraction value does ', & 'not reflect the real void fraction in the cloud since the ', & 'cells
-                ! which do not have bubbles are not accounted'
             else
                 open (LAG_VOID_ID, FILE=trim(file_loc), form='formatted', position='append')
             end if
diff --git a/src/simulation/m_collisions.fpp b/src/simulation/m_collisions.fpp
index fadaa7ece..c60420242 100644
--- a/src/simulation/m_collisions.fpp
+++ b/src/simulation/m_collisions.fpp
@@ -90,8 +90,6 @@ contains
 
         if (num_considered_collisions == 0) return
 
-        ! print *, "Checking Collisions: ", num_considered_collisions, " on rank ", proc_rank
-
         ! Iterate over all collisions detected
         $:GPU_PARALLEL_LOOP(private='[i, l, encoded_pid1, encoded_pid2, xp1, xp2, yp1, yp2, zp1, zp2, pid1, pid2, centroid_1, &
                             & centroid_2, normal_vector, overlap_distance, effective_mass, k, eta, normal_velocity, &
@@ -105,7 +103,6 @@ contains
             pid1 = collision_lookup(i, 1)
             pid2 = collision_lookup(i, 2)
 
-            ! call s_get_neighborhood_idx(pid1, pid1) ! global patch ID -> local index call s_get_neighborhood_idx(pid2, pid2)
             if (pid1 <= 0 .or. pid2 <= 0) cycle
 
             centroid_1(1) = patch_ib(pid1)%x_centroid + real(xp1, wp)*(glb_bounds(1)%end - glb_bounds(1)%beg)
diff --git a/src/simulation/m_compute_cbc.fpp b/src/simulation/m_compute_cbc.fpp
index c2c415b1d..a890130c3 100644
--- a/src/simulation/m_compute_cbc.fpp
+++ b/src/simulation/m_compute_cbc.fpp
@@ -51,7 +51,6 @@ contains
         real(wp), intent(in) :: dpres_ds
         integer              :: i
 
-        ! $:GPU_LOOP(parallelism='[seq]')
         do i = 2, eqn_idx%mom%beg
             L(i) = lambda_factor*lambda2*(c*c*dalpha_rho_ds(i - 1) - mf(i - 1)*dpres_ds)
         end do
@@ -75,7 +74,6 @@ contains
         real(wp), intent(in) :: lambda_factor, lambda2
         integer              :: i
 
-        ! $:GPU_LOOP(parallelism='[seq]')
         do i = eqn_idx%mom%beg + 1, eqn_idx%mom%end
             L(i) = lambda_factor*lambda2*dvel_ds(dir_idx(i - eqn_idx%cont%end))
         end do
@@ -99,7 +97,6 @@ contains
         real(wp), intent(in) :: lambda_factor, lambda2
         integer              :: i
 
-        ! $:GPU_LOOP(parallelism='[seq]')
         do i = eqn_idx%E, eqn_idx%adv%end - 1
             L(i) = lambda_factor*lambda2*dadv_ds(i - eqn_idx%mom%end)
         end do
@@ -125,7 +122,6 @@ contains
 
         if (.not. chemistry) return
 
-        ! $:GPU_LOOP(parallelism='[seq]')
         do i = eqn_idx%species%beg, eqn_idx%species%end
             L(i) = lambda_factor*lambda2*dYs_ds(i - eqn_idx%species%beg + 1)
         end do
diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp
index 8bfc807d7..f91dd84e7 100644
--- a/src/simulation/m_global_parameters.fpp
+++ b/src/simulation/m_global_parameters.fpp
@@ -16,7 +16,6 @@ module m_global_parameters
     use m_helper_basic
     ! Shared state: generated_decls, generated_case_opt_decls, sys_size, eqn_idx, chemistry, shear_*
     use m_global_parameters_common
-    ! $:USE_GPU_MODULE()
 
     implicit none
 
@@ -86,7 +85,6 @@ module m_global_parameters
     logical                :: bodyForces
     real(wp), dimension(3) :: accel_bf
     $:GPU_DECLARE(create='[accel_bf]')
-    ! $:GPU_DECLARE(create='[k_x,w_x,p_x,g_x,k_y,w_y,p_y,g_y,k_z,w_z,p_z,g_z]')
 
     !> Source fields for the spatially supported body force. `spatial_bf` and
     !> `bf_spatial_support` are auto-generated in generated_decls.fpp.
diff --git a/src/simulation/m_riemann_state.fpp b/src/simulation/m_riemann_state.fpp
index abb1733b7..c1daa2204 100644
--- a/src/simulation/m_riemann_state.fpp
+++ b/src/simulation/m_riemann_state.fpp
@@ -941,7 +941,6 @@ contains
                     end if
 
                     if (shear_stress) then
-                        ! current_tau_shear = 0.0_wp
                         call s_calculate_shear_stress_tensor(vel_grad_avg, Re_shear, divergence_v, current_tau_shear)
 
                         do i_dim = 1, num_dims
@@ -955,7 +954,6 @@ contains
                     end if
 
                     if (bulk_stress) then
-                        ! current_tau_bulk = 0.0_wp
                         call s_calculate_bulk_stress_tensor(Re_bulk, divergence_v, current_tau_bulk)
 
                         do i_dim = 1, num_dims
diff --git a/src/simulation/m_weno.fpp b/src/simulation/m_weno.fpp
index 34e8c9b5a..ef18917ab 100644
--- a/src/simulation/m_weno.fpp
+++ b/src/simulation/m_weno.fpp
@@ -10,8 +10,6 @@ module m_weno
     use m_derived_types
     use m_global_parameters
     use m_variables_conversion
-    ! $:USE_GPU_MODULE()
-
     use m_mpi_proxy
     use m_thinc, only: s_thinc_compression
     use m_nvtx