diff --git a/src/pre_process/m_assign_variables.fpp b/src/pre_process/m_assign_variables.fpp
index 5e250310b7..667e360079 100644
--- a/src/pre_process/m_assign_variables.fpp
+++ b/src/pre_process/m_assign_variables.fpp
@@ -97,15 +97,15 @@ contains
 
         do i = 1, eqn_idx%E - eqn_idx%mom%beg
             q_prim_vf(i + 1)%sf(j, k, l) = 1._wp/q_prim_vf(1)%sf(j, k, &
-                      & l)*(eta*patch_icpp(patch_id)%rho*patch_icpp(patch_id)%vel(i) + (1._wp - eta)*patch_icpp(smooth_patch_id) &
-                      & %rho*patch_icpp(smooth_patch_id)%vel(i))
+                      & l)*(eta*patch_icpp(patch_id)%rho*patch_icpp(patch_id)%vel(i) + (1._wp - eta) &
+                      & *patch_icpp(smooth_patch_id)%rho*patch_icpp(smooth_patch_id)%vel(i))
         end do
 
         q_prim_vf(eqn_idx%gamma)%sf(j, k, l) = eta*patch_icpp(patch_id)%gamma + (1._wp - eta)*patch_icpp(smooth_patch_id)%gamma
 
         q_prim_vf(eqn_idx%E)%sf(j, k, l) = 1._wp/q_prim_vf(eqn_idx%gamma)%sf(j, k, &
-                  & l)*(eta*patch_icpp(patch_id)%gamma*patch_icpp(patch_id)%pres + (1._wp - eta)*patch_icpp(smooth_patch_id) &
-                  & %gamma*patch_icpp(smooth_patch_id)%pres)
+                  & l)*(eta*patch_icpp(patch_id)%gamma*patch_icpp(patch_id)%pres + (1._wp - eta) &
+                  & *patch_icpp(smooth_patch_id)%gamma*patch_icpp(smooth_patch_id)%pres)
 
         q_prim_vf(eqn_idx%pi_inf)%sf(j, k, l) = eta*patch_icpp(patch_id)%pi_inf + (1._wp - eta)*patch_icpp(smooth_patch_id)%pi_inf
 
diff --git a/src/pre_process/m_check_ib_patches.fpp b/src/pre_process/m_check_ib_patches.fpp
index 844d051716..db1b01543d 100644
--- a/src/pre_process/m_check_ib_patches.fpp
+++ b/src/pre_process/m_check_ib_patches.fpp
@@ -93,9 +93,9 @@ contains
 
         call s_int_to_str(patch_id, iStr)
 
-        @:PROHIBIT(n == 0 .or. p > 0 .or. patch_ib(patch_id)%length_x <= 0._wp .or. patch_ib(patch_id) &
-                   & %length_y <= 0._wp .or. f_is_default(patch_ib(patch_id)%x_centroid) .or. f_is_default(patch_ib(patch_id) &
-                   & %y_centroid), 'in ellipse IB patch ' // trim(iStr))
+        @:PROHIBIT(n == 0 .or. p > 0 .or. patch_ib(patch_id)%length_x <= 0._wp .or. patch_ib(patch_id)%length_y <= 0._wp &
+                   & .or. f_is_default(patch_ib(patch_id)%x_centroid) .or. f_is_default(patch_ib(patch_id)%y_centroid), &
+                   & 'in ellipse IB patch ' // trim(iStr))
 
     end subroutine s_check_ellipse_ib_patch_geometry
 
@@ -107,9 +107,10 @@ contains
 
         call s_int_to_str(patch_id, iStr)
 
-        @:PROHIBIT(n == 0 .or. p > 0 .or. patch_ib(patch_id)%c <= 0._wp .or. patch_ib(patch_id)%p <= 0._wp .or. patch_ib(patch_id) &
-                   & %t <= 0._wp .or. patch_ib(patch_id)%m <= 0._wp .or. f_is_default(patch_ib(patch_id)%x_centroid) &
-                   & .or. f_is_default(patch_ib(patch_id)%y_centroid), 'in airfoil IB patch ' // trim(iStr))
+        @:PROHIBIT(n == 0 .or. p > 0 .or. patch_ib(patch_id)%c <= 0._wp .or. patch_ib(patch_id)%p <= 0._wp &
+                   & .or. patch_ib(patch_id)%t <= 0._wp .or. patch_ib(patch_id)%m <= 0._wp &
+                   & .or. f_is_default(patch_ib(patch_id)%x_centroid) .or. f_is_default(patch_ib(patch_id)%y_centroid), &
+                   & 'in airfoil IB patch ' // trim(iStr))
 
     end subroutine s_check_airfoil_ib_patch_geometry
 
@@ -121,9 +122,9 @@ contains
 
         call s_int_to_str(patch_id, iStr)
 
-        @:PROHIBIT(n == 0 .or. p == 0 .or. patch_ib(patch_id)%c <= 0._wp .or. patch_ib(patch_id) &
-                   & %p <= 0._wp .or. patch_ib(patch_id)%t <= 0._wp .or. patch_ib(patch_id) &
-                   & %m <= 0._wp .or. f_is_default(patch_ib(patch_id)%x_centroid) .or. f_is_default(patch_ib(patch_id)%y_centroid) &
+        @:PROHIBIT(n == 0 .or. p == 0 .or. patch_ib(patch_id)%c <= 0._wp .or. patch_ib(patch_id)%p <= 0._wp &
+                   & .or. patch_ib(patch_id)%t <= 0._wp .or. patch_ib(patch_id)%m <= 0._wp &
+                   & .or. f_is_default(patch_ib(patch_id)%x_centroid) .or. f_is_default(patch_ib(patch_id)%y_centroid) &
                    & .or. f_is_default(patch_ib(patch_id)%z_centroid) .or. f_is_default(patch_ib(patch_id)%length_z), &
                    & 'in 3d airfoil IB patch ' // trim(iStr))
 
@@ -137,9 +138,9 @@ contains
 
         call s_int_to_str(patch_id, iStr)
 
-        @:PROHIBIT(n == 0 .or. p > 0 .or. f_is_default(patch_ib(patch_id)%x_centroid) .or. f_is_default(patch_ib(patch_id) &
-                   & %y_centroid) .or. patch_ib(patch_id)%length_x <= 0._wp .or. patch_ib(patch_id)%length_y <= 0._wp, &
-                   & 'in rectangle IB patch ' // trim(iStr))
+        @:PROHIBIT(n == 0 .or. p > 0 .or. f_is_default(patch_ib(patch_id)%x_centroid) &
+                   & .or. f_is_default(patch_ib(patch_id)%y_centroid) .or. patch_ib(patch_id)%length_x <= 0._wp &
+                   & .or. patch_ib(patch_id)%length_y <= 0._wp, 'in rectangle IB patch ' // trim(iStr))
 
     end subroutine s_check_rectangle_ib_patch_geometry
 
@@ -151,9 +152,9 @@ contains
 
         call s_int_to_str(patch_id, iStr)
 
-        @:PROHIBIT(n == 0 .or. p == 0 .or. f_is_default(patch_ib(patch_id)%x_centroid) .or. f_is_default(patch_ib(patch_id) &
-                   & %y_centroid) .or. f_is_default(patch_ib(patch_id)%z_centroid) .or. patch_ib(patch_id)%radius <= 0._wp, &
-                   & 'in sphere IB patch ' // trim(iStr))
+        @:PROHIBIT(n == 0 .or. p == 0 .or. f_is_default(patch_ib(patch_id)%x_centroid) &
+                   & .or. f_is_default(patch_ib(patch_id)%y_centroid) .or. f_is_default(patch_ib(patch_id)%z_centroid) &
+                   & .or. patch_ib(patch_id)%radius <= 0._wp, 'in sphere IB patch ' // trim(iStr))
 
     end subroutine s_check_sphere_ib_patch_geometry
 
@@ -165,10 +166,10 @@ contains
 
         call s_int_to_str(patch_id, iStr)
 
-        @:PROHIBIT(n == 0 .or. p == 0 .or. f_is_default(patch_ib(patch_id)%x_centroid) .or. f_is_default(patch_ib(patch_id) &
-                   & %y_centroid) .or. f_is_default(patch_ib(patch_id)%z_centroid) .or. patch_ib(patch_id) &
-                   & %length_x <= 0._wp .or. patch_ib(patch_id)%length_y <= 0._wp .or. patch_ib(patch_id)%length_z <= 0._wp, &
-                   & 'in cuboid IB patch ' // trim(iStr))
+        @:PROHIBIT(n == 0 .or. p == 0 .or. f_is_default(patch_ib(patch_id)%x_centroid) &
+                   & .or. f_is_default(patch_ib(patch_id)%y_centroid) .or. f_is_default(patch_ib(patch_id)%z_centroid) &
+                   & .or. patch_ib(patch_id)%length_x <= 0._wp .or. patch_ib(patch_id)%length_y <= 0._wp &
+                   & .or. patch_ib(patch_id)%length_z <= 0._wp, 'in cuboid IB patch ' // trim(iStr))
 
     end subroutine s_check_cuboid_ib_patch_geometry
 
@@ -181,15 +182,15 @@ contains
         call s_int_to_str(patch_id, iStr)
 
         @:PROHIBIT(p == 0 .or. f_is_default(patch_ib(patch_id)%x_centroid) .or. f_is_default(patch_ib(patch_id)%y_centroid) &
-                   & .or. f_is_default(patch_ib(patch_id)%z_centroid) .or. (patch_ib(patch_id) &
-                   & %length_x <= 0._wp .and. patch_ib(patch_id)%length_y <= 0._wp .and. patch_ib(patch_id)%length_z <= 0._wp) &
+                   & .or. f_is_default(patch_ib(patch_id)%z_centroid) .or. (patch_ib(patch_id)%length_x <= 0._wp &
+                   & .and. patch_ib(patch_id)%length_y <= 0._wp .and. patch_ib(patch_id)%length_z <= 0._wp) &
                    & .or. patch_ib(patch_id)%radius <= 0._wp, 'in cylinder IB patch ' // trim(iStr))
 
         @:PROHIBIT((patch_ib(patch_id)%length_x > 0._wp .and. ((.not. f_is_default(patch_ib(patch_id)%length_y)) &
-                   & .or. (.not. f_is_default(patch_ib(patch_id)%length_z)))) .or. (patch_ib(patch_id) &
-                   & %length_y > 0._wp .and. ((.not. f_is_default(patch_ib(patch_id)%length_x)) &
-                   & .or. (.not. f_is_default(patch_ib(patch_id)%length_z)))) .or. (patch_ib(patch_id) &
-                   & %length_z > 0._wp .and. ((.not. f_is_default(patch_ib(patch_id)%length_x)) &
+                   & .or. (.not. f_is_default(patch_ib(patch_id)%length_z)))) .or. (patch_ib(patch_id)%length_y > 0._wp &
+                   & .and. ((.not. f_is_default(patch_ib(patch_id)%length_x)) &
+                   & .or. (.not. f_is_default(patch_ib(patch_id)%length_z)))) .or. (patch_ib(patch_id)%length_z > 0._wp &
+                   & .and. ((.not. f_is_default(patch_ib(patch_id)%length_x)) &
                    & .or. (.not. f_is_default(patch_ib(patch_id)%length_y)))), 'in cylinder IB patch ' // trim(iStr))
 
     end subroutine s_check_cylinder_ib_patch_geometry
@@ -204,8 +205,8 @@ contains
 
         @:PROHIBIT(patch_ib(patch_id)%model_filepath == dflt_char, 'Empty model file path for patch '//trim(iStr))
 
-        @:PROHIBIT(patch_ib(patch_id)%model_scale(1) <= 0._wp .or. patch_ib(patch_id)%model_scale(2) &
-                   & <= 0._wp .or. patch_ib(patch_id)%model_scale(3) <= 0._wp, 'Negative scale in model IB patch ' // trim(iStr))
+        @:PROHIBIT(patch_ib(patch_id)%model_scale(1) <= 0._wp .or. patch_ib(patch_id)%model_scale(2) <= 0._wp &
+                   & .or. patch_ib(patch_id)%model_scale(3) <= 0._wp, 'Negative scale in model IB patch ' // trim(iStr))
 
     end subroutine s_check_model_ib_patch_geometry
 
@@ -217,8 +218,8 @@ contains
         call s_int_to_str(patch_id, iStr)
 
         @:PROHIBIT((.not. f_is_default(patch_ib(patch_id)%x_centroid)) .or. (.not. f_is_default(patch_ib(patch_id)%y_centroid)) &
-                   & .or. (.not. f_is_default(patch_ib(patch_id)%z_centroid)) .or. (.not. f_is_default(patch_ib(patch_id) &
-                   & %length_x)) .or. (.not. f_is_default(patch_ib(patch_id)%length_y)) &
+                   & .or. (.not. f_is_default(patch_ib(patch_id)%z_centroid)) &
+                   & .or. (.not. f_is_default(patch_ib(patch_id)%length_x)) .or. (.not. f_is_default(patch_ib(patch_id)%length_y)) &
                    & .or. (.not. f_is_default(patch_ib(patch_id)%length_z)) .or. (.not. f_is_default(patch_ib(patch_id)%radius)), &
                    & 'in inactive IB patch ' // trim(iStr))
 
diff --git a/src/pre_process/m_check_patches.fpp b/src/pre_process/m_check_patches.fpp
index adcb233f97..ddad7e5323 100644
--- a/src/pre_process/m_check_patches.fpp
+++ b/src/pre_process/m_check_patches.fpp
@@ -104,10 +104,10 @@ contains
 
         ! Constraints on smoothing initial condition patch parameters
         do i = 1, num_patches
-            if (i > 1 .and. (patch_icpp(i)%geometry == 2 .or. patch_icpp(i)%geometry == 3 .or. patch_icpp(i) &
-                & %geometry == 4 .or. patch_icpp(i)%geometry == 5 .or. patch_icpp(i)%geometry == 8 .or. patch_icpp(i) &
-                & %geometry == 9 .or. patch_icpp(i)%geometry == 10 .or. patch_icpp(i)%geometry == 11 .or. patch_icpp(i) &
-                & %geometry == 12 .or. patch_icpp(i)%geometry == 13 .or. patch_icpp(i)%geometry == 14)) then
+            if (i > 1 .and. (patch_icpp(i)%geometry == 2 .or. patch_icpp(i)%geometry == 3 .or. patch_icpp(i)%geometry == 4 &
+                & .or. patch_icpp(i)%geometry == 5 .or. patch_icpp(i)%geometry == 8 .or. patch_icpp(i)%geometry == 9 &
+                & .or. patch_icpp(i)%geometry == 10 .or. patch_icpp(i)%geometry == 11 .or. patch_icpp(i)%geometry == 12 &
+                & .or. patch_icpp(i)%geometry == 13 .or. patch_icpp(i)%geometry == 14)) then
                 call s_check_supported_patch_smoothing(i)
             else
                 call s_check_unsupported_patch_smoothing(i)
@@ -462,11 +462,13 @@ contains
         call s_int_to_str(patch_id, iStr)
 
         @:PROHIBIT(f_is_default(patch_icpp(patch_id)%vel(1)), "Patch "//trim(iStr)//": vel(1) must be set")
-        @:PROHIBIT(n == 0 .and. (.not. f_is_default(patch_icpp(patch_id)%vel(2))) .and. (.not. f_approx_equal(patch_icpp(patch_id) &
-                   & %vel(2), 0._wp)) .and. (.not. mhd), "Patch " // trim(iStr) // ": vel(2) must not be set when n = 0")
+        @:PROHIBIT(n == 0 .and. (.not. f_is_default(patch_icpp(patch_id)%vel(2))) &
+                   & .and. (.not. f_approx_equal(patch_icpp(patch_id)%vel(2), 0._wp)) .and. (.not. mhd), &
+                   & "Patch " // trim(iStr) // ": vel(2) must not be set when n = 0")
         @:PROHIBIT(n > 0 .and. f_is_default(patch_icpp(patch_id)%vel(2)), "Patch "//trim(iStr)//": vel(2) must be set when n > 0")
-        @:PROHIBIT(p == 0 .and. (.not. f_is_default(patch_icpp(patch_id)%vel(3))) .and. (.not. f_approx_equal(patch_icpp(patch_id) &
-                   & %vel(3), 0._wp)) .and. (.not. mhd), "Patch " // trim(iStr) // ": vel(3) must not be set when p = 0")
+        @:PROHIBIT(p == 0 .and. (.not. f_is_default(patch_icpp(patch_id)%vel(3))) &
+                   & .and. (.not. f_approx_equal(patch_icpp(patch_id)%vel(3), 0._wp)) .and. (.not. mhd), &
+                   & "Patch " // trim(iStr) // ": vel(3) must not be set when p = 0")
         @:PROHIBIT(p > 0 .and. f_is_default(patch_icpp(patch_id)%vel(3)), "Patch "//trim(iStr)//": vel(3) must be set when p > 0")
         @:PROHIBIT(mhd .and. (f_is_default(patch_icpp(patch_id)%vel(2)) .or. f_is_default(patch_icpp(patch_id)%vel(3))), &
                    & "Patch " // trim(iStr) // ": All velocities (vel(1:3)) must be set when mhd = true")
diff --git a/src/pre_process/m_icpp_patches.fpp b/src/pre_process/m_icpp_patches.fpp
index 75104dda37..42e9332c53 100644
--- a/src/pre_process/m_icpp_patches.fpp
+++ b/src/pre_process/m_icpp_patches.fpp
@@ -316,8 +316,9 @@ contains
                                & dy)*(sqrt((x_cc(i) - x_centroid)**2 + (y_cc(j) - y_centroid)**2) - radius))*(-0.5_wp) + 0.5_wp
                 end if
 
-                if (((x_cc(i) - x_centroid)**2 + (y_cc(j) - y_centroid)**2 <= radius**2 .and. patch_icpp(patch_id) &
-                    & %alter_patch(patch_id_fp(i, j, 0))) .or. patch_id_fp(i, j, 0) == smooth_patch_id) then
+                if (((x_cc(i) - x_centroid)**2 + (y_cc(j) - y_centroid)**2 <= radius**2 &
+                    & .and. patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, 0))) .or. patch_id_fp(i, j, &
+                    & 0) == smooth_patch_id) then
                     call s_assign_patch_primitive_variables(patch_id, i, j, 0, eta, q_prim_vf, patch_id_fp)
 
                     @:analytical()
@@ -367,8 +368,8 @@ contains
             do i = 0, m
                 myr = sqrt((x_cc(i) - x_centroid)**2 + (y_cc(j) - y_centroid)**2)
 
-                if (myr <= radius + thickness/2._wp .and. myr >= radius - thickness/2._wp .and. patch_icpp(patch_id) &
-                    & %alter_patch(patch_id_fp(i, j, 0))) then
+                if (myr <= radius + thickness/2._wp .and. myr >= radius - thickness/2._wp &
+                    & .and. patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, 0))) then
                     call s_assign_patch_primitive_variables(patch_id, i, j, 0, eta, q_prim_vf, patch_id_fp)
 
                     @:analytical()
@@ -429,8 +430,8 @@ contains
                 do i = 0, m
                     myr = sqrt((x_cc(i) - x_centroid)**2 + (y_cc(j) - y_centroid)**2)
 
-                    if (myr <= radius + thickness/2._wp .and. myr >= radius - thickness/2._wp .and. patch_icpp(patch_id) &
-                        & %alter_patch(patch_id_fp(i, j, k))) then
+                    if (myr <= radius + thickness/2._wp .and. myr >= radius - thickness/2._wp &
+                        & .and. patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, k))) then
                         call s_assign_patch_primitive_variables(patch_id, i, j, k, eta, q_prim_vf, patch_id_fp)
 
                         @:analytical()
@@ -489,8 +490,9 @@ contains
                                & + 0.5_wp
                 end if
 
-                if ((((x_cc(i) - x_centroid)/a)**2 + ((y_cc(j) - y_centroid)/b)**2 <= 1._wp .and. patch_icpp(patch_id) &
-                    & %alter_patch(patch_id_fp(i, j, 0))) .or. patch_id_fp(i, j, 0) == smooth_patch_id) then
+                if ((((x_cc(i) - x_centroid)/a)**2 + ((y_cc(j) - y_centroid)/b)**2 <= 1._wp &
+                    & .and. patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, 0))) .or. patch_id_fp(i, j, &
+                    & 0) == smooth_patch_id) then
                     call s_assign_patch_primitive_variables(patch_id, i, j, 0, eta, q_prim_vf, patch_id_fp)
 
                     @:analytical()
@@ -558,8 +560,8 @@ contains
                                    & - z_centroid)/c)**2) - 1._wp))*(-0.5_wp) + 0.5_wp
                     end if
 
-                    if ((((x_cc(i) - x_centroid)/a)**2 + ((cart_y - y_centroid)/b)**2 + ((cart_z - z_centroid)/c) &
-                        & **2 <= 1._wp .and. patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, k))) .or. patch_id_fp(i, j, &
+                    if ((((x_cc(i) - x_centroid)/a)**2 + ((cart_y - y_centroid)/b)**2 + ((cart_z - z_centroid)/c)**2 <= 1._wp &
+                        & .and. patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, k))) .or. patch_id_fp(i, j, &
                         & k) == smooth_patch_id) then
                         call s_assign_patch_primitive_variables(patch_id, i, j, k, eta, q_prim_vf, patch_id_fp)
 
@@ -634,8 +636,8 @@ contains
                         if ((q_prim_vf(1)%sf(i, j, 0) < 1.e-10) .and. (model_eqns == 4)) then
                             ! zero density, reassign according to Tait EOS
                             q_prim_vf(1)%sf(i, j, 0) = (((q_prim_vf(eqn_idx%E)%sf(i, j, &
-                                      & 0) + pi_inf)/(pref + pi_inf))**(1._wp/lit_gamma))*rhoref*(1._wp - q_prim_vf(eqn_idx%alf) &
-                                      & %sf(i, j, 0))
+                                      & 0) + pi_inf)/(pref + pi_inf))**(1._wp/lit_gamma))*rhoref*(1._wp &
+                                      & - q_prim_vf(eqn_idx%alf)%sf(i, j, 0))
                         end if
 
                         ! Updating the patch identities bookkeeping variable
@@ -1067,8 +1069,8 @@ contains
                         cart_z = z_cc(k)
                     end if
 
-                    if (x_boundary%beg <= x_cc(i) .and. x_boundary%end >= x_cc(i) &
-                        & .and. y_boundary%beg <= cart_y .and. y_boundary%end >= cart_y .and. z_boundary%beg <= cart_z .and. z_boundary%end >= cart_z) then
+                    if (x_boundary%beg <= x_cc(i) .and. x_boundary%end >= x_cc(i) .and. y_boundary%beg <= cart_y &
+                        & .and. y_boundary%end >= cart_y .and. z_boundary%beg <= cart_z .and. z_boundary%end >= cart_z) then
                         if (patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, k))) then
                             call s_assign_patch_primitive_variables(patch_id, i, j, k, eta, q_prim_vf, patch_id_fp)
 
@@ -1158,12 +1160,11 @@ contains
                         end if
                     end if
 
-                    if (((.not. f_is_default(length_x) .and. (cart_y - y_centroid)**2 + (cart_z - z_centroid) &
-                        & **2 <= radius**2 .and. x_boundary%beg <= x_cc(i) .and. x_boundary%end >= x_cc(i)) &
-                        & .or. (.not. f_is_default(length_y) .and. (x_cc(i) - x_centroid)**2 + (cart_z - z_centroid) &
-                        & **2 <= radius**2 .and. y_boundary%beg <= cart_y .and. y_boundary%end >= cart_y) &
-                        & .or. (.not. f_is_default(length_z) .and. (x_cc(i) - x_centroid)**2 + (cart_y - y_centroid) &
-                        & **2 <= radius**2 .and. z_boundary%beg <= cart_z .and. z_boundary%end >= cart_z) &
+                    if (((.not. f_is_default(length_x) .and. (cart_y - y_centroid)**2 + (cart_z - z_centroid)**2 <= radius**2 &
+                        & .and. x_boundary%beg <= x_cc(i) .and. x_boundary%end >= x_cc(i)) .or. (.not. f_is_default(length_y) &
+                        & .and. (x_cc(i) - x_centroid)**2 + (cart_z - z_centroid)**2 <= radius**2 .and. y_boundary%beg <= cart_y &
+                        & .and. y_boundary%end >= cart_y) .or. (.not. f_is_default(length_z) .and. (x_cc(i) - x_centroid)**2 &
+                        & + (cart_y - y_centroid)**2 <= radius**2 .and. z_boundary%beg <= cart_z .and. z_boundary%end >= cart_z) &
                         & .and. patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, k))) .or. patch_id_fp(i, j, &
                         & k) == smooth_patch_id) then
                         call s_assign_patch_primitive_variables(patch_id, i, j, k, eta, q_prim_vf, patch_id_fp)
diff --git a/src/simulation/m_acoustic_src.fpp b/src/simulation/m_acoustic_src.fpp
index 9b2a969aa4..f5bc800e54 100644
--- a/src/simulation/m_acoustic_src.fpp
+++ b/src/simulation/m_acoustic_src.fpp
@@ -69,8 +69,7 @@ contains
                    & gauss_sigma_dist(1:num_source), gauss_sigma_time(1:num_source), foc_length(1:num_source), &
                    & aperture(1:num_source), npulse(1:num_source), pulse(1:num_source), dir(1:num_source), delay(1:num_source), &
                    & element_polygon_ratio(1:num_source), rotate_angle(1:num_source), element_spacing_angle(1:num_source), &
-                   & num_elements(1:num_source), element_on(1:num_source), bb_num_freq(1:num_source), bb_bandwidth(1:num_source), &
-                                  & bb_lowest_freq(1:num_source))
+                   & num_elements(1:num_source), element_on(1:num_source), bb_num_freq(1:num_source), bb_bandwidth(1:num_source), bb_lowest_freq(1:num_source))
 
         do i = 1, num_source
             do j = 1, 3
@@ -526,8 +525,8 @@ contains
         else if (support(ai) == 2 .or. support(ai) == 3) then  ! 2D or 3D
             ! If we let unit vector e = (cos(dir), sin(dir)),
             dist = r(1)*cos(dir(ai)) + r(2)*sin(dir(ai))  ! dot(r,e)
-            if ((r(1) - dist*cos(dir(ai)))**2._wp + (r(2) - dist*sin(dir(ai)))**2._wp < 0.25_wp*length(ai)**2._wp) &
-                & then  ! |r - dist*e| < length/2
+            ! |r - dist*e| < length/2
+            if ((r(1) - dist*cos(dir(ai)))**2._wp + (r(2) - dist*sin(dir(ai)))**2._wp < 0.25_wp*length(ai)**2._wp) then
                 if (support(ai) /= 3 .or. abs(r(3)) < 0.25_wp*height(ai)) then  ! additional height constraint for 3D
                     source = 1._wp/(sqrt(2._wp*pi)*sig/2._wp)*exp(-0.5_wp*(dist/(sig/2._wp))**2._wp)
                 end if
diff --git a/src/simulation/m_cbc.fpp b/src/simulation/m_cbc.fpp
index d703a3c1e3..e6e14aeb5d 100644
--- a/src/simulation/m_cbc.fpp
+++ b/src/simulation/m_cbc.fpp
@@ -725,15 +725,15 @@ contains
 
                         Ma = vel(dir_idx(1))/c
 
-                        if ((cbc_loc == -1 .and. bc_${XYZ}$%beg == BC_CHAR_SLIP_WALL) &
-                            & .or. (cbc_loc == 1 .and. bc_${XYZ}$%end == BC_CHAR_SLIP_WALL)) then
+                        if ((cbc_loc == -1 .and. bc_${XYZ}$%beg == BC_CHAR_SLIP_WALL) .or. (cbc_loc == 1 &
+                            & .and. bc_${XYZ}$%end == BC_CHAR_SLIP_WALL)) then
                             call s_compute_slip_wall_L(lambda, L, rho, c, dpres_ds, dvel_ds)
-                        else if ((cbc_loc == -1 .and. bc_${XYZ}$%beg == BC_CHAR_NR_SUB_BUFFER) &
-                                 & .or. (cbc_loc == 1 .and. bc_${XYZ}$%end == BC_CHAR_NR_SUB_BUFFER)) then
+                        else if ((cbc_loc == -1 .and. bc_${XYZ}$%beg == BC_CHAR_NR_SUB_BUFFER) .or. (cbc_loc == 1 &
+                                 & .and. bc_${XYZ}$%end == BC_CHAR_NR_SUB_BUFFER)) then
                             call s_compute_nonreflecting_subsonic_buffer_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, &
                                 & dvel_ds, dadv_ds, dYs_ds)
-                        else if ((cbc_loc == -1 .and. bc_${XYZ}$%beg == BC_CHAR_NR_SUB_INFLOW) &
-                                 & .or. (cbc_loc == 1 .and. bc_${XYZ}$%end == BC_CHAR_NR_SUB_INFLOW)) then
+                        else if ((cbc_loc == -1 .and. bc_${XYZ}$%beg == BC_CHAR_NR_SUB_INFLOW) .or. (cbc_loc == 1 &
+                                 & .and. bc_${XYZ}$%end == BC_CHAR_NR_SUB_INFLOW)) then
                             call s_compute_nonreflecting_subsonic_inflow_L(lambda, L, rho, c, dpres_ds, dvel_ds)
                             ! Add GRCBC for Subsonic Inflow
                             if (bc_${XYZ}$%grcbc_in) then
@@ -759,8 +759,8 @@ contains
                                   & dir_idx(1))*sign(1, &
                                   & cbc_loc))/Del_in(${CBC_DIR}$) + c*(1._wp + Ma)*(pres - pres_in(${CBC_DIR}$))/Del_in(${CBC_DIR}$)
                             end if
-                        else if ((cbc_loc == -1 .and. bc_${XYZ}$%beg == BC_CHAR_NR_SUB_OUTFLOW) &
-                                 & .or. (cbc_loc == 1 .and. bc_${XYZ}$%end == BC_CHAR_NR_SUB_OUTFLOW)) then
+                        else if ((cbc_loc == -1 .and. bc_${XYZ}$%beg == BC_CHAR_NR_SUB_OUTFLOW) .or. (cbc_loc == 1 &
+                                 & .and. bc_${XYZ}$%end == BC_CHAR_NR_SUB_OUTFLOW)) then
                             call s_compute_nonreflecting_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, &
                                 & dvel_ds, dadv_ds, dYs_ds)
                             ! Add GRCBC for Subsonic Outflow (Pressure)
@@ -773,19 +773,19 @@ contains
                                       & + vel_out(${CBC_DIR}$, dir_idx(1))*sign(1, cbc_loc))/Del_out(${CBC_DIR}$)
                                 end if
                             end if
-                        else if ((cbc_loc == -1 .and. bc_${XYZ}$%beg == BC_CHAR_FF_SUB_OUTFLOW) &
-                                 & .or. (cbc_loc == 1 .and. bc_${XYZ}$%end == BC_CHAR_FF_SUB_OUTFLOW)) then
+                        else if ((cbc_loc == -1 .and. bc_${XYZ}$%beg == BC_CHAR_FF_SUB_OUTFLOW) .or. (cbc_loc == 1 &
+                                 & .and. bc_${XYZ}$%end == BC_CHAR_FF_SUB_OUTFLOW)) then
                             call s_compute_force_free_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, &
                                 & dadv_ds)
-                        else if ((cbc_loc == -1 .and. bc_${XYZ}$%beg == BC_CHAR_CP_SUB_OUTFLOW) &
-                                 & .or. (cbc_loc == 1 .and. bc_${XYZ}$%end == BC_CHAR_CP_SUB_OUTFLOW)) then
+                        else if ((cbc_loc == -1 .and. bc_${XYZ}$%beg == BC_CHAR_CP_SUB_OUTFLOW) .or. (cbc_loc == 1 &
+                                 & .and. bc_${XYZ}$%end == BC_CHAR_CP_SUB_OUTFLOW)) then
                             call s_compute_constant_pressure_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, &
                                 & dvel_ds, dadv_ds)
-                        else if ((cbc_loc == -1 .and. bc_${XYZ}$%beg == BC_CHAR_SUP_INFLOW) &
-                                 & .or. (cbc_loc == 1 .and. bc_${XYZ}$%end == BC_CHAR_SUP_INFLOW)) then
+                        else if ((cbc_loc == -1 .and. bc_${XYZ}$%beg == BC_CHAR_SUP_INFLOW) .or. (cbc_loc == 1 &
+                                 & .and. bc_${XYZ}$%end == BC_CHAR_SUP_INFLOW)) then
                             call s_compute_supersonic_inflow_L(L)
-                        else if ((cbc_loc == -1 .and. bc_${XYZ}$%beg == BC_CHAR_SUP_OUTFLOW) &
-                                 & .or. (cbc_loc == 1 .and. bc_${XYZ}$%end == BC_CHAR_SUP_OUTFLOW)) then
+                        else if ((cbc_loc == -1 .and. bc_${XYZ}$%beg == BC_CHAR_SUP_OUTFLOW) .or. (cbc_loc == 1 &
+                                 & .and. bc_${XYZ}$%end == BC_CHAR_SUP_OUTFLOW)) then
                             call s_compute_supersonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds, &
                                                                 & dYs_ds)
                         end if
diff --git a/src/simulation/m_compute_levelset.fpp b/src/simulation/m_compute_levelset.fpp
index 96fdebb357..9e7519790c 100644
--- a/src/simulation/m_compute_levelset.fpp
+++ b/src/simulation/m_compute_levelset.fpp
@@ -378,8 +378,8 @@ contains
         xy_local = matmul(inverse_rotation, xy_local)
 
         normal_vector = xy_local
-        normal_vector(2) = normal_vector(2)*(ellipse_coeffs(1)/ellipse_coeffs(2)) &
-                      & **2._wp  ! get the normal direction via the coordinate transformation method
+        ! get the normal direction via the coordinate transformation method
+        normal_vector(2) = normal_vector(2)*(ellipse_coeffs(1)/ellipse_coeffs(2))**2._wp
         normal_vector = normal_vector/sqrt(dot_product(normal_vector, normal_vector))  ! normalize the vector
         gp%levelset_norm = matmul(rotation, normal_vector)  ! save after rotating the vector to the global frame
 
diff --git a/src/simulation/m_hypoelastic.fpp b/src/simulation/m_hypoelastic.fpp
index 471c155f03..390a873377 100644
--- a/src/simulation/m_hypoelastic.fpp
+++ b/src/simulation/m_hypoelastic.fpp
@@ -312,8 +312,7 @@ contains
                         ! S_rr -= rho * v/r * (tau_rr + 2/3*G)
                         rhs_vf(eqn_idx%stress%beg + 2)%sf(k, l, q) = rhs_vf(eqn_idx%stress%beg + 2)%sf(k, l, q) - rho_K_field(k, &
                                & l, q)*q_prim_vf(eqn_idx%mom%beg + 1)%sf(k, l, &
-                               & q)/y_cc(l)*(q_prim_vf(eqn_idx%stress%beg + 2)%sf(k, l, q) + (2._wp/3._wp)*G_K_field(k, l, &
-                               & q))  ! tau_rr + 2/3*G
+                               & q)/y_cc(l)*(q_prim_vf(eqn_idx%stress%beg + 2)%sf(k, l, q) + (2._wp/3._wp)*G_K_field(k, l, q))  ! tau_rr + 2/3*G
 
                         ! S_thetatheta += rho * ( -(tau_thetatheta + 2/3*G)*(du/dx + dv/dr + v/r) + 2*(tau_thetatheta + G)*v/r )
                         rhs_vf(eqn_idx%stress%beg + 3)%sf(k, l, q) = rhs_vf(eqn_idx%stress%beg + 3)%sf(k, l, q) + rho_K_field(k, &
diff --git a/src/simulation/m_ib_patches.fpp b/src/simulation/m_ib_patches.fpp
index 0d3145ddf5..b21483e097 100644
--- a/src/simulation/m_ib_patches.fpp
+++ b/src/simulation/m_ib_patches.fpp
@@ -390,8 +390,8 @@ contains
         do l = ll, lr
             do j = jl, jr
                 do i = il, ir
-                    xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), &
-                                      & z_cc(l) - center(3)]  ! get coordinate frame centered on IB
+                    ! get coordinate frame centered on IB
+                    xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(l) - center(3)]
                     xyz_local = matmul(inverse_rotation, xyz_local)  ! rotate the frame into the IB's coordinates
                     xyz_local = xyz_local - offset  ! airfoils are a patch that require a centroid offset
 
diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp
index 5f865cffeb..39ec6b1470 100644
--- a/src/simulation/m_ibm.fpp
+++ b/src/simulation/m_ibm.fpp
@@ -544,8 +544,7 @@ contains
         if (p == 0) gp_layers_z = 0
 
         $:GPU_PARALLEL_LOOP(private='[i, j, k, ii, jj, kk, is_gp, local_idx, patch_id, encoded_patch_id, xp, yp, zp]', &
-                            & copyin='[count, count_i, x_domain, y_domain, z_domain]', firstprivate='[gp_layers, gp_layers_z]', &
-                                & collapse=3)
+                            & copyin='[count, count_i, x_domain, y_domain, z_domain]', firstprivate='[gp_layers, gp_layers_z]', collapse=3)
         do i = 0, m
             do j = 0, n
                 do k = 0, p
@@ -888,10 +887,10 @@ contains
         type(physical_parameters), dimension(1:num_fluids), intent(in) :: fluid_pp
         integer                                                        :: i, j, k, l, encoded_ib_idx, ib_idx, fluid_idx
         real(wp), dimension(num_ibs, 3)                                :: forces, torques
-        real(wp), dimension(1:3,1:3)                                   :: viscous_stress_div, viscous_stress_div_1, &
-             & viscous_stress_div_2  ! viscous stress tensor with temp vectors to hold divergence calculations
-        real(wp), dimension(1:3) :: local_force_contribution, radial_vector, local_torque_contribution
-        real(wp)                 :: cell_volume, dx, dy, dz, dynamic_viscosity
+        ! viscous stress tensor with temp vectors to hold divergence calculations
+        real(wp), dimension(1:3,1:3) :: viscous_stress_div, viscous_stress_div_1, viscous_stress_div_2
+        real(wp), dimension(1:3)     :: local_force_contribution, radial_vector, local_torque_contribution
+        real(wp)                     :: cell_volume, dx, dy, dz, dynamic_viscosity
 
         #:if not MFC_CASE_OPTIMIZATION and USING_AMD
             real(wp), dimension(3) :: dynamic_viscosities
@@ -941,17 +940,16 @@ contains
                             ! Get the pressure contribution to force via a finite difference to compute the 2D components of the
                             ! gradient of the pressure and cell volume
                             local_force_contribution(1) = local_force_contribution(1) - (q_prim_vf(eqn_idx%E + fluid_idx)%sf(i &
-                                                     & + 1, j, k) - q_prim_vf(eqn_idx%E + fluid_idx)%sf(i - 1, j, &
-                                                     & k))/(2._wp*dx)  ! force is the negative pressure gradient
+                                                     & + 1, j, k) - q_prim_vf(eqn_idx%E + fluid_idx)%sf(i - 1, j, k))/(2._wp*dx)  ! force is the negative pressure gradient
                             local_force_contribution(2) = local_force_contribution(2) - (q_prim_vf(eqn_idx%E + fluid_idx)%sf(i, &
                                                      & j + 1, k) - q_prim_vf(eqn_idx%E + fluid_idx)%sf(i, j - 1, k))/(2._wp*dy)
                             cell_volume = abs(dx*dy)
                             ! add the 3D component of the pressure gradient, if we are working in 3 dimensions
                             if (num_dims == 3) then
                                 dz = z_cc(k + 1) - z_cc(k)
-                                local_force_contribution(3) = local_force_contribution(3) - (q_prim_vf(eqn_idx%E + fluid_idx) &
-                                                         & %sf(i, j, k + 1) - q_prim_vf(eqn_idx%E + fluid_idx)%sf(i, j, &
-                                                         & k - 1))/(2._wp*dz)
+                                local_force_contribution(3) = local_force_contribution(3) - (q_prim_vf(eqn_idx%E &
+                                                         & + fluid_idx)%sf(i, j, k + 1) - q_prim_vf(eqn_idx%E + fluid_idx)%sf(i, &
+                                                         & j, k - 1))/(2._wp*dz)
                                 cell_volume = abs(cell_volume*dz)
                             end if
                         end do
@@ -969,27 +967,27 @@ contains
                             ! get the linear force components first
                             call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i - 1, j, k)
                             call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i + 1, j, k)
-                            viscous_stress_div(1,1:3) = (viscous_stress_div_2(1,1:3) - viscous_stress_div_1(1, &
-                                               & 1:3))/(2._wp*dx)  ! get x derivative of the first-row of viscous stress tensor
-                            local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(1, &
-                                                     & 1:3)  ! add the x components of the divergence to the force
+                            ! get x derivative of the first-row of viscous stress tensor
+                            viscous_stress_div(1,1:3) = (viscous_stress_div_2(1,1:3) - viscous_stress_div_1(1,1:3))/(2._wp*dx)
+                            ! add the x components of the divergence to the force
+                            local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(1,1:3)
 
                             call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i, j - 1, k)
                             call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i, j + 1, k)
-                            viscous_stress_div(2,1:3) = (viscous_stress_div_2(2,1:3) - viscous_stress_div_1(2, &
-                                               & 1:3))/(2._wp*dy)  ! get y derivative of the second-row of viscous stress tensor
-                            local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(2, &
-                                                     & 1:3)  ! add the y components of the divergence to the force
+                            ! get y derivative of the second-row of viscous stress tensor
+                            viscous_stress_div(2,1:3) = (viscous_stress_div_2(2,1:3) - viscous_stress_div_1(2,1:3))/(2._wp*dy)
+                            ! add the y components of the divergence to the force
+                            local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(2,1:3)
 
                             if (num_dims == 3) then
                                 call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i, j, &
                                                                      & k - 1)
                                 call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i, j, &
                                                                      & k + 1)
-                                viscous_stress_div(3,1:3) = (viscous_stress_div_2(3,1:3) - viscous_stress_div_1(3, &
-                                                   & 1:3))/(2._wp*dz)  ! get z derivative of the third-row of viscous stress tensor
-                                local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(3, &
-                                                         & 1:3)  ! add the z components of the divergence to the force
+                                ! get z derivative of the third-row of viscous stress tensor
+                                viscous_stress_div(3,1:3) = (viscous_stress_div_2(3,1:3) - viscous_stress_div_1(3,1:3))/(2._wp*dz)
+                                ! add the z components of the divergence to the force
+                                local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(3,1:3)
                             end if
                         end if
 
@@ -1060,8 +1058,8 @@ contains
 
         ! Offset only needs to be computes for specific geometries
 
-        if (patch_ib(ib_marker)%geometry == 4 .or. patch_ib(ib_marker)%geometry == 5 .or. patch_ib(ib_marker) &
-            & %geometry == 11 .or. patch_ib(ib_marker)%geometry == 12) then
+        if (patch_ib(ib_marker)%geometry == 4 .or. patch_ib(ib_marker)%geometry == 5 .or. patch_ib(ib_marker)%geometry == 11 &
+            & .or. patch_ib(ib_marker)%geometry == 12) then
             center_of_mass_local = [0._wp, 0._wp, 0._wp]
             num_cells_local = 0
 
@@ -1130,11 +1128,11 @@ contains
         if (patch_ib(ib_marker)%geometry == 2) then  ! circle
             patch_ib(ib_marker)%moment = 0.5_wp*patch_ib(ib_marker)%mass*(patch_ib(ib_marker)%radius)**2
         else if (patch_ib(ib_marker)%geometry == 3) then  ! rectangle
-            patch_ib(ib_marker)%moment = patch_ib(ib_marker)%mass*(patch_ib(ib_marker)%length_x**2 + patch_ib(ib_marker) &
-                     & %length_y**2)/6._wp
+            patch_ib(ib_marker)%moment = patch_ib(ib_marker)%mass*(patch_ib(ib_marker)%length_x**2 &
+                     & + patch_ib(ib_marker)%length_y**2)/6._wp
         else if (patch_ib(ib_marker)%geometry == 6) then  ! ellipse
-            patch_ib(ib_marker)%moment = 0.0625_wp*patch_ib(ib_marker)%mass*(patch_ib(ib_marker)%length_x**2 + patch_ib(ib_marker) &
-                     & %length_y**2)
+            patch_ib(ib_marker)%moment = 0.0625_wp*patch_ib(ib_marker)%mass*(patch_ib(ib_marker)%length_x**2 &
+                     & + patch_ib(ib_marker)%length_y**2)
         else if (patch_ib(ib_marker)%geometry == 8) then  ! sphere
             patch_ib(ib_marker)%moment = 0.4*patch_ib(ib_marker)%mass*(patch_ib(ib_marker)%radius)**2
         else  ! we do not have an analytic moment of inertia calculation and need to approximate it directly via a sum
diff --git a/src/simulation/m_qbmm.fpp b/src/simulation/m_qbmm.fpp
index 9f0fceecdf..8465d4d349 100644
--- a/src/simulation/m_qbmm.fpp
+++ b/src/simulation/m_qbmm.fpp
@@ -402,12 +402,11 @@ contains
         type(scalar_field), dimension(sys_size), intent(in)                                        :: flux_n_vf
         real(stp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), intent(inout) :: pb
 
-        real(wp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), &
-             & intent(inout) :: rhs_pb  ! TODO :: I think that this should be stp as well.
-
-        integer  :: i, j, k, l, q
+        ! TODO :: I think that this should be stp as well.
+        real(wp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), intent(inout) :: rhs_pb
+        integer :: i, j, k, l, q
         real(wp) :: nb_q, nb_dot, R, R2, nR, nR2, nR_dot, nR2_dot, var, AX
-        logical  :: is_axisym
+        logical :: is_axisym
 
         select case (idir)
         case (1)
@@ -841,8 +840,8 @@ contains
                                                                             & momrhs(:,i1, i2, j, q))
                                                 end if
                                             case (2)
-                                                if ((j >= 7 .and. j <= 9) .or. (j >= 22 .and. j <= 23) &
-                                                    & .or. (j >= 10 .and. j <= 11) .or. (j == 26)) then
+                                                if ((j >= 7 .and. j <= 9) .or. (j >= 22 .and. j <= 23) .or. (j >= 10 &
+                                                    & .and. j <= 11) .or. (j == 26)) then
                                                     momsum = momsum + coeff(j, i1, i2)*(R0(q)**momrhs(3, i1, i2, j, &
                                                                             & q))*f_quad2D(abscX(:,q), abscY(:,q), wght_pb(:,q), &
                                                                             & momrhs(:,i1, i2, j, q))
diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp
index 60235675ea..4316c9278c 100644
--- a/src/simulation/m_rhs.fpp
+++ b/src/simulation/m_rhs.fpp
@@ -517,9 +517,8 @@ contains
         type(scalar_field), dimension(sys_size), intent(inout)                                     :: rhs_vf
         real(stp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), intent(inout) :: pb_in
 
-        real(wp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), &
-             & intent(inout) &
-             & :: rhs_pb  ! TODO :: I think these other two variables need to be stp as well, but it doesn't compile like that right now
+        ! TODO :: I think these other two variables need to be stp as well, but it doesn't compile like that right now
+        real(wp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), intent(inout) :: rhs_pb
         real(stp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), intent(inout) :: mv_in
         real(wp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), intent(inout) :: rhs_mv
         integer, intent(in) :: t_step
@@ -530,6 +529,7 @@ contains
         integer(kind=8) :: i, j, k, l, q  !< Generic loop iterators
 
         ! RHS: halo exchange -> reconstruct -> Riemann solve -> flux difference -> source terms
+
         call nvtxStartRange("COMPUTE-RHS")
 
         call cpu_time(t_start)
diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp
index b5e618438d..78a00759ad 100644
--- a/src/simulation/m_riemann_solvers.fpp
+++ b/src/simulation/m_riemann_solvers.fpp
@@ -423,11 +423,11 @@ contains
                                     pres_mag%R = 0.5_wp*(B%R(1)**2._wp + B%R(2)**2._wp + B%R(3)**2._wp)
                                 #:endif
                                 E_L = gamma_L*pres_L + pi_inf_L + 0.5_wp*rho_L*vel_L_rms + qv_L + pres_mag%L
-                                E_R = gamma_R*pres_R + pi_inf_R + 0.5_wp*rho_R*vel_R_rms + qv_R &
-                                    & + pres_mag%R  ! includes magnetic energy
+                                ! includes magnetic energy
+                                E_R = gamma_R*pres_R + pi_inf_R + 0.5_wp*rho_R*vel_R_rms + qv_R + pres_mag%R
                                 H_L = (E_L + pres_L - pres_mag%L)/rho_L
-                                H_R = (E_R + pres_R - pres_mag%R) &
-                                       & /rho_R  ! stagnation enthalpy here excludes magnetic energy (only used to find speed of sound)
+                                ! stagnation enthalpy here excludes magnetic energy (only used to find speed of sound)
+                                H_R = (E_R + pres_R - pres_mag%R)/rho_R
                             else
                                 E_L = gamma_L*pres_L + pi_inf_L + 5.e-1*rho_L*vel_L_rms + qv_L
                                 E_R = gamma_R*pres_R + pi_inf_R + 5.e-1*rho_R*vel_R_rms + qv_R
@@ -754,9 +754,8 @@ contains
                                                           & eqn_idx%psi) - qR_prim_rs${XYZ}$_vf(j + 1, k, l, &
                                                           & eqn_idx%psi)))/(s_M - s_P)
                                     else
-                                        flux_rs${XYZ}$_vf(j, k, l, &
-                                                          & eqn_idx%B%beg + norm_dir - 1) &
-                                                          & = 0._wp  ! Without hyperbolic cleaning, make sure flux of B_normal is identically zero
+                                        ! Without hyperbolic cleaning, make sure flux of B_normal is identically zero
+                                        flux_rs${XYZ}$_vf(j, k, l, eqn_idx%B%beg + norm_dir - 1) = 0._wp
                                     end if
                                 end if
                                 flux_src_rs${XYZ}$_vf(j, k, l, eqn_idx%adv%beg) = 0._wp
@@ -1113,11 +1112,11 @@ contains
                                 pres_mag%L = 0.5_wp*(B%L(1)**2._wp + B%L(2)**2._wp + B%L(3)**2._wp)
                                 pres_mag%R = 0.5_wp*(B%R(1)**2._wp + B%R(2)**2._wp + B%R(3)**2._wp)
                                 E_L = gamma_L*pres_L + pi_inf_L + 0.5_wp*rho_L*vel_L_rms + qv_L + pres_mag%L
-                                E_R = gamma_R*pres_R + pi_inf_R + 0.5_wp*rho_R*vel_R_rms + qv_R &
-                                    & + pres_mag%R  ! includes magnetic energy
+                                ! includes magnetic energy
+                                E_R = gamma_R*pres_R + pi_inf_R + 0.5_wp*rho_R*vel_R_rms + qv_R + pres_mag%R
                                 H_L = (E_L + pres_L - pres_mag%L)/rho_L
-                                H_R = (E_R + pres_R - pres_mag%R) &
-                                       & /rho_R  ! stagnation enthalpy here excludes magnetic energy (only used to find speed of sound)
+                                ! stagnation enthalpy here excludes magnetic energy (only used to find speed of sound)
+                                H_R = (E_R + pres_R - pres_mag%R)/rho_R
                             else
                                 E_L = gamma_L*pres_L + pi_inf_L + 5.e-1*rho_L*vel_L_rms + qv_L
                                 E_R = gamma_R*pres_R + pi_inf_R + 5.e-1*rho_R*vel_R_rms + qv_R
@@ -3474,8 +3473,8 @@ contains
                             E%L = gamma%L*pres%L + pi_inf%L + 0.5_wp*rho%L*vel_rms%L + qv%L + pres_mag%L
                             E%R = gamma%R*pres%R + pi_inf%R + 0.5_wp*rho%R*vel_rms%R + qv%R + pres_mag%R  ! includes magnetic energy
                             H_no_mag%L = (E%L + pres%L - pres_mag%L)/rho%L
-                            H_no_mag%R = (E%R + pres%R - pres_mag%R) &
-                                          & /rho%R  ! stagnation enthalpy here excludes magnetic energy (only used to find speed of sound)
+                            ! stagnation enthalpy here excludes magnetic energy (only used to find speed of sound)
+                            H_no_mag%R = (E%R + pres%R - pres_mag%R)/rho%R
 
                             ! (2) Compute fast wave speeds
                             call s_compute_speed_of_sound(pres%L, rho%L, gamma%L, pi_inf%L, H_no_mag%L, alpha_L, vel_rms%L, &
diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp
index ce0ae48c4b..130053fce5 100644
--- a/src/simulation/m_time_steppers.fpp
+++ b/src/simulation/m_time_steppers.fpp
@@ -130,8 +130,8 @@ contains
         pool_dims(4) = sys_size
         pool_starts(4) = 1
 #ifdef MFC_MIXED_PRECISION
-        pool_size = 1_8*(idwbuff(1)%end - idwbuff(1)%beg + 1)*(idwbuff(2)%end - idwbuff(2)%beg + 1)*(idwbuff(3)%end - idwbuff(3) &
-                         & %beg + 1)*sys_size
+        pool_size = 1_8*(idwbuff(1)%end - idwbuff(1)%beg + 1)*(idwbuff(2)%end - idwbuff(2)%beg + 1)*(idwbuff(3)%end &
+                         & - idwbuff(3)%beg + 1)*sys_size
         call hipCheck(hipMalloc_(cptr_device, pool_size*2_8))
         call c_f_pointer(cptr_device, q_cons_ts_pool_device, shape=pool_dims)
         q_cons_ts_pool_device(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:) => q_cons_ts_pool_device
@@ -748,8 +748,8 @@ contains
                              & 3)*dt*patch_ib(i)%torque/rk_coef(s, 4))  ! add the torque to the angular momentum
                     call s_compute_moment_of_inertia(i, patch_ib(i)%angular_vel)
                     ! update the moment of inertia to be based on the direction of the angular momentum
-                    patch_ib(i)%angular_vel = patch_ib(i)%angular_vel/patch_ib(i) &
-                             & %moment  ! convert back to angular velocity with the new moment of inertia
+                    ! convert back to angular velocity with the new moment of inertia
+                    patch_ib(i)%angular_vel = patch_ib(i)%angular_vel/patch_ib(i)%moment
                 end if
 
                 ! Update the angle of the IB
diff --git a/src/simulation/m_weno.fpp b/src/simulation/m_weno.fpp
index 0e6f83ae2f..5ccd85756f 100644
--- a/src/simulation/m_weno.fpp
+++ b/src/simulation/m_weno.fpp
@@ -1204,8 +1204,8 @@ contains
                                         delta(:) = 0._wp
                                         beta(:) = weno_eps
 
-                                        if (teno) v = v_rs_ws_${XYZ}$ (j - 3:j + 3,k, l, &
-                                            & i)  ! temporary field value array for clarity
+                                        ! temporary field value array for clarity
+                                        if (teno) v = v_rs_ws_${XYZ}$ (j - 3:j + 3,k, l, i)
 
                                         if (.not. teno) then
                                             dvd(2) = v_rs_ws_${XYZ}$ (j + 3, k, l, i) - v_rs_ws_${XYZ}$ (j + 2, k, l, i)
@@ -1305,8 +1305,8 @@ contains
                                             tau = abs(beta(3) - beta(0))  ! Equation 50
                                             $:GPU_LOOP(parallelism='[seq]')
                                             do q = 0, weno_num_stencils
-                                                alpha(q) = d_cbL_${XYZ}$ (q, &
-                                                      & j)*(1._wp + (tau/beta(q))**wenoz_q)  ! wenoz_q = 2,3,4 for stability
+                                                ! wenoz_q = 2,3,4 for stability
+                                                alpha(q) = d_cbL_${XYZ}$ (q, j)*(1._wp + (tau/beta(q))**wenoz_q)
                                             end do
                                         else if (teno) then
                                             #:if not MFC_CASE_OPTIMIZATION or weno_num_stencils > 3
@@ -1379,8 +1379,8 @@ contains
                                         else if (wenoz) then
                                             $:GPU_LOOP(parallelism='[seq]')
                                             do q = 0, weno_num_stencils
-                                                alpha(q) = d_cbR_${XYZ}$ (q, &
-                                                      & j)*(1._wp + (tau/beta(q))**wenoz_q)  ! wenoz_q = 2,3,4 for stability
+                                                ! wenoz_q = 2,3,4 for stability
+                                                alpha(q) = d_cbR_${XYZ}$ (q, j)*(1._wp + (tau/beta(q))**wenoz_q)
                                             end do
                                         else if (teno) then
                                             $:GPU_LOOP(parallelism='[seq]')
diff --git a/toolchain/pyproject.toml b/toolchain/pyproject.toml
index fe9b7b7fee..0d85419f9c 100644
--- a/toolchain/pyproject.toml
+++ b/toolchain/pyproject.toml
@@ -24,7 +24,7 @@ dependencies = [
     # Code Health
     "typos",
     "ruff",
-    "ffmt",
+    "ffmt==0.4.0",
     "ansi2txt",
 
     # Profiling
