diff --git a/src/simulation/m_body_forces.fpp b/src/simulation/m_body_forces.fpp
index 97d66f6ce6..490f0f45bb 100644
--- a/src/simulation/m_body_forces.fpp
+++ b/src/simulation/m_body_forces.fpp
@@ -19,9 +19,10 @@ module m_body_forces
 
     implicit none
 
-    private; public :: s_compute_body_forces_rhs, &
- s_initialize_body_forces_module, &
- s_finalize_body_forces_module
+    private; 
+    public :: s_compute_body_forces_rhs, &
+              s_initialize_body_forces_module, &
+              s_finalize_body_forces_module
 
 #ifdef CRAY_ACC_WAR
     @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), rhoM)
@@ -35,7 +36,7 @@ contains
 
     !> This subroutine inializes the module global array of mixture
     !! densities in each grid cell
-    subroutine s_initialize_body_forces_module()
+    subroutine s_initialize_body_forces_module
 
         ! Simulation is at least 2D
         if (n > 0) then
@@ -62,7 +63,7 @@ contains
     !> This subroutine computes the acceleration at time t
     subroutine s_compute_acceleration(t)
 
-        real(kind(0d0)) :: t
+        real(kind(0d0)), intent(in) :: t
 
         if (m > 0) then
             accel_bf(1) = g_x + k_x*sin(w_x*t - p_x)
@@ -80,9 +81,10 @@ contains
 
     !> This subroutine calculates the mixture density at each cell
     !! center
+    !! param q_cons_vf Conservative variable
     subroutine s_compute_mixture_density(q_cons_vf)
 
-        type(scalar_field), dimension(sys_size), intent(IN) :: q_cons_vf
+        type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf
         integer :: i, j, k, l !< standard iterators
 
         !$acc parallel loop collapse(3) gang vector default(present)
@@ -102,11 +104,13 @@ contains
 
     !> This subroutine calculates the source term due to body forces
     !! so the system can be advanced in time
+    !! @param q_cons_vf Conservative variables
+    !! @param q_prim_vf Primitive variables
     subroutine s_compute_body_forces_rhs(q_cons_vf, q_prim_vf, rhs_vf)
 
-        type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf
-        type(scalar_field), dimension(sys_size), intent(IN) :: q_cons_vf
-        type(scalar_field), dimension(sys_size), intent(INOUT) :: rhs_vf
+        type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
+        type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf
+        type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf
 
         integer :: i, j, k, l !< Loop variables
 
@@ -172,7 +176,7 @@ contains
 
     end subroutine s_compute_body_forces_rhs
 
-    subroutine s_finalize_body_forces_module()
+    subroutine s_finalize_body_forces_module
 
         @:DEALLOCATE_GLOBAL(rhoM)
 
diff --git a/src/simulation/m_boundary_conditions.fpp b/src/simulation/m_boundary_conditions.fpp
index 14b41f8545..ca77632dc1 100644
--- a/src/simulation/m_boundary_conditions.fpp
+++ b/src/simulation/m_boundary_conditions.fpp
@@ -18,18 +18,21 @@ module m_boundary_conditions
 
     implicit none
 
-    private; public :: s_populate_primitive_variables_buffers, &
- s_populate_capillary_buffers
+    private; 
+    public :: s_populate_primitive_variables_buffers, &
+              s_populate_capillary_buffers
 
 contains
 
     !>  The purpose of this procedure is to populate the buffers
-    !!      of the conservative variables, depending on the selected
+    !!      of the primitive variables, depending on the selected
     !!      boundary conditions.
+    !! @param q_prim_vf Primitive variable
     subroutine s_populate_primitive_variables_buffers(q_prim_vf, pb, mv)
 
-        type(scalar_field), dimension(sys_size) :: q_prim_vf
-        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(INOUT) :: pb, mv
+        type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
+        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: pb, mv
+
         integer :: bc_loc, bc_dir
 
         ! Population of Buffers in x-direction =============================
@@ -214,9 +217,9 @@ contains
 
     subroutine s_ghost_cell_extrapolation(q_prim_vf, pb, mv, bc_dir, bc_loc)
 
-        type(scalar_field), dimension(sys_size) :: q_prim_vf
-        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(INOUT) :: pb, mv
-        integer :: bc_dir, bc_loc
+        type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
+        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: pb, mv
+        integer, intent(in) :: bc_dir, bc_loc
         integer :: j, k, l, q, i
 
         !< x-direction =========================================================
@@ -325,9 +328,10 @@ contains
 
     subroutine s_symmetry(q_prim_vf, pb, mv, bc_dir, bc_loc)
 
-        type(scalar_field), dimension(sys_size) :: q_prim_vf
-        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(INOUT) :: pb, mv
-        integer :: bc_dir, bc_loc
+        type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
+        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: pb, mv
+        integer, intent(in) :: bc_dir, bc_loc
+
         integer :: j, k, l, q, i
 
         !< x-direction =========================================================
@@ -606,9 +610,10 @@ contains
 
     subroutine s_periodic(q_prim_vf, pb, mv, bc_dir, bc_loc)
 
-        type(scalar_field), dimension(sys_size) :: q_prim_vf
-        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(INOUT) :: pb, mv
-        integer :: bc_dir, bc_loc
+        type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
+        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: pb, mv
+        integer, intent(in) :: bc_dir, bc_loc
+
         integer :: j, k, l, q, i
 
         !< x-direction =========================================================
@@ -825,9 +830,10 @@ contains
 
     subroutine s_axis(q_prim_vf, pb, mv, bc_dir, bc_loc)
 
-        type(scalar_field), dimension(sys_size) :: q_prim_vf
-        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(INOUT) :: pb, mv
-        integer :: bc_dir, bc_loc
+        type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
+        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: pb, mv
+        integer, intent(in) :: bc_dir, bc_loc
+
         integer :: j, k, l, q, i
 
         !$acc parallel loop collapse(3) gang vector default(present)
@@ -897,9 +903,10 @@ contains
 
     subroutine s_slip_wall(q_prim_vf, pb, mv, bc_dir, bc_loc)
 
-        type(scalar_field), dimension(sys_size) :: q_prim_vf
-        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(INOUT) :: pb, mv
-        integer :: bc_dir, bc_loc
+        type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
+        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: pb, mv
+        integer, intent(in) :: bc_dir, bc_loc
+
         integer :: j, k, l, q, i
 
         !< x-direction =========================================================
@@ -1038,9 +1045,10 @@ contains
 
     subroutine s_no_slip_wall(q_prim_vf, pb, mv, bc_dir, bc_loc)
 
-        type(scalar_field), dimension(sys_size) :: q_prim_vf
-        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(INOUT) :: pb, mv
-        integer :: bc_dir, bc_loc
+        type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
+        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: pb, mv
+        integer, intent(in) :: bc_dir, bc_loc
+
         integer :: j, k, l, q, i
 
         !< x-direction =========================================================
@@ -1215,8 +1223,9 @@ contains
 
     subroutine s_qbmm_extrapolation(pb, mv, bc_dir, bc_loc)
 
-        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(INOUT) :: pb, mv
-        integer :: bc_dir, bc_loc
+        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: pb, mv
+        integer, intent(in) :: bc_dir, bc_loc
+
         integer :: j, k, l, q, i
 
         !< x-direction =========================================================
@@ -1348,7 +1357,7 @@ contains
 
     subroutine s_populate_capillary_buffers(c_divs)
 
-        type(scalar_field), dimension(num_dims + 1) :: c_divs
+        type(scalar_field), dimension(num_dims + 1), intent(inout) :: c_divs
         integer :: i, j, k, l
 
         ! x - direction
diff --git a/src/simulation/m_bubbles.fpp b/src/simulation/m_bubbles.fpp
index 279719bbc3..4e36ca6ae1 100644
--- a/src/simulation/m_bubbles.fpp
+++ b/src/simulation/m_bubbles.fpp
@@ -55,7 +55,7 @@ module m_bubbles
 
 contains
 
-    subroutine s_initialize_bubbles_module()
+    subroutine s_initialize_bubbles_module
 
         integer :: i, j, k, l, q
         type(int_bounds_info) :: ix, iy, iz
@@ -103,7 +103,7 @@ contains
     ! Compute the bubble volume fraction alpha from the bubble number density n
         !! @param q_cons_vf is the conservative variable
     subroutine s_comp_alpha_from_n(q_cons_vf)
-        type(scalar_field), dimension(sys_size), intent(INOUT) :: q_cons_vf
+        type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
         real(kind(0d0)) :: nR3bar
         integer(kind(0d0)) :: i, j, k, l
 
@@ -125,8 +125,9 @@ contains
 
     subroutine s_compute_bubbles_rhs(idir, q_prim_vf)
 
-        type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf
-        integer, intent(IN) :: idir
+        integer, intent(in) :: idir
+        type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
+
         integer :: i, j, k, l, q
 
         if (idir == 1) then
@@ -182,17 +183,12 @@ contains
         !!      that are needed for the bubble modeling
         !!  @param q_prim_vf Primitive variables
         !!  @param q_cons_vf Conservative variables
-        !!  @param divu Divergence of velocity
-        !!  @param bub_adv_src Advection equation source due to bubble compression/expansion
-        !!  @param bub_r_src   Bubble radius equation source
-        !!  @param bub_v_src   Bubble velocity equation source
-        !!  @param bub_p_src   Bubble pressure equation source
-        !!  @param bub_m_src   Bubble mass equation source
     subroutine s_compute_bubble_source(q_cons_vf, q_prim_vf, t_step, rhs_vf)
 
-        type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf
-        type(scalar_field), dimension(sys_size), intent(INOUT) :: rhs_vf, q_cons_vf
-        integer, intent(IN) :: t_step
+        type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
+        type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
+        integer, intent(in) :: t_step
+        type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf
 
         real(kind(0d0)) :: tmp1, tmp2, tmp3, tmp4, &
                            c_gas, c_liquid, &
@@ -497,7 +493,7 @@ contains
         !!  @param fpb Internal bubble pressure
     function f_cpbw(fR0, fR, fV, fpb)
         !$acc routine seq
-        real(kind(0d0)), intent(IN) :: fR0, fR, fV, fpb
+        real(kind(0d0)), intent(in) :: fR0, fR, fV, fpb
 
         real(kind(0d0)) :: f_cpbw
 
@@ -516,7 +512,7 @@ contains
         !!  @param fBtait Tait EOS parameter
     function f_H(fCpbw, fCpinf, fntait, fBtait)
         !$acc routine seq
-        real(kind(0d0)), intent(IN) :: fCpbw, fCpinf, fntait, fBtait
+        real(kind(0d0)), intent(in) :: fCpbw, fCpinf, fntait, fBtait
 
         real(kind(0d0)) :: tmp1, tmp2, tmp3
         real(kind(0d0)) :: f_H
@@ -536,7 +532,7 @@ contains
         !! @param fH Bubble enthalpy
     function f_cgas(fCpinf, fntait, fBtait, fH)
         !$acc routine seq
-        real(kind(0d0)), intent(IN) :: fCpinf, fntait, fBtait, fH
+        real(kind(0d0)), intent(in) :: fCpinf, fntait, fBtait, fH
 
         real(kind(0d0)) :: tmp
         real(kind(0d0)) :: f_cgas
@@ -559,7 +555,7 @@ contains
         !!  @param divu Divergence of velocity
     function f_cpinfdot(fRho, fP, falf, fntait, fBtait, advsrc, divu)
         !$acc routine seq
-        real(kind(0d0)), intent(IN) :: fRho, fP, falf, fntait, fBtait, advsrc, divu
+        real(kind(0d0)), intent(in) :: fRho, fP, falf, fntait, fBtait, advsrc, divu
 
         real(kind(0d0)) :: c2_liquid
         real(kind(0d0)) :: f_cpinfdot
@@ -583,14 +579,14 @@ contains
         !!  @param fCpinf_dot Time derivative of the driving pressure
         !!  @param fntait Tait EOS parameter
         !!  @param fBtait Tait EOS parameter
-        !!  @param fR0 Equilibrium bubble radius
         !!  @param fR Current bubble radius
         !!  @param fV Current bubble velocity
+        !!  @param fR0 Equilibrium bubble radius
         !!  @param fpbdot Time derivative of the internal bubble pressure
     function f_Hdot(fCpbw, fCpinf, fCpinf_dot, fntait, fBtait, fR, fV, fR0, fpbdot)
         !$acc routine seq
-        real(kind(0d0)), intent(IN) :: fCpbw, fCpinf, fCpinf_dot, fntait, fBtait
-        real(kind(0d0)), intent(IN) :: fR, fV, fR0, fpbdot
+        real(kind(0d0)), intent(in) :: fCpbw, fCpinf, fCpinf_dot, fntait, fBtait
+        real(kind(0d0)), intent(in) :: fR, fV, fR0, fpbdot
 
         real(kind(0d0)) :: tmp1, tmp2
         real(kind(0d0)) :: f_Hdot
@@ -631,8 +627,9 @@ contains
         !!  @param f_divu Divergence of velocity
     function f_rddot(fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, fntait, fBtait, f_bub_adv_src, f_divu)
         !$acc routine seq
-        real(kind(0d0)), intent(IN) :: fRho, fP, fR, fV, fR0, fpb, fpbdot, alf
-        real(kind(0d0)), intent(IN) :: fntait, fBtait, f_bub_adv_src, f_divu
+        real(kind(0d0)), intent(in) :: fRho, fP, fR, fV, fR0, fpb, fpbdot, alf
+        real(kind(0d0)), intent(in) :: fntait, fBtait, f_bub_adv_src, f_divu
+
         real(kind(0d0)) :: fCpbw, fCpinf, fCpinf_dot, fH, fHdot, c_gas, c_liquid
         real(kind(0d0)) :: f_rddot
 
@@ -668,7 +665,8 @@ contains
         !!  @param fCpbw Boundary wall pressure
     function f_rddot_RP(fCp, fRho, fR, fV, fR0, fCpbw)
         !$acc routine seq
-        real(kind(0d0)), intent(IN) :: fCp, fRho, fR, fV, fR0, fCpbw
+        real(kind(0d0)), intent(in) :: fCp, fRho, fR, fV, fR0, fCpbw
+
         real(kind(0d0)) :: f_rddot_RP
 
             !! rddot = (1/r) (  -3/2 rdot^2 + ((r0/r)^3\gamma - Cp)/rho )
@@ -690,8 +688,8 @@ contains
         !!  @param fBtait Tait EOS parameter
     function f_rddot_G(fCpbw, fR, fV, fH, fHdot, fcgas, fntait, fBtait)
         !$acc routine seq
-        real(kind(0d0)), intent(IN) :: fCpbw, fR, fV, fH, fHdot
-        real(kind(0d0)), intent(IN) :: fcgas, fntait, fBtait
+        real(kind(0d0)), intent(in) :: fCpbw, fR, fV, fH, fHdot
+        real(kind(0d0)), intent(in) :: fcgas, fntait, fBtait
 
         real(kind(0d0)) :: tmp1, tmp2, tmp3
         real(kind(0d0)) :: f_rddot_G
@@ -713,7 +711,8 @@ contains
         !!  @param fpb Internal bubble pressure
     function f_cpbw_KM(fR0, fR, fV, fpb)
         !$acc routine seq
-        real(kind(0d0)), intent(IN) :: fR0, fR, fV, fpb
+        real(kind(0d0)), intent(in) :: fR0, fR, fV, fpb
+
         real(kind(0d0)) :: f_cpbw_KM
 
         if (polytropic) then
@@ -740,8 +739,9 @@ contains
         !!  @param fC Current sound speed
     function f_rddot_KM(fpbdot, fCp, fCpbw, fRho, fR, fV, fR0, fC)
         !$acc routine seq
-        real(kind(0d0)), intent(IN) :: fpbdot, fCp, fCpbw
-        real(kind(0d0)), intent(IN) :: fRho, fR, fV, fR0, fC
+        real(kind(0d0)), intent(in) :: fpbdot, fCp, fCpbw
+        real(kind(0d0)), intent(in) :: fRho, fR, fV, fR0, fC
+
         real(kind(0d0)) :: tmp1, tmp2, cdot_star
         real(kind(0d0)) :: f_rddot_KM
 
@@ -771,11 +771,11 @@ contains
 
     !>  Subroutine that computes bubble wall properties for vapor bubbles
         !!  @param pb Internal bubble pressure
-    !!  @param iR0 Current bubble size index
+        !!  @param iR0 Current bubble size index
     subroutine s_bwproperty(pb, iR0)
         !$acc routine seq
-        real(kind(0.d0)), intent(IN) :: pb
-        integer, intent(IN) :: iR0
+        real(kind(0.d0)), intent(in) :: pb
+        integer, intent(in) :: iR0
 
         real(kind(0.d0)) :: x_vw
 
@@ -797,10 +797,10 @@ contains
         !!  @param iR0 Bubble size index
     function f_vflux(fR, fV, fmass_v, iR0)
         !$acc routine seq
-        real(kind(0.d0)), intent(IN) :: fR
-        real(kind(0.d0)), intent(IN) :: fV
-        real(kind(0.d0)), intent(IN) :: fmass_v
-        integer, intent(IN) :: iR0
+        real(kind(0.d0)), intent(in) :: fR
+        real(kind(0.d0)), intent(in) :: fV
+        real(kind(0.d0)), intent(in) :: fmass_v
+        integer, intent(in) :: iR0
 
         real(kind(0.d0)) :: chi_bar
         real(kind(0.d0)) :: grad_chi
@@ -828,12 +828,12 @@ contains
         !!  @param iR0 Bubble size index
     function f_bpres_dot(fvflux, fR, fV, fpb, fmass_v, iR0)
         !$acc routine seq
-        real(kind(0.d0)), intent(IN) :: fvflux
-        real(kind(0.d0)), intent(IN) :: fR
-        real(kind(0.d0)), intent(IN) :: fV
-        real(kind(0.d0)), intent(IN) :: fpb
-        real(kind(0.d0)), intent(IN) :: fmass_v
-        integer, intent(IN) :: iR0
+        real(kind(0.d0)), intent(in) :: fvflux
+        real(kind(0.d0)), intent(in) :: fR
+        real(kind(0.d0)), intent(in) :: fV
+        real(kind(0.d0)), intent(in) :: fpb
+        real(kind(0.d0)), intent(in) :: fmass_v
+        integer, intent(in) :: iR0
 
         real(kind(0.d0)) :: T_bar
         real(kind(0.d0)) :: grad_T
diff --git a/src/simulation/m_cbc.fpp b/src/simulation/m_cbc.fpp
index f1c1bb1131..175de967d8 100644
--- a/src/simulation/m_cbc.fpp
+++ b/src/simulation/m_cbc.fpp
@@ -440,7 +440,7 @@ contains
         !              provided the CBC coordinate direction and location.
 
         ! CBC coordinate direction and location
-        integer, intent(IN) :: cbc_dir_in, cbc_loc_in
+        integer, intent(in) :: cbc_dir_in, cbc_loc_in
 
         ! Cell-boundary locations in the s-direction
         real(kind(0d0)), dimension(0:buff_size + 1) :: s_cb
@@ -544,7 +544,7 @@ contains
     !!  @param cbc_loc_in CBC coordinate location
     subroutine s_associate_cbc_coefficients_pointers(cbc_dir_in, cbc_loc_in)
 
-        integer, intent(IN) :: cbc_dir_in, cbc_loc_in
+        integer, intent(in) :: cbc_dir_in, cbc_loc_in
 
         integer :: i !< Generic loop iterator
 
@@ -622,15 +622,15 @@ contains
 
         type(scalar_field), &
             dimension(sys_size), &
-            intent(IN) :: q_prim_vf
+            intent(in) :: q_prim_vf
 
         type(scalar_field), &
             dimension(sys_size), &
-            intent(INOUT) :: flux_vf, flux_src_vf
+            intent(inout) :: flux_vf, flux_src_vf
 
-        integer, intent(IN) :: cbc_dir_norm, cbc_loc_norm
+        integer, intent(in) :: cbc_dir_norm, cbc_loc_norm
 
-        type(int_bounds_info), intent(IN) :: ix, iy, iz
+        type(int_bounds_info), intent(in) :: ix, iy, iz
 
         ! First-order time derivatives of the partial densities, density,
         ! velocity, pressure, advection variables, and the specific heat
@@ -1024,13 +1024,13 @@ contains
 
         type(scalar_field), &
             dimension(sys_size), &
-            intent(IN) :: q_prim_vf
+            intent(in) :: q_prim_vf
 
         type(scalar_field), &
             dimension(sys_size), &
-            intent(IN) :: flux_vf, flux_src_vf
+            intent(in) :: flux_vf, flux_src_vf
 
-        type(int_bounds_info), intent(IN) :: ix, iy, iz
+        type(int_bounds_info), intent(in) :: ix, iy, iz
 
         integer :: i, j, k, r !< Generic loop iterators
 
@@ -1305,9 +1305,9 @@ contains
 
         type(scalar_field), &
             dimension(sys_size), &
-            intent(INOUT) :: flux_vf, flux_src_vf
+            intent(inout) :: flux_vf, flux_src_vf
 
-        type(int_bounds_info), intent(IN) :: ix, iy, iz
+        type(int_bounds_info), intent(in) :: ix, iy, iz
 
         integer :: i, j, k, r !< Generic loop iterators
 
diff --git a/src/simulation/m_checker.fpp b/src/simulation/m_checker.fpp
index 7ada911842..6426a26463 100644
--- a/src/simulation/m_checker.fpp
+++ b/src/simulation/m_checker.fpp
@@ -19,7 +19,7 @@ module m_checker
 
 contains
 
-    subroutine s_check_inputs()
+    subroutine s_check_inputs
 
         character(len=5) :: iStr
         character(len=5) :: jStr
diff --git a/src/simulation/m_compute_cbc.fpp b/src/simulation/m_compute_cbc.fpp
index ed045104ac..c8386c81a1 100644
--- a/src/simulation/m_compute_cbc.fpp
+++ b/src/simulation/m_compute_cbc.fpp
@@ -33,11 +33,13 @@ contains
 #else
         !$acc routine seq
 #endif
-        real(kind(0d0)), dimension(3), intent(IN) :: lambda
-        real(kind(0d0)), dimension(num_fluids), intent(IN) :: mf, dalpha_rho_ds, dadv_ds
-        real(kind(0d0)), dimension(num_dims), intent(IN) :: dvel_ds
-        real(kind(0d0)), intent(IN) :: rho, c, dpres_ds
-        real(kind(0d0)), dimension(sys_size), intent(INOUT) :: L
+        real(kind(0d0)), dimension(3), intent(in) :: lambda
+        real(kind(0d0)), dimension(sys_size), intent(inout) :: L
+        real(kind(0d0)), intent(in) :: rho, c
+        real(kind(0d0)), dimension(num_fluids), intent(in) :: mf, dalpha_rho_ds
+        real(kind(0d0)), intent(in) :: dpres_ds
+        real(kind(0d0)), dimension(num_dims), intent(in) :: dvel_ds
+        real(kind(0d0)), dimension(num_fluids), intent(in) :: dadv_ds
 
         integer :: i
 
@@ -61,11 +63,13 @@ contains
 #else
         !$acc routine seq
 #endif
-        real(kind(0d0)), dimension(3), intent(IN) :: lambda
-        real(kind(0d0)), dimension(num_fluids), intent(IN) :: mf, dalpha_rho_ds, dadv_ds
-        real(kind(0d0)), dimension(num_dims), intent(IN) :: dvel_ds
-        real(kind(0d0)), intent(IN) :: rho, c, dpres_ds
-        real(kind(0d0)), dimension(sys_size), intent(INOUT) :: L
+        real(kind(0d0)), dimension(3), intent(in) :: lambda
+        real(kind(0d0)), dimension(sys_size), intent(inout) :: L
+        real(kind(0d0)), intent(in) :: rho, c
+        real(kind(0d0)), dimension(num_fluids), intent(in) :: mf, dalpha_rho_ds
+        real(kind(0d0)), intent(in) :: dpres_ds
+        real(kind(0d0)), dimension(num_dims), intent(in) :: dvel_ds
+        real(kind(0d0)), dimension(num_fluids), intent(in) :: dadv_ds
 
         integer :: i !< Generic loop iterator
 
@@ -101,11 +105,13 @@ contains
 #else
         !$acc routine seq
 #endif
-        real(kind(0d0)), dimension(3), intent(IN) :: lambda
-        real(kind(0d0)), dimension(num_fluids), intent(IN) :: mf, dalpha_rho_ds, dadv_ds
-        real(kind(0d0)), dimension(num_dims), intent(IN) :: dvel_ds
-        real(kind(0d0)), intent(IN) :: rho, c, dpres_ds
-        real(kind(0d0)), dimension(sys_size), intent(INOUT) :: L
+        real(kind(0d0)), dimension(3), intent(in) :: lambda
+        real(kind(0d0)), dimension(sys_size), intent(inout) :: L
+        real(kind(0d0)), intent(in) :: rho, c
+        real(kind(0d0)), dimension(num_fluids), intent(in) :: mf, dalpha_rho_ds
+        real(kind(0d0)), intent(in) :: dpres_ds
+        real(kind(0d0)), dimension(num_dims), intent(in) :: dvel_ds
+        real(kind(0d0)), dimension(num_fluids), intent(in) :: dadv_ds
 
         integer :: i
 
@@ -127,11 +133,13 @@ contains
 #else
         !$acc routine seq
 #endif
-        real(kind(0d0)), dimension(3), intent(IN) :: lambda
-        real(kind(0d0)), dimension(num_fluids), intent(IN) :: mf, dalpha_rho_ds, dadv_ds
-        real(kind(0d0)), dimension(num_dims), intent(IN) :: dvel_ds
-        real(kind(0d0)), intent(IN) :: rho, c, dpres_ds
-        real(kind(0d0)), dimension(sys_size), intent(INOUT) :: L
+        real(kind(0d0)), dimension(3), intent(in) :: lambda
+        real(kind(0d0)), dimension(sys_size), intent(inout) :: L
+        real(kind(0d0)), intent(in) :: rho, c
+        real(kind(0d0)), dimension(num_fluids), intent(in) :: mf, dalpha_rho_ds
+        real(kind(0d0)), intent(in) :: dpres_ds
+        real(kind(0d0)), dimension(num_dims), intent(in) :: dvel_ds
+        real(kind(0d0)), dimension(num_fluids), intent(in) :: dadv_ds
 
         integer :: i !> Generic loop iterator
 
@@ -167,11 +175,13 @@ contains
 #else
         !$acc routine seq
 #endif
-        real(kind(0d0)), dimension(3), intent(IN) :: lambda
-        real(kind(0d0)), dimension(num_fluids), intent(IN) :: mf, dalpha_rho_ds, dadv_ds
-        real(kind(0d0)), dimension(num_dims), intent(IN) :: dvel_ds
-        real(kind(0d0)), intent(IN) :: rho, c, dpres_ds
-        real(kind(0d0)), dimension(sys_size), intent(INOUT) :: L
+        real(kind(0d0)), dimension(3), intent(in) :: lambda
+        real(kind(0d0)), dimension(sys_size), intent(inout) :: L
+        real(kind(0d0)), intent(in) :: rho, c
+        real(kind(0d0)), dimension(num_fluids), intent(in) :: mf, dalpha_rho_ds
+        real(kind(0d0)), intent(in) :: dpres_ds
+        real(kind(0d0)), dimension(num_dims), intent(in) :: dvel_ds
+        real(kind(0d0)), dimension(num_fluids), intent(in) :: dadv_ds
 
         integer :: i !> Generic loop iterator
 
@@ -203,11 +213,13 @@ contains
 #else
         !$acc routine seq
 #endif
-        real(kind(0d0)), dimension(3), intent(IN) :: lambda
-        real(kind(0d0)), dimension(num_fluids), intent(IN) :: mf, dalpha_rho_ds, dadv_ds
-        real(kind(0d0)), dimension(num_dims), intent(IN) :: dvel_ds
-        real(kind(0d0)), intent(IN) :: rho, c, dpres_ds
-        real(kind(0d0)), dimension(sys_size), intent(INOUT) :: L
+        real(kind(0d0)), dimension(3), intent(in) :: lambda
+        real(kind(0d0)), dimension(sys_size), intent(inout) :: L
+        real(kind(0d0)), intent(in) :: rho, c
+        real(kind(0d0)), dimension(num_fluids), intent(in) :: mf, dalpha_rho_ds
+        real(kind(0d0)), intent(in) :: dpres_ds
+        real(kind(0d0)), dimension(num_dims), intent(in) :: dvel_ds
+        real(kind(0d0)), dimension(num_fluids), intent(in) :: dadv_ds
 
         integer :: i !> Generic loop iterator
 
@@ -240,12 +252,13 @@ contains
 #else
         !$acc routine seq
 #endif
-        real(kind(0d0)), dimension(3), intent(IN) :: lambda
-        real(kind(0d0)), dimension(num_fluids), intent(IN) :: mf, dalpha_rho_ds, dadv_ds
-        real(kind(0d0)), dimension(num_dims), intent(IN) :: dvel_ds
-        real(kind(0d0)), intent(IN) :: rho, c, dpres_ds
-        real(kind(0d0)), dimension(sys_size), intent(INOUT) :: L
-
+        real(kind(0d0)), dimension(3), intent(in) :: lambda
+        real(kind(0d0)), dimension(sys_size), intent(inout) :: L
+        real(kind(0d0)), intent(in) :: rho, c
+        real(kind(0d0)), dimension(num_fluids), intent(in) :: mf, dalpha_rho_ds
+        real(kind(0d0)), intent(in) :: dpres_ds
+        real(kind(0d0)), dimension(num_dims), intent(in) :: dvel_ds
+        real(kind(0d0)), dimension(num_fluids), intent(in) :: dadv_ds
         integer :: i
 
         do i = 1, advxe
@@ -264,11 +277,13 @@ contains
 #else
         !$acc routine seq
 #endif
-        real(kind(0d0)), dimension(3), intent(IN) :: lambda
-        real(kind(0d0)), dimension(num_fluids), intent(IN) :: mf, dalpha_rho_ds, dadv_ds
-        real(kind(0d0)), dimension(num_dims), intent(IN) :: dvel_ds
-        real(kind(0d0)), intent(IN) :: rho, c, dpres_ds
-        real(kind(0d0)), dimension(sys_size), intent(INOUT) :: L
+        real(kind(0d0)), dimension(3), intent(in) :: lambda
+        real(kind(0d0)), dimension(sys_size), intent(inout) :: L
+        real(kind(0d0)), intent(in) :: rho, c
+        real(kind(0d0)), dimension(num_fluids), intent(in) :: mf, dalpha_rho_ds
+        real(kind(0d0)), intent(in) :: dpres_ds
+        real(kind(0d0)), dimension(num_dims), intent(in) :: dvel_ds
+        real(kind(0d0)), dimension(num_fluids), intent(in) :: dadv_ds
 
         integer :: i !< Generic loop iterator
 
diff --git a/src/simulation/m_compute_levelset.fpp b/src/simulation/m_compute_levelset.fpp
index 4ea3546a81..0af0ca429e 100644
--- a/src/simulation/m_compute_levelset.fpp
+++ b/src/simulation/m_compute_levelset.fpp
@@ -41,9 +41,9 @@ contains
     !>  Initialize IBM module
     subroutine s_compute_circle_levelset(levelset, levelset_norm, ib_patch_id)
 
-        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs), intent(INOUT) :: levelset
-        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs, 3), intent(INOUT) :: levelset_norm
-        integer, intent(IN) :: ib_patch_id
+        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs), intent(inout) :: levelset
+        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs, 3), intent(inout) :: levelset_norm
+        integer, intent(in) :: ib_patch_id
 
         real(kind(0d0)) :: radius, dist
         real(kind(0d0)) :: x_centroid, y_centroid
@@ -77,9 +77,9 @@ contains
 
     subroutine s_compute_airfoil_levelset(levelset, levelset_norm, ib_patch_id)
 
-        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs), intent(INOUT) :: levelset
-        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs, 3), intent(INOUT) :: levelset_norm
-        integer, intent(IN) :: ib_patch_id
+        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs), intent(inout) :: levelset
+        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs, 3), intent(inout) :: levelset_norm
+        integer, intent(in) :: ib_patch_id
 
         real(kind(0d0)) :: radius, dist, global_dist
         integer :: global_id
@@ -160,9 +160,9 @@ contains
 
     subroutine s_compute_3D_airfoil_levelset(levelset, levelset_norm, ib_patch_id)
 
-        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs), intent(INOUT) :: levelset
-        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs, 3), intent(INOUT) :: levelset_norm
-        integer, intent(IN) :: ib_patch_id
+        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs), intent(inout) :: levelset
+        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs, 3), intent(inout) :: levelset_norm
+        integer, intent(in) :: ib_patch_id
 
         real(kind(0d0)) :: radius, dist, dist_surf, dist_side, global_dist
         integer :: global_id
@@ -262,10 +262,10 @@ contains
     !>  Initialize IBM module
     subroutine s_compute_rectangle_levelset(levelset, levelset_norm, ib_patch_id)
 
-        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs), intent(INOUT) :: levelset
-        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs, 3), intent(INOUT) :: levelset_norm
+        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs), intent(inout) :: levelset
+        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs, 3), intent(inout) :: levelset_norm
+        integer, intent(in) :: ib_patch_id
 
-        integer :: ib_patch_id
         real(kind(0d0)) :: top_right(2), bottom_left(2)
         real(kind(0d0)) :: x, y, min_dist
         real(kind(0d0)) :: side_dists(4)
@@ -347,9 +347,9 @@ contains
 
     subroutine s_compute_sphere_levelset(levelset, levelset_norm, ib_patch_id)
 
-        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs), intent(INOUT) :: levelset
-        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs, 3), intent(INOUT) :: levelset_norm
-        integer, intent(IN) :: ib_patch_id
+        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs), intent(inout) :: levelset
+        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs, 3), intent(inout) :: levelset_norm
+        integer, intent(in) :: ib_patch_id
 
         real(kind(0d0)) :: radius, dist
         real(kind(0d0)) :: x_centroid, y_centroid, z_centroid
@@ -384,9 +384,9 @@ contains
 
     subroutine s_compute_cylinder_levelset(levelset, levelset_norm, ib_patch_id)
 
-        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs), intent(INOUT) :: levelset
-        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs, 3), intent(INOUT) :: levelset_norm
-        integer, intent(IN) :: ib_patch_id
+        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs), intent(inout) :: levelset
+        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs, 3), intent(inout) :: levelset_norm
+        integer, intent(in) :: ib_patch_id
 
         real(kind(0d0)) :: radius, dist
         real(kind(0d0)) :: x_centroid, y_centroid, z_centroid
diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp
index 94336b6a70..d11f16d323 100644
--- a/src/simulation/m_data_output.fpp
+++ b/src/simulation/m_data_output.fpp
@@ -34,22 +34,24 @@ module m_data_output
 
     implicit none
 
-    private; public :: s_initialize_data_output_module, &
- s_open_run_time_information_file, &
- s_open_probe_files, &
- s_write_run_time_information, &
- s_write_data_files, &
- s_write_serial_data_files, &
- s_write_parallel_data_files, &
- s_write_probe_files, &
- s_close_run_time_information_file, &
- s_close_probe_files, &
- s_finalize_data_output_module
+    private; 
+    public :: s_initialize_data_output_module, &
+              s_open_run_time_information_file, &
+              s_open_probe_files, &
+              s_write_run_time_information, &
+              s_write_data_files, &
+              s_write_serial_data_files, &
+              s_write_parallel_data_files, &
+              s_write_probe_files, &
+              s_close_run_time_information_file, &
+              s_close_probe_files, &
+              s_finalize_data_output_module
 
     abstract interface ! ===================================================
 
         !> Write data files
         !! @param q_cons_vf Conservative variables
+        !! @param q_prim_vf Primitive variables
         !! @param t_step Current time step
         subroutine s_write_abstract_data_files(q_cons_vf, q_prim_vf, t_step)
 
@@ -57,13 +59,13 @@ module m_data_output
 
             type(scalar_field), &
                 dimension(sys_size), &
-                intent(IN) :: q_cons_vf
+                intent(in) :: q_cons_vf
 
             type(scalar_field), &
                 dimension(sys_size), &
-                intent(INOUT) :: q_prim_vf
+                intent(inout) :: q_prim_vf
 
-            integer, intent(IN) :: t_step
+            integer, intent(in) :: t_step
 
         end subroutine s_write_abstract_data_files
     end interface ! ========================================================
@@ -455,13 +457,13 @@ contains
     !>  The goal of this subroutine is to output the grid and
         !!      conservative variables data files for given time-step.
         !!  @param q_cons_vf Cell-average conservative variables
+        !!  @param q_prim_vf Cell-average primitive variables
         !!  @param t_step Current time-step
     subroutine s_write_serial_data_files(q_cons_vf, q_prim_vf, t_step)
 
-        type(scalar_field), dimension(sys_size), intent(IN) :: q_cons_vf
-        type(scalar_field), dimension(sys_size), intent(INOUT) :: q_prim_vf
-
-        integer, intent(IN) :: t_step
+        type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf
+        type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
+        integer, intent(in) :: t_step
 
         character(LEN=path_len + 2*name_len) :: t_step_dir !<
             !! Relative path to the current time-step directory
@@ -826,18 +828,13 @@ contains
     !>  The goal of this subroutine is to output the grid and
         !!      conservative variables data files for given time-step.
         !!  @param q_cons_vf Cell-average conservative variables
+        !!  @param q_prim_vf Cell-average primitive variables
         !!  @param t_step Current time-step
     subroutine s_write_parallel_data_files(q_cons_vf, q_prim_vf, t_step)
 
-        type(scalar_field), &
-            dimension(sys_size), &
-            intent(IN) :: q_cons_vf
-
-        type(scalar_field), &
-            dimension(sys_size), &
-            intent(INOUT) :: q_prim_vf
-
-        integer, intent(IN) :: t_step
+        type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf
+        type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
+        integer, intent(in) :: t_step
 
 #ifdef MFC_MPI
 
@@ -1014,9 +1011,9 @@ contains
         !!  @param accel_mag Acceleration magnitude information
     subroutine s_write_probe_files(t_step, q_cons_vf, accel_mag)
 
-        integer, intent(IN) :: t_step
-        type(scalar_field), dimension(sys_size), intent(IN) :: q_cons_vf
-        real(kind(0d0)), dimension(0:m, 0:n, 0:p), intent(IN) :: accel_mag
+        integer, intent(in) :: t_step
+        type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf
+        real(kind(0d0)), dimension(0:m, 0:n, 0:p), intent(in) :: accel_mag
 
         real(kind(0d0)), dimension(-1:m) :: distx
         real(kind(0d0)), dimension(-1:n) :: disty
@@ -1659,7 +1656,7 @@ contains
     end subroutine s_close_run_time_information_file
 
     !> Closes probe files
-    subroutine s_close_probe_files()
+    subroutine s_close_probe_files
 
         integer :: i !< Generic loop iterator
 
@@ -1713,7 +1710,7 @@ contains
     end subroutine s_initialize_data_output_module
 
     !> Module deallocation and/or disassociation procedures
-    subroutine s_finalize_data_output_module()
+    subroutine s_finalize_data_output_module
 
         integer :: i !< Generic loop iterator
 
diff --git a/src/simulation/m_derived_variables.f90 b/src/simulation/m_derived_variables.f90
index f27f68284e..4cbe11ab3e 100644
--- a/src/simulation/m_derived_variables.f90
+++ b/src/simulation/m_derived_variables.f90
@@ -114,7 +114,7 @@ end subroutine s_initialize_derived_variables
         !!  @param t_step Current time-step
     subroutine s_compute_derived_variables(t_step)
 
-        integer, intent(IN) :: t_step
+        integer, intent(in) :: t_step
 
         integer :: i, j, k !< Generic loop iterators
 
@@ -178,14 +178,14 @@ end subroutine s_compute_derived_variables
     subroutine s_derive_acceleration_component(i, q_prim_vf0, q_prim_vf1, &
                                                q_prim_vf2, q_prim_vf3, q_sf)
         !DIR$ INLINEALWAYS s_derive_acceleration_component
-        integer, intent(IN) :: i
+        integer, intent(in) :: i
 
-        type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf0
-        type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf1
-        type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf2
-        type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf3
+        type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf0
+        type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf1
+        type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf2
+        type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf3
 
-        real(kind(0d0)), dimension(0:m, 0:n, 0:p), intent(OUT) :: q_sf
+        real(kind(0d0)), dimension(0:m, 0:n, 0:p), intent(out) :: q_sf
 
         integer :: j, k, l, r !< Generic loop iterators
 
diff --git a/src/simulation/m_fftw.fpp b/src/simulation/m_fftw.fpp
index 937ab1dabd..7650e89dce 100644
--- a/src/simulation/m_fftw.fpp
+++ b/src/simulation/m_fftw.fpp
@@ -140,7 +140,7 @@ contains
         !! @param q_cons_vf Conservative variables
     subroutine s_apply_fourier_filter(q_cons_vf)
 
-        type(scalar_field), dimension(sys_size), intent(INOUT) :: q_cons_vf
+        type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
 
         integer :: i, j, k, l !< Generic loop iterators
 
diff --git a/src/simulation/m_hypoelastic.fpp b/src/simulation/m_hypoelastic.fpp
index d01ce4ffaf..e3bb7ec08f 100644
--- a/src/simulation/m_hypoelastic.fpp
+++ b/src/simulation/m_hypoelastic.fpp
@@ -78,9 +78,9 @@ contains
         !!  @param rhs_vf rhs variables
     subroutine s_compute_hypoelastic_rhs(idir, q_prim_vf, rhs_vf)
 
-        type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf
-        type(scalar_field), dimension(sys_size), intent(INOUT) :: rhs_vf
-        integer, intent(IN) :: idir
+        integer, intent(in) :: idir
+        type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
+        type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf
 
         real(kind(0d0)) :: rho_K, G_K
 
diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp
index deef838d60..8b5ed59145 100644
--- a/src/simulation/m_ibm.fpp
+++ b/src/simulation/m_ibm.fpp
@@ -68,7 +68,7 @@ module m_ibm
 contains
 
     !>  Initialize IBM module
-    subroutine s_initialize_ibm_module()
+    subroutine s_initialize_ibm_module
 
         gp_layers = 3
 
@@ -89,7 +89,7 @@ contains
 
     end subroutine s_initialize_ibm_module
 
-    subroutine s_ibm_setup()
+    subroutine s_ibm_setup
 
         integer :: i, j, k
 
@@ -121,18 +121,19 @@ contains
     end subroutine s_ibm_setup
 
     !>  Subroutine that updates the conservative variables at the ghost points
-        !!  @param q_cons_vf
+        !!  @param q_cons_vf Conservative variables
+        !!  @param q_prim_vf Primitive variables
     subroutine s_ibm_correct_state(q_cons_vf, q_prim_vf, pb, mv)
 
         type(scalar_field), &
             dimension(sys_size), &
-            intent(INOUT) :: q_cons_vf !< Primitive Variables
+            intent(inout) :: q_cons_vf !< Conservative Variables
 
         type(scalar_field), &
             dimension(sys_size), &
-            intent(INOUT) :: q_prim_vf !< Primitive Variables
+            intent(inout) :: q_prim_vf !< Primitive Variables
 
-        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), optional, intent(INOUT) :: pb, mv
+        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), optional, intent(inout) :: pb, mv
 
         integer :: i, j, k, l, q, r!< Iterator variables
         integer :: patch_id !< Patch ID of ghost point
@@ -333,16 +334,12 @@ contains
 
     end subroutine s_ibm_correct_state
 
-    !>  Function that computes that bubble wall pressure for Gilmore bubbles
-        !!  @param fR0 Equilibrium bubble radius
-        !!  @param fR Current bubble radius
-        !!  @param fV Current bubble velocity
-        !!  @param fpb Internal bubble pressure
+    !>  Subroutine that computes that bubble wall pressure for Gilmore bubbles
     subroutine s_compute_image_points(ghost_points, levelset, levelset_norm)
 
-        type(ghost_point), dimension(num_gps), intent(INOUT) :: ghost_points
-        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs), intent(IN) :: levelset
-        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs, 3), intent(IN) :: levelset_norm
+        type(ghost_point), dimension(num_gps), intent(inout) :: ghost_points
+        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs), intent(in) :: levelset
+        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs, 3), intent(in) :: levelset_norm
 
         real(kind(0d0)) :: dist
         real(kind(0d0)), dimension(3) :: norm
@@ -445,7 +442,7 @@ contains
 
     end subroutine s_compute_image_points
 
-    subroutine s_find_num_ghost_points()
+    subroutine s_find_num_ghost_points
         integer, dimension(2*gp_layers + 1, 2*gp_layers + 1) &
             :: subsection_2D
         integer, dimension(2*gp_layers + 1, 2*gp_layers + 1, 2*gp_layers + 1) &
@@ -487,8 +484,9 @@ contains
 
     subroutine s_find_ghost_points(ghost_points, inner_points)
 
-        type(ghost_point), dimension(num_gps), intent(INOUT) :: ghost_points
-        type(ghost_point), dimension(num_inner_gps), intent(INOUT) :: inner_points
+        type(ghost_point), dimension(num_gps), intent(inout) :: ghost_points
+        type(ghost_point), dimension(num_inner_gps), intent(inout) :: inner_points
+
         integer, dimension(2*gp_layers + 1, 2*gp_layers + 1) &
             :: subsection_2D
         integer, dimension(2*gp_layers + 1, 2*gp_layers + 1, 2*gp_layers + 1) &
@@ -643,7 +641,7 @@ contains
         !!  @param fpb Internal bubble pressure
     subroutine s_compute_interpolation_coeffs(ghost_points)
 
-        type(ghost_point), dimension(num_gps), intent(INOUT) :: ghost_points
+        type(ghost_point), dimension(num_gps), intent(inout) :: ghost_points
 
         real(kind(0d0)), dimension(2, 2, 2) :: dist
         real(kind(0d0)), dimension(2, 2, 2) :: alpha
@@ -795,18 +793,15 @@ contains
 
     subroutine s_interpolate_image_point(q_prim_vf, gp, alpha_rho_IP, alpha_IP, pres_IP, vel_IP, r_IP, v_IP, pb_IP, mv_IP, nmom_IP, pb, mv, presb_IP, massv_IP)
         !$acc routine seq
-        type(scalar_field), &
-            dimension(sys_size), &
-            intent(IN) :: q_prim_vf !< Primitive Variables
-        real(kind(0d0)), optional, dimension(startx:, starty:, startz:, 1:, 1:), intent(INOUT) :: pb, mv
-
-        type(ghost_point), intent(IN) :: gp
-        real(kind(0d0)), intent(INOUT) :: pres_IP
-        real(kind(0d0)), dimension(3), intent(INOUT) :: vel_IP
-        real(kind(0d0)), dimension(num_fluids), intent(INOUT) :: alpha_IP, alpha_rho_IP
-        real(kind(0d0)), optional, dimension(:), intent(INOUT) :: r_IP, v_IP, pb_IP, mv_IP
-        real(kind(0d0)), optional, dimension(:), intent(INOUT) :: nmom_IP
-        real(kind(0d0)), optional, dimension(:), intent(INOUT) :: presb_IP, massv_IP
+        type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf !< Primitive Variables
+        type(ghost_point), intent(in) :: gp
+        real(kind(0d0)), dimension(num_fluids), intent(inout) :: alpha_IP, alpha_rho_IP
+        real(kind(0d0)), intent(inout) :: pres_IP
+        real(kind(0d0)), dimension(3), intent(inout) :: vel_IP
+        real(kind(0d0)), optional, dimension(:), intent(inout) :: r_IP, v_IP, pb_IP, mv_IP
+        real(kind(0d0)), optional, dimension(:), intent(inout) :: nmom_IP
+        real(kind(0d0)), optional, dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: pb, mv
+        real(kind(0d0)), optional, dimension(:), intent(inout) :: presb_IP, massv_IP
 
         integer :: i, j, k, l, q !< Iterator variables
         integer :: i1, i2, j1, j2, k1, k2 !< Iterator variables
@@ -905,15 +900,11 @@ contains
 
     end subroutine s_interpolate_image_point
 
-    !>  Function that computes that bubble wall pressure for Gilmore bubbles
-        !!  @param fR0 Equilibrium bubble radius
-        !!  @param fR Current bubble radius
-        !!  @param fV Current bubble velocity
-        !!  @param fpb Internal bubble pressure
+    !>  Subroutine that computes that bubble wall pressure for Gilmore bubbles
     subroutine s_compute_levelset(levelset, levelset_norm)
 
-        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs), intent(INOUT) :: levelset
-        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs, 3), intent(INOUT) :: levelset_norm
+        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs), intent(inout) :: levelset
+        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs, 3), intent(inout) :: levelset_norm
         integer :: i !< Iterator variables
         integer :: geometry
 
@@ -936,12 +927,8 @@ contains
 
     end subroutine s_compute_levelset
 
-    !>  Function that computes that bubble wall pressure for Gilmore bubbles
-        !!  @param fR0 Equilibrium bubble radius
-        !!  @param fR Current bubble radius
-        !!  @param fV Current bubble velocity
-        !!  @param fpb Internal bubble pressure
-    subroutine s_finalize_ibm_module()
+    !>  Subroutine that computes that bubble wall pressure for Gilmore bubbles
+    subroutine s_finalize_ibm_module
 
         @:DEALLOCATE(ib_markers%sf)
         @:DEALLOCATE_GLOBAL(levelset)
diff --git a/src/simulation/m_monopole.fpp b/src/simulation/m_monopole.fpp
index 3c316d07e7..d8fe839175 100644
--- a/src/simulation/m_monopole.fpp
+++ b/src/simulation/m_monopole.fpp
@@ -63,7 +63,7 @@ module m_monopole
 
 contains
 
-    subroutine s_initialize_monopole_module()
+    subroutine s_initialize_monopole_module
         integer :: i, j !< generic loop variables
 
         @:ALLOCATE_GLOBAL(mag(1:num_mono), support(1:num_mono), length(1:num_mono), npulse(1:num_mono), pulse(1:num_mono), dir(1:num_mono), delay(1:num_mono), loc_mono(1:3, 1:num_mono), foc_length(1:num_mono), aperture(1:num_mono), support_width(1:num_mono))
@@ -91,7 +91,7 @@ contains
 
     end subroutine
 
-    subroutine s_compute_monopole_rhs()
+    subroutine s_compute_monopole_rhs
 
     end subroutine s_compute_monopole_rhs
 
@@ -108,9 +108,9 @@ contains
         !! are calculated from the conservative variables and gradient magnitude (GM)
         !! of the volume fractions, q_cons_qp and gm_alpha_qp, respectively.
 
-        type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf
+        integer, intent(in) :: t_step, id
 
-        integer, intent(IN) :: t_step, id
+        type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf
 
         real(kind(0d0)) :: myR, myV, alf, myP, myRho, R2Vav
 
@@ -297,12 +297,13 @@ contains
         !! @param mysos Alternative speed of sound for testing
     function f_g(the_time, sos, mysos, nm, term_index)
         !$acc routine seq
-        real(kind(0d0)), intent(IN) :: the_time, sos, mysos
-        integer, intent(IN) :: nm
+        real(kind(0d0)), intent(in) :: the_time, sos, mysos
+        integer, intent(in) :: nm
+        integer, intent(in) :: term_index
+
         real(kind(0d0)) :: period, t0, sigt, pa
         real(kind(0d0)) :: offset
         real(kind(0d0)) :: f_g
-        integer :: term_index
 
         offset = 0d0
         if (delay(nm) /= dflt_real) offset = delay(nm)
@@ -345,10 +346,12 @@ contains
     function f_delta(j, k, l, mono_loc, mono_leng, nm, angle, angle_z)
 
         !$acc routine seq
-        real(kind(0d0)), dimension(3), intent(IN) :: mono_loc
-        integer, intent(IN) :: nm
-        real(kind(0d0)), intent(IN) :: mono_leng
         integer, intent(in) :: j, k, l
+        real(kind(0d0)), dimension(3), intent(in) :: mono_loc
+        real(kind(0d0)), intent(in) :: mono_leng
+        integer, intent(in) :: nm
+        real(kind(0d0)), intent(out) :: angle
+        real(kind(0d0)), intent(out) :: angle_z
 
         integer :: q
         real(kind(0d0)) :: h, hx, hy, hz
@@ -357,8 +360,6 @@ contains
         real(kind(0d0)) :: hxnew_cyl, hynew_cyl
         real(kind(0d0)) :: sig
         real(kind(0d0)) :: f_delta
-        real(kind(0d0)) :: angle
-        real(kind(0d0)) :: angle_z
 
         if (n == 0) then
             sig = dx(j)
diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp
index 94272f15e1..8b5b9d2c1f 100644
--- a/src/simulation/m_mpi_proxy.fpp
+++ b/src/simulation/m_mpi_proxy.fpp
@@ -648,8 +648,8 @@ contains
         !!  @param pbc_loc Processor boundary condition (PBC) location
     subroutine s_mpi_sendrecv_grid_variables_buffers(mpi_dir, pbc_loc)
 
-        integer, intent(IN) :: mpi_dir
-        integer, intent(IN) :: pbc_loc
+        integer, intent(in) :: mpi_dir
+        integer, intent(in) :: pbc_loc
 
         integer :: dst_proc(1:3)
 
@@ -834,10 +834,9 @@ contains
                                                 mpi_dir, &
                                                 pbc_loc)
 
-        type(scalar_field), dimension(sys_size), intent(INOUT) :: q_cons_vf
-        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(INOUT) :: pb, mv
-
-        integer, intent(IN) :: mpi_dir, pbc_loc
+        type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
+        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: pb, mv
+        integer, intent(in) :: mpi_dir, pbc_loc
 
         integer :: i, j, k, l, r, q !< Generic loop iterators
 
@@ -1247,14 +1246,10 @@ contains
     !>  The goal of this procedure is to populate the buffers of
         !!      the cell-average conservative variables by communicating
         !!      with the neighboring processors.
-        !!  @param q_cons_vf Cell-average conservative variables
-        !!  @param mpi_dir MPI communication coordinate direction
-        !!  @param pbc_loc Processor boundary condition (PBC) location
     subroutine s_mpi_sendrecv_ib_buffers(ib_markers, gp_layers)
 
-        type(integer_field), intent(INOUT) :: ib_markers
-
-        integer, intent(IN) :: gp_layers
+        type(integer_field), intent(inout) :: ib_markers
+        integer, intent(in) :: gp_layers
 
         integer :: i, j, k, l, r !< Generic loop iterators
 
@@ -2038,9 +2033,8 @@ contains
 
     subroutine s_mpi_sendrecv_capilary_variables_buffers(c_divs_vf, mpi_dir, pbc_loc)
 
-        type(scalar_field), dimension(num_dims + 1), intent(INOUT) :: c_divs_vf
-
-        integer, intent(IN) :: mpi_dir, pbc_loc
+        type(scalar_field), dimension(num_dims + 1), intent(inout) :: c_divs_vf
+        integer, intent(in) :: mpi_dir, pbc_loc
 
         integer :: i, j, k, l, r, q !< Generic loop iterators
 
diff --git a/src/simulation/m_qbmm.fpp b/src/simulation/m_qbmm.fpp
index 526ca03dd2..c44785da54 100644
--- a/src/simulation/m_qbmm.fpp
+++ b/src/simulation/m_qbmm.fpp
@@ -55,7 +55,7 @@ module m_qbmm
 
 contains
 
-    subroutine s_initialize_qbmm_module()
+    subroutine s_initialize_qbmm_module
 
         integer :: i1, i2, q, i, j
 
@@ -425,11 +425,13 @@ contains
 
     subroutine s_compute_qbmm_rhs(idir, q_cons_vf, q_prim_vf, rhs_vf, flux_n_vf, pb, rhs_pb, mv, rhs_mv)
 
-        type(scalar_field), dimension(sys_size) :: q_cons_vf, q_prim_vf, rhs_vf, flux_n_vf
-        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(INOUT) :: pb, mv
-        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(INOUT) :: rhs_pb, rhs_mv
+        integer, intent(in) :: idir
+        type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf, q_prim_vf
+        type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf
+        type(scalar_field), dimension(sys_size), intent(in) :: flux_n_vf
+        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: pb, rhs_pb
+        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: mv, rhs_mv
 
-        integer :: idir
         integer :: i, j, k, l, q
 
         real(kind(0d0)) :: nb_q, nb_dot, R, R2, nR, nR2, nR_dot, nR2_dot, var, AX
@@ -695,8 +697,9 @@ contains
 #else
         !$acc routine seq
 #endif
-        real(kind(0.d0)), intent(IN) :: pres, rho, c
-        real(kind(0.d0)), dimension(nterms, 0:2, 0:2), intent(OUT) :: coeffs
+        real(kind(0.d0)), intent(in) :: pres, rho, c
+        real(kind(0.d0)), dimension(nterms, 0:2, 0:2), intent(out) :: coeffs
+
         integer :: i1, i2, q
 
         coeffs = 0d0
@@ -768,8 +771,9 @@ contains
         !$acc routine seq
 #endif
 
-        real(kind(0.d0)), intent(INOUT) :: pres, rho, c
-        real(kind(0.d0)), dimension(nterms, 0:2, 0:2), intent(OUT) :: coeffs
+        real(kind(0.d0)), intent(inout) :: pres, rho, c
+        real(kind(0.d0)), dimension(nterms, 0:2, 0:2), intent(out) :: coeffs
+
         integer :: i1, i2, q
 
         coeffs = 0d0
@@ -825,13 +829,13 @@ contains
 
     subroutine s_mom_inv(q_cons_vf, q_prim_vf, momsp, moms3d, pb, rhs_pb, mv, rhs_mv, ix, iy, iz, nbub_sc)
 
-        type(scalar_field), dimension(:), intent(INOUT) :: q_prim_vf, q_cons_vf
-        type(scalar_field), dimension(:), intent(INOUT) :: momsp
-        type(scalar_field), dimension(0:, 0:, :), intent(INOUT) :: moms3d
-        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(INOUT) :: pb, mv
-        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(INOUT) :: rhs_pb, rhs_mv
-        real(kind(0d0)), dimension(startx:, starty:, startz:) :: nbub_sc
-        type(int_bounds_info), intent(IN) :: ix, iy, iz
+        type(scalar_field), dimension(:), intent(inout) :: q_cons_vf, q_prim_vf
+        type(scalar_field), dimension(:), intent(inout) :: momsp
+        type(scalar_field), dimension(0:, 0:, :), intent(inout) :: moms3d
+        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: pb, rhs_pb
+        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: mv, rhs_mv
+        type(int_bounds_info), intent(in) :: ix, iy, iz
+        real(kind(0d0)), dimension(startx:, starty:, startz:) :: nbub_sc !> Unused Variable not sure what to put as intent
 
         real(kind(0d0)), dimension(nmom) :: moms, msum
         real(kind(0d0)), dimension(nnode, nb) :: wght, abscX, abscY, wght_pb, wght_mv, wght_ht, ht
@@ -1045,8 +1049,8 @@ contains
 #else
         !$acc routine seq
 #endif
-        real(kind(0d0)), dimension(nnode), intent(INOUT) :: wght, abscX, abscY
-        real(kind(0d0)), dimension(nmom), intent(IN) :: momin
+        real(kind(0d0)), dimension(nmom), intent(in) :: momin
+        real(kind(0d0)), dimension(nnode), intent(inout) :: wght, abscX, abscY
 
         real(kind(0d0)), dimension(0:2, 0:2) :: moms
         real(kind(0d0)), dimension(3) :: M1, M3
@@ -1112,8 +1116,9 @@ contains
 #else
         !$acc routine seq
 #endif
-        real(kind(0d0)), dimension(2), intent(INOUT) :: frho, fup
-        real(kind(0d0)), dimension(3), intent(IN) :: fmom
+        real(kind(0d0)), dimension(2), intent(inout) :: frho, fup
+        real(kind(0d0)), dimension(3), intent(in) :: fmom
+
         real(kind(0d0)) :: bu, d2, c2
 
         bu = fmom(2)/fmom(1)
@@ -1129,8 +1134,9 @@ contains
 
     function f_quad(abscX, abscY, wght_in, q, r, s)
         !$acc routine seq
-        real(kind(0.d0)), dimension(nnode, nb), intent(IN) :: abscX, abscY, wght_in
-        real(kind(0.d0)), intent(IN) :: q, r, s
+        real(kind(0.d0)), dimension(nnode, nb), intent(in) :: abscX, abscY, wght_in
+        real(kind(0.d0)), intent(in) :: q, r, s
+
         real(kind(0.d0)) :: f_quad_RV, f_quad
         integer :: i
 
@@ -1144,9 +1150,9 @@ contains
 
     function f_quad2D(abscX, abscY, wght_in, pow)
         !$acc routine seq
-        real(kind(0.d0)), dimension(nnode), intent(IN) :: abscX, abscY, wght_in
+        real(kind(0.d0)), dimension(nnode), intent(in) :: abscX, abscY, wght_in
+        real(kind(0.d0)), dimension(3), intent(in) :: pow
 
-        real(kind(0.d0)), dimension(3), intent(IN) :: pow
         real(kind(0.d0)) :: f_quad2D
 
         f_quad2D = sum(wght_in(:)*(abscX(:)**pow(1))*(abscY(:)**pow(2)))
diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp
index 8f451c3cfa..cb7e30255c 100644
--- a/src/simulation/m_rhs.fpp
+++ b/src/simulation/m_rhs.fpp
@@ -698,15 +698,16 @@ contains
 
     subroutine s_compute_rhs(q_cons_vf, q_prim_vf, rhs_vf, pb, rhs_pb, mv, rhs_mv, t_step, time_avg)
 
-        type(scalar_field), dimension(sys_size), intent(INOUT) :: q_cons_vf
-        type(scalar_field), dimension(sys_size), intent(INOUT) :: q_prim_vf
-        type(scalar_field), dimension(sys_size), intent(INOUT) :: rhs_vf
-        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(INOUT) :: pb, mv
-        real(kind(0d0)), intent(INOUT) :: time_avg
+        type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
+        type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
+        type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf
+        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: pb, rhs_pb
+        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: mv, rhs_mv
+        integer, intent(in) :: t_step
+        real(kind(0d0)), intent(inout) :: time_avg
+
         real(kind(0d0)) :: t_start, t_finish
         real(kind(0d0)) :: gp_sum
-        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(INOUT) :: rhs_pb, rhs_mv
-        integer, intent(IN) :: t_step
 
         real(kind(0d0)) :: top, bottom  !< Numerator and denominator when evaluating flux limiter function
         real(kind(0d0)), dimension(num_fluids) :: myalpha_rho, myalpha
@@ -1041,12 +1042,12 @@ contains
 
     subroutine s_compute_advection_source_term(idir, rhs_vf, q_cons_vf, q_prim_vf, flux_src_n_vf)
 
-        type(vector_field), intent(INOUT) :: q_cons_vf
-        type(vector_field), intent(INOUT) :: q_prim_vf
-        type(scalar_field), dimension(sys_size), intent(INOUT) :: rhs_vf
-        type(vector_field), intent(INOUT) :: flux_src_n_vf
-
         integer, intent(in) :: idir
+        type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf
+        type(vector_field), intent(inout) :: q_cons_vf
+        type(vector_field), intent(inout) :: q_prim_vf
+        type(vector_field), intent(inout) :: flux_src_n_vf
+
         integer :: i, j, k, l, q
 
         if (alt_soundspeed) then
@@ -1597,14 +1598,13 @@ contains
     subroutine s_compute_additional_physics_rhs(idir, q_prim_vf, rhs_vf, flux_src_n, &
                                                 dq_prim_dx_vf, dq_prim_dy_vf, dq_prim_dz_vf, ixt, iyt, izt)
 
-        type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf, &
-                                                               flux_src_n, &
-                                                               dq_prim_dx_vf, &
-                                                               dq_prim_dy_vf, &
-                                                               dq_prim_dz_vf
-        type(scalar_field), dimension(sys_size), intent(INOUT) :: rhs_vf
-        type(int_bounds_info) :: ixt, iyt, izt
-        integer, intent(IN) :: idir
+        integer, intent(in) :: idir
+        type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
+        type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf
+        type(scalar_field), dimension(sys_size), intent(in) :: flux_src_n
+        type(scalar_field), dimension(sys_size), intent(in) :: dq_prim_dx_vf, dq_prim_dy_vf, dq_prim_dz_vf
+        type(int_bounds_info), intent(in) :: ixt, iyt, izt
+
         integer :: i, j, k, l, q
 
         if (idir == 1) then ! x-direction
@@ -1837,7 +1837,7 @@ contains
         !!  @param q_cons_vf Cell-average conservative variables
     subroutine s_pressure_relaxation_procedure(q_cons_vf)
 
-        type(scalar_field), dimension(sys_size), intent(INOUT) :: q_cons_vf
+        type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
 
         !> @name Relaxed pressure, initial partial pressures, function f(p) and its partial
             !! derivative df(p), isentropic partial density, sum of volume fractions,
@@ -2094,11 +2094,10 @@ contains
     subroutine s_reconstruct_cell_boundary_values(v_vf, vL_x, vL_y, vL_z, vR_x, vR_y, vR_z, &
                                                   norm_dir)
 
-        type(scalar_field), dimension(iv%beg:iv%end), intent(IN) :: v_vf
-
-        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(INOUT) :: vL_x, vL_y, vL_z, vR_x, vR_y, vR_z
-
-        integer, intent(IN) :: norm_dir
+        type(scalar_field), dimension(iv%beg:iv%end), intent(in) :: v_vf
+        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(inout) :: vL_x, vL_y, vL_z
+        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(inout) :: vR_x, vR_y, vR_z
+        integer, intent(in) :: norm_dir
 
         integer :: weno_dir !< Coordinate direction of the WENO reconstruction
 
@@ -2149,11 +2148,10 @@ contains
     subroutine s_reconstruct_cell_boundary_values_first_order(v_vf, vL_x, vL_y, vL_z, vR_x, vR_y, vR_z, &
                                                               norm_dir)
 
-        type(scalar_field), dimension(iv%beg:iv%end), intent(IN) :: v_vf
-
-        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(INOUT) :: vL_x, vL_y, vL_z, vR_x, vR_y, vR_z
-
-        integer, intent(IN) :: norm_dir
+        type(scalar_field), dimension(iv%beg:iv%end), intent(in) :: v_vf
+        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(inout) :: vL_x, vL_y, vL_z
+        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(inout) :: vR_x, vR_y, vR_z
+        integer, intent(in) :: norm_dir
 
         integer :: recon_dir !< Coordinate direction of the WENO reconstruction
 
diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp
index 58489cd47b..492a7de8e0 100644
--- a/src/simulation/m_riemann_solvers.fpp
+++ b/src/simulation/m_riemann_solvers.fpp
@@ -94,24 +94,24 @@ module m_riemann_solvers
 
             import :: scalar_field, int_bounds_info, sys_size, startx, starty, startz
 
-            real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(INOUT) :: qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, qR_prim_rsx_vf, qR_prim_rsy_vf, qR_prim_rsz_vf
-            type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf
+            real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(inout) :: qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, qR_prim_rsx_vf, qR_prim_rsy_vf, qR_prim_rsz_vf
+            type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
 
-            type(scalar_field), allocatable, dimension(:), intent(INOUT) :: qL_prim_vf, qR_prim_vf
+            type(scalar_field), allocatable, dimension(:), intent(inout) :: qL_prim_vf, qR_prim_vf
 
             type(scalar_field), &
                 allocatable, dimension(:), &
-                intent(INOUT) :: dqL_prim_dx_vf, dqR_prim_dx_vf, &
+                intent(inout) :: dqL_prim_dx_vf, dqR_prim_dx_vf, &
                                  dqL_prim_dy_vf, dqR_prim_dy_vf, &
                                  dqL_prim_dz_vf, dqR_prim_dz_vf
 
             type(scalar_field), &
                 dimension(sys_size), &
-                intent(INOUT) :: flux_vf, flux_src_vf, flux_gsrc_vf
+                intent(inout) :: flux_vf, flux_src_vf, flux_gsrc_vf
 
-            integer, intent(IN) :: norm_dir
+            integer, intent(in) :: norm_dir
 
-            type(int_bounds_info), intent(IN) :: ix, iy, iz
+            type(int_bounds_info), intent(in) :: ix, iy, iz
 
         end subroutine s_abstract_riemann_solver
 
@@ -136,18 +136,18 @@ module m_riemann_solvers
 
             type(scalar_field), &
                 dimension(num_dims), &
-                intent(IN) :: velL_vf, velR_vf, &
+                intent(in) :: velL_vf, velR_vf, &
                               dvelL_dx_vf, dvelR_dx_vf, &
                               dvelL_dy_vf, dvelR_dy_vf, &
                               dvelL_dz_vf, dvelR_dz_vf
 
             type(scalar_field), &
                 dimension(sys_size), &
-                intent(INOUT) :: flux_src_vf
+                intent(inout) :: flux_src_vf
 
-            integer, intent(IN) :: norm_dir
+            integer, intent(in) :: norm_dir
 
-            type(int_bounds_info), intent(IN) :: ix, iy, iz
+            type(int_bounds_info), intent(in) :: ix, iy, iz
 
         end subroutine s_compute_abstract_viscous_source_flux
 
@@ -279,24 +279,24 @@ contains
                                     flux_gsrc_vf, &
                                     norm_dir, ix, iy, iz)
 
-        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(INOUT) :: qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, qR_prim_rsx_vf, qR_prim_rsy_vf, qR_prim_rsz_vf
-        type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf
+        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(inout) :: qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, qR_prim_rsx_vf, qR_prim_rsy_vf, qR_prim_rsz_vf
+        type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
 
-        type(scalar_field), allocatable, dimension(:), intent(INOUT) :: qL_prim_vf, qR_prim_vf
+        type(scalar_field), allocatable, dimension(:), intent(inout) :: qL_prim_vf, qR_prim_vf
 
         type(scalar_field), &
             allocatable, dimension(:), &
-            intent(INOUT) :: dqL_prim_dx_vf, dqR_prim_dx_vf, &
+            intent(inout) :: dqL_prim_dx_vf, dqR_prim_dx_vf, &
                              dqL_prim_dy_vf, dqR_prim_dy_vf, &
                              dqL_prim_dz_vf, dqR_prim_dz_vf
 
         ! Intercell fluxes
         type(scalar_field), &
             dimension(sys_size), &
-            intent(INOUT) :: flux_vf, flux_src_vf, flux_gsrc_vf
+            intent(inout) :: flux_vf, flux_src_vf, flux_gsrc_vf
 
-        integer, intent(IN) :: norm_dir
-        type(int_bounds_info), intent(IN) :: ix, iy, iz
+        integer, intent(in) :: norm_dir
+        type(int_bounds_info), intent(in) :: ix, iy, iz
 
         real(kind(0d0)), dimension(num_fluids) :: alpha_rho_L, alpha_rho_R
         real(kind(0d0)) :: rho_L, rho_R
@@ -820,24 +820,24 @@ contains
                                      flux_gsrc_vf, &
                                      norm_dir, ix, iy, iz)
 
-        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(INOUT) :: qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, qR_prim_rsx_vf, qR_prim_rsy_vf, qR_prim_rsz_vf
-        type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf
+        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(inout) :: qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, qR_prim_rsx_vf, qR_prim_rsy_vf, qR_prim_rsz_vf
+        type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
 
-        type(scalar_field), allocatable, dimension(:), intent(INOUT) :: qL_prim_vf, qR_prim_vf
+        type(scalar_field), allocatable, dimension(:), intent(inout) :: qL_prim_vf, qR_prim_vf
 
         type(scalar_field), &
             allocatable, dimension(:), &
-            intent(INOUT) :: dqL_prim_dx_vf, dqR_prim_dx_vf, &
+            intent(inout) :: dqL_prim_dx_vf, dqR_prim_dx_vf, &
                              dqL_prim_dy_vf, dqR_prim_dy_vf, &
                              dqL_prim_dz_vf, dqR_prim_dz_vf
 
         ! Intercell fluxes
         type(scalar_field), &
             dimension(sys_size), &
-            intent(INOUT) :: flux_vf, flux_src_vf, flux_gsrc_vf
+            intent(inout) :: flux_vf, flux_src_vf, flux_gsrc_vf
 
-        integer, intent(IN) :: norm_dir
-        type(int_bounds_info), intent(IN) :: ix, iy, iz
+        integer, intent(in) :: norm_dir
+        type(int_bounds_info), intent(in) :: ix, iy, iz
 
         real(kind(0d0)), dimension(num_fluids) :: alpha_rho_L, alpha_rho_R
         real(kind(0d0)) :: rho_L, rho_R
@@ -2492,18 +2492,17 @@ contains
         qR_prim_vf, &
         norm_dir, ix, iy, iz)
 
-        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(INOUT) :: qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, qR_prim_rsx_vf, qR_prim_rsy_vf, qR_prim_rsz_vf
+        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(inout) :: qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, qR_prim_rsx_vf, qR_prim_rsy_vf, qR_prim_rsz_vf
 
         type(scalar_field), &
             allocatable, dimension(:), &
-            intent(INOUT) :: dqL_prim_dx_vf, dqR_prim_dx_vf, &
+            intent(inout) :: dqL_prim_dx_vf, dqR_prim_dx_vf, &
                              dqL_prim_dy_vf, dqR_prim_dy_vf, &
                              dqL_prim_dz_vf, dqR_prim_dz_vf, &
                              qL_prim_vf, qR_prim_vf
 
-        integer, intent(IN) :: norm_dir
-
-        type(int_bounds_info), intent(IN) :: ix, iy, iz
+        integer, intent(in) :: norm_dir
+        type(int_bounds_info), intent(in) :: ix, iy, iz
 
         integer :: i, j, k, l !< Generic loop iterator
 
@@ -2870,15 +2869,13 @@ contains
         flux_gsrc_vf, &
         norm_dir, ix, iy, iz)
 
-        type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf
-
+        type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
         type(scalar_field), &
             dimension(sys_size), &
-            intent(INOUT) :: flux_vf, flux_src_vf, flux_gsrc_vf
-
-        integer, intent(IN) :: norm_dir
+            intent(inout) :: flux_vf, flux_src_vf, flux_gsrc_vf
 
-        type(int_bounds_info), intent(IN) :: ix, iy, iz
+        integer, intent(in) :: norm_dir
+        type(int_bounds_info), intent(in) :: ix, iy, iz
 
         integer :: i, j, k, l ! Generic loop iterators
 
@@ -3012,18 +3009,18 @@ contains
 
         type(scalar_field), &
             dimension(num_dims), &
-            intent(IN) :: velL_vf, velR_vf, &
+            intent(in) :: velL_vf, velR_vf, &
                           dvelL_dx_vf, dvelR_dx_vf, &
                           dvelL_dy_vf, dvelR_dy_vf, &
                           dvelL_dz_vf, dvelR_dz_vf
 
         type(scalar_field), &
             dimension(sys_size), &
-            intent(INOUT) :: flux_src_vf
+            intent(inout) :: flux_src_vf
 
-        integer, intent(IN) :: norm_dir
+        integer, intent(in) :: norm_dir
 
-        type(int_bounds_info), intent(IN) :: ix, iy, iz
+        type(int_bounds_info), intent(in) :: ix, iy, iz
 
         ! Arithmetic mean of the left and right, WENO-reconstructed, cell-
         ! boundary values of cell-average first-order spatial derivatives
@@ -3539,18 +3536,17 @@ contains
 
         type(scalar_field), &
             dimension(num_dims), &
-            intent(IN) :: velL_vf, velR_vf, &
+            intent(in) :: velL_vf, velR_vf, &
                           dvelL_dx_vf, dvelR_dx_vf, &
                           dvelL_dy_vf, dvelR_dy_vf, &
                           dvelL_dz_vf, dvelR_dz_vf
 
         type(scalar_field), &
             dimension(sys_size), &
-            intent(INOUT) :: flux_src_vf
+            intent(inout) :: flux_src_vf
 
-        integer, intent(IN) :: norm_dir
-
-        type(int_bounds_info), intent(IN) :: ix, iy, iz
+        integer, intent(in) :: norm_dir
+        type(int_bounds_info), intent(in) :: ix, iy, iz
 
         ! Arithmetic mean of the left and right, WENO-reconstructed, cell-
         ! boundary values of cell-average first-order spatial derivatives
@@ -4017,11 +4013,10 @@ contains
 
         type(scalar_field), &
             dimension(sys_size), &
-            intent(INOUT) :: flux_vf, flux_src_vf, flux_gsrc_vf
-
-        integer, intent(IN) :: norm_dir
+            intent(inout) :: flux_vf, flux_src_vf, flux_gsrc_vf
 
-        type(int_bounds_info), intent(IN) :: ix, iy, iz
+        integer, intent(in) :: norm_dir
+        type(int_bounds_info), intent(in) :: ix, iy, iz
 
         integer :: i, j, k, l !< Generic loop iterators
 
diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp
index e8911ad290..9bb39173f6 100644
--- a/src/simulation/m_start_up.fpp
+++ b/src/simulation/m_start_up.fpp
@@ -97,7 +97,7 @@ module m_start_up
 
             type(scalar_field), &
                 dimension(sys_size), &
-                intent(INOUT) :: q_cons_vf
+                intent(inout) :: q_cons_vf
 
         end subroutine s_read_abstract_data_files
 
@@ -1030,7 +1030,8 @@ contains
         !! @param v_vf conservative variables
     subroutine s_initialize_internal_energy_equations(v_vf)
 
-        type(scalar_field), dimension(sys_size), intent(INOUT) :: v_vf
+        type(scalar_field), dimension(sys_size), intent(inout) :: v_vf
+
         real(kind(0d0)) :: rho
         real(kind(0d0)) :: dyn_pres
         real(kind(0d0)) :: gamma
@@ -1066,17 +1067,17 @@ contains
             end do
         end do
 
-    end subroutine s_initialize_internal_energy_equations !-----------------
+    end subroutine s_initialize_internal_energy_equations
 
     subroutine s_perform_time_step(t_step, time_avg, time_final, io_time_avg, io_time_final, proc_time, io_proc_time, file_exists, start, finish, nt)
-        integer, intent(INOUT) :: t_step
-        real(kind(0d0)), intent(INOUT) :: time_avg, time_final
-        real(kind(0d0)), intent(INOUT) :: io_time_avg, io_time_final
-        real(kind(0d0)), dimension(:), intent(INOUT) :: proc_time
-        real(kind(0d0)), dimension(:), intent(INOUT) :: io_proc_time
-        logical, intent(INOUT) :: file_exists
-        real(kind(0d0)), intent(INOUT) :: start, finish
-        integer, intent(INOUT) :: nt
+        integer, intent(inout) :: t_step
+        real(kind(0d0)), intent(inout) :: time_avg, time_final
+        real(kind(0d0)), intent(inout) :: io_time_avg, io_time_final
+        real(kind(0d0)), dimension(:), intent(inout) :: proc_time
+        real(kind(0d0)), dimension(:), intent(inout) :: io_proc_time
+        logical, intent(inout) :: file_exists
+        real(kind(0d0)), intent(inout) :: start, finish
+        integer, intent(inout) :: nt
 
         integer :: i, j, k, l
 
@@ -1120,14 +1121,14 @@ contains
 
     subroutine s_save_performance_metrics(t_step, time_avg, time_final, io_time_avg, io_time_final, proc_time, io_proc_time, file_exists, start, finish, nt)
 
-        integer, intent(INOUT) :: t_step
-        real(kind(0d0)), intent(INOUT) :: time_avg, time_final
-        real(kind(0d0)), intent(INOUT) :: io_time_avg, io_time_final
-        real(kind(0d0)), dimension(:), intent(INOUT) :: proc_time
-        real(kind(0d0)), dimension(:), intent(INOUT) :: io_proc_time
-        logical, intent(INOUT) :: file_exists
-        real(kind(0d0)), intent(INOUT) :: start, finish
-        integer, intent(INOUT) :: nt
+        integer, intent(inout) :: t_step
+        real(kind(0d0)), intent(inout) :: time_avg, time_final
+        real(kind(0d0)), intent(inout) :: io_time_avg, io_time_final
+        real(kind(0d0)), dimension(:), intent(inout) :: proc_time
+        real(kind(0d0)), dimension(:), intent(inout) :: io_proc_time
+        logical, intent(inout) :: file_exists
+        real(kind(0d0)), intent(inout) :: start, finish
+        integer, intent(inout) :: nt
 
         call s_mpi_barrier()
 
@@ -1175,8 +1176,10 @@ contains
     end subroutine s_save_performance_metrics
 
     subroutine s_save_data(t_step, start, finish, io_time_avg, nt)
-        real(kind(0d0)), intent(INOUT) :: start, finish, io_time_avg
-        integer, intent(INOUT) :: t_step, nt
+        integer, intent(inout) :: t_step
+        real(kind(0d0)), intent(inout) :: start, finish, io_time_avg
+        integer, intent(inout) :: nt
+        
         integer :: i, j, k, l
 
         if (mod(t_step - t_step_start, t_step_save) == 0 .or. t_step == t_step_stop) then
@@ -1215,7 +1218,7 @@ contains
 
     end subroutine s_save_data
 
-    subroutine s_initialize_modules()
+    subroutine s_initialize_modules
         call s_initialize_global_parameters_module()
         !Quadrature weights and nodes for polydisperse simulations
         if (bubbles .and. nb > 1 .and. R0_type == 1) then
@@ -1310,7 +1313,7 @@ contains
 
     end subroutine s_initialize_modules
 
-    subroutine s_initialize_mpi_domain()
+    subroutine s_initialize_mpi_domain
         integer :: ierr
 #ifdef MFC_OpenACC
         real(kind(0d0)) :: starttime, endtime
@@ -1382,7 +1385,7 @@ contains
 
     end subroutine s_initialize_mpi_domain
 
-    subroutine s_initialize_gpu_vars()
+    subroutine s_initialize_gpu_vars
         integer :: i
         !Update GPU DATA
         do i = 1, sys_size
@@ -1411,7 +1414,7 @@ contains
 
     end subroutine s_initialize_gpu_vars
 
-    subroutine s_finalize_modules()
+    subroutine s_finalize_modules
         ! Disassociate pointers for serial and parallel I/O
         s_read_data_files => null()
         s_write_data_files => null()
diff --git a/src/simulation/m_surface_tension.fpp b/src/simulation/m_surface_tension.fpp
index de1d567942..2e9ca8fedf 100644
--- a/src/simulation/m_surface_tension.fpp
+++ b/src/simulation/m_surface_tension.fpp
@@ -57,7 +57,7 @@ module m_surface_tension
 
 contains
 
-    subroutine s_initialize_surface_tension_module()
+    subroutine s_initialize_surface_tension_module
 
         ! Configuring Coordinate Direction Indexes =========================
         ix%beg = -buff_size; iy%beg = 0; iz%beg = 0
@@ -91,15 +91,15 @@ contains
                                               flux_src_vf, &
                                               id, isx, isy, isz)
 
-        type(int_bounds_info) :: isx, isy, isz
-        type(scalar_field), dimension(sys_size) :: q_prim_vf
-        real(kind(0d0)), dimension(-1:, 0:, 0:, 1:), intent(IN) :: vSrc_rsx_vf
-        real(kind(0d0)), dimension(-1:, 0:, 0:, 1:), intent(IN) :: vSrc_rsy_vf
-        real(kind(0d0)), dimension(-1:, 0:, 0:, 1:), intent(IN) :: vSrc_rsz_vf
+        type(scalar_field), dimension(sys_size) :: q_prim_vf !> unused so unsure what intent to give it
+        real(kind(0d0)), dimension(-1:, 0:, 0:, 1:), intent(in) :: vSrc_rsx_vf
+        real(kind(0d0)), dimension(-1:, 0:, 0:, 1:), intent(in) :: vSrc_rsy_vf
+        real(kind(0d0)), dimension(-1:, 0:, 0:, 1:), intent(in) :: vSrc_rsz_vf
         type(scalar_field), &
             dimension(sys_size), &
-            intent(INOUT) :: flux_src_vf
-        integer :: id
+            intent(inout) :: flux_src_vf
+        integer, intent(in) :: id
+        type(int_bounds_info), intent(in) :: isx, isy, isz
 
         real(kind(0d0)), dimension(num_dims, num_dims) :: Omega
         real(kind(0d0)) :: w1L, w1R, w2L, w2R, w3L, w3R, w1, w2, w3
@@ -248,7 +248,8 @@ contains
 
     subroutine s_get_capilary(q_prim_vf)
 
-        type(scalar_field), dimension(sys_size) :: q_prim_vf
+        type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
+
         type(int_bounds_info) :: isx, isy, isz
 
         isx%beg = -1; isy%beg = 0; isz%beg = 0
@@ -321,11 +322,11 @@ contains
     subroutine s_reconstruct_cell_boundary_values_capillary(v_vf, vL_x, vL_y, vL_z, vR_x, vR_y, vR_z, &
                                                             norm_dir)
 
-        type(scalar_field), dimension(iv%beg:iv%end), intent(IN) :: v_vf
-
-        real(kind(0d0)), dimension(startx:, starty:, startz:, iv%beg:), intent(OUT) :: vL_x, vL_y, vL_z, vR_x, vR_y, vR_z
+        type(scalar_field), dimension(iv%beg:iv%end), intent(in) :: v_vf
 
-        integer, intent(IN) :: norm_dir
+        real(kind(0d0)), dimension(startx:, starty:, startz:, iv%beg:), intent(out) :: vL_x, vL_y, vL_z
+        real(kind(0d0)), dimension(startx:, starty:, startz:, iv%beg:), intent(out) :: vR_x, vR_y, vR_z
+        integer, intent(in) :: norm_dir
 
         integer :: recon_dir !< Coordinate direction of the WENO reconstruction
 
@@ -395,7 +396,7 @@ contains
 
     end subroutine s_reconstruct_cell_boundary_values_capillary
 
-    subroutine s_finalize_surface_tension_module()
+    subroutine s_finalize_surface_tension_module
 
         do j = 1, num_dims
             @:DEALLOCATE(c_divs(j)%sf)
diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp
index f600093c1b..3c10777bff 100644
--- a/src/simulation/m_time_steppers.fpp
+++ b/src/simulation/m_time_steppers.fpp
@@ -291,8 +291,8 @@ contains
         !! @param t_step Current time step
     subroutine s_1st_order_tvd_rk(t_step, time_avg)
 
-        integer, intent(IN) :: t_step
-        real(kind(0d0)), intent(INOUT) :: time_avg
+        integer, intent(in) :: t_step
+        real(kind(0d0)), intent(inout) :: time_avg
 
         integer :: i, j, k, l, q!< Generic loop iterator
         real(kind(0d0)) :: nR3bar
@@ -405,8 +405,8 @@ contains
         !! @param t_step Current time-step
     subroutine s_2nd_order_tvd_rk(t_step, time_avg)
 
-        integer, intent(IN) :: t_step
-        real(kind(0d0)), intent(INOUT) :: time_avg
+        integer, intent(in) :: t_step
+        real(kind(0d0)), intent(inout) :: time_avg
 
         integer :: i, j, k, l, q!< Generic loop iterator
         real(kind(0d0)) :: start, finish
@@ -592,9 +592,9 @@ contains
         !! @param t_step Current time-step
     subroutine s_3rd_order_tvd_rk(t_step, time_avg, dt_in)
 
-        integer, intent(IN) :: t_step
-        real(kind(0d0)), intent(INOUT) :: time_avg
-        real(kind(0d0)), intent(IN) :: dt_in
+        integer, intent(in) :: t_step
+        real(kind(0d0)), intent(inout) :: time_avg
+        real(kind(0d0)), intent(in) :: dt_in
 
         integer :: i, j, k, l, q !< Generic loop iterator
         real(kind(0d0)) :: ts_error, denom, error_fraction, time_step_factor !< Generic loop iterator
@@ -863,8 +863,8 @@ contains
         !! @param t_step Current time-step
     subroutine s_strang_splitting(t_step, time_avg)
 
-        integer, intent(IN) :: t_step
-        real(kind(0d0)), intent(INOUT) :: time_avg
+        integer, intent(in) :: t_step
+        real(kind(0d0)), intent(inout) :: time_avg
 
         integer :: i, j, k, l !< Generic loop iterator
         real(kind(0d0)) :: start, finish
@@ -893,10 +893,10 @@ contains
     end subroutine s_strang_splitting
 
     !> Bubble source part in Strang operator splitting scheme
-        !! @param q_cons_vf conservative variables
+        !! @param t_step Current time-step
     subroutine s_adaptive_dt_bubble(t_step)
 
-        integer, intent(IN) :: t_step
+        integer, intent(in) :: t_step
 
         type(int_bounds_info) :: ix, iy, iz
         type(vector_field) :: gm_alpha_qp
@@ -919,8 +919,11 @@ contains
         !! Runge-Kutta stage
     subroutine s_apply_bodyforces(q_cons_vf, q_prim_vf, rhs_vf, ldt)
 
-        type(scalar_field), dimension(1:sys_size) :: q_cons_vf, q_prim_vf, rhs_vf
-        real(kind(0d0)) :: ldt !< local dt
+        type(scalar_field), dimension(1:sys_size), intent(inout) :: q_cons_vf
+        type(scalar_field), dimension(1:sys_size), intent(in) :: q_prim_vf
+        type(scalar_field), dimension(1:sys_size), intent(inout) :: rhs_vf
+
+        real(kind(0d0)), intent(in) :: ldt !< local dt
 
         integer :: i, j, k, l
 
@@ -945,7 +948,7 @@ contains
         !! @param t_step current time-step
     subroutine s_time_step_cycling(t_step)
 
-        integer, intent(IN) :: t_step
+        integer, intent(in) :: t_step
 
         integer :: i !< Generic loop iterator
 
diff --git a/src/simulation/m_viscous.fpp b/src/simulation/m_viscous.fpp
index a6f1892b55..0019779e55 100644
--- a/src/simulation/m_viscous.fpp
+++ b/src/simulation/m_viscous.fpp
@@ -36,7 +36,7 @@ module m_viscous
 
 contains
 
-    subroutine s_initialize_viscous_module()
+    subroutine s_initialize_viscous_module
 
         integer :: i, j !< generic loop iterators
         type(int_bounds_info) :: ix, iy, iz
@@ -74,10 +74,10 @@ contains
                                                tau_Re_vf, &
                                                ix, iy, iz)
 
-        type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf
-        type(scalar_field), dimension(num_dims), intent(IN) :: grad_x_vf, grad_y_vf, grad_z_vf
-
-        type(scalar_field), dimension(1:sys_size) :: tau_Re_vf
+        type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
+        type(scalar_field), dimension(num_dims), intent(in) :: grad_x_vf, grad_y_vf, grad_z_vf
+        type(scalar_field), dimension(1:sys_size), intent(inout) :: tau_Re_vf
+        type(int_bounds_info), intent(in) :: ix, iy, iz
 
         real(kind(0d0)) :: rho_visc, gamma_visc, pi_inf_visc, alpha_visc_sum  !< Mixture variables
         real(kind(0d0)), dimension(2) :: Re_visc
@@ -87,8 +87,6 @@ contains
 
         integer :: i, j, k, l, q !< Generic loop iterator
 
-        type(int_bounds_info) :: ix, iy, iz
-
         is1_viscous = ix; is2_viscous = iy; is3_viscous = iz
 
         !$acc update device(is1_viscous, is2_viscous, is3_viscous)
@@ -539,22 +537,21 @@ contains
                              ix, iy, iz)
 
         real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), &
-            intent(INOUT) :: qL_prim_rsx_vf, qR_prim_rsx_vf, &
+            intent(inout) :: qL_prim_rsx_vf, qR_prim_rsx_vf, &
                              qL_prim_rsy_vf, qR_prim_rsy_vf, &
                              qL_prim_rsz_vf, qR_prim_rsz_vf
 
-        type(vector_field), dimension(num_dims), intent(INOUT) :: qL_prim, qR_prim
+        type(vector_field), dimension(num_dims), intent(inout) :: qL_prim, qR_prim
 
-        type(vector_field) :: q_prim_qp
+        type(vector_field), intent(in) :: q_prim_qp
 
         type(vector_field), dimension(1:num_dims), &
-            intent(INOUT) :: dqL_prim_dx_n, dqR_prim_dx_n, &
+            intent(inout) :: dqL_prim_dx_n, dqR_prim_dx_n, &
                              dqL_prim_dy_n, dqR_prim_dy_n, &
                              dqL_prim_dz_n, dqR_prim_dz_n
 
-        type(vector_field), dimension(1), intent(INOUT) :: dq_prim_dx_qp, dq_prim_dy_qp, dq_prim_dz_qp
-
-        type(int_bounds_info), intent(IN) :: ix, iy, iz
+        type(vector_field), dimension(1), intent(inout) :: dq_prim_dx_qp, dq_prim_dy_qp, dq_prim_dz_qp
+        type(int_bounds_info), intent(inout) :: ix, iy, iz
 
         integer :: i, j, k, l
 
@@ -984,18 +981,17 @@ contains
     subroutine s_reconstruct_cell_boundary_values_visc(v_vf, vL_x, vL_y, vL_z, vR_x, vR_y, vR_z, &
                                                        norm_dir, vL_prim_vf, vR_prim_vf, ix, iy, iz)
 
-        type(scalar_field), dimension(iv%beg:iv%end), intent(IN) :: v_vf
-        type(scalar_field), dimension(iv%beg:iv%end), intent(INOUT) :: vL_prim_vf, vR_prim_vf
+        type(scalar_field), dimension(iv%beg:iv%end), intent(in) :: v_vf
+        type(scalar_field), dimension(iv%beg:iv%end), intent(inout) :: vL_prim_vf, vR_prim_vf
 
-        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(INOUT) :: vL_x, vL_y, vL_z, vR_x, vR_y, vR_z
-
-        integer, intent(IN) :: norm_dir
+        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(inout) :: vL_x, vL_y, vL_z, vR_x, vR_y, vR_z
+        integer, intent(in) :: norm_dir
+        type(int_bounds_info), intent(in) :: ix, iy, iz
 
         integer :: weno_dir !< Coordinate direction of the WENO reconstruction
 
         integer :: i, j, k, l
 
-        type(int_bounds_info) :: ix, iy, iz
         ! Reconstruction in s1-direction ===================================
 
         if (norm_dir == 1) then
@@ -1087,12 +1083,10 @@ contains
     subroutine s_reconstruct_cell_boundary_values_visc_deriv(v_vf, vL_x, vL_y, vL_z, vR_x, vR_y, vR_z, &
                                                              norm_dir, vL_prim_vf, vR_prim_vf, ix, iy, iz)
 
-        type(scalar_field), dimension(iv%beg:iv%end), intent(IN) :: v_vf
-        type(scalar_field), dimension(iv%beg:iv%end), intent(INOUT) :: vL_prim_vf, vR_prim_vf
-
-        type(int_bounds_info) :: ix, iy, iz
-
-        real(kind(0d0)), dimension(startx:, starty:, startz:, iv%beg:), intent(INOUT) :: vL_x, vL_y, vL_z, vR_x, vR_y, vR_z
+        type(scalar_field), dimension(iv%beg:iv%end), intent(in) :: v_vf
+        real(kind(0d0)), dimension(startx:, starty:, startz:, iv%beg:), intent(inout) :: vL_x, vL_y, vL_z, vR_x, vR_y, vR_z
+        type(scalar_field), dimension(iv%beg:iv%end), intent(inout) :: vL_prim_vf, vR_prim_vf
+        type(int_bounds_info), intent(in) :: ix, iy, iz
 
         integer, intent(IN) :: norm_dir
 
@@ -1201,21 +1195,19 @@ contains
                                                  ix, iy, iz, iv_in, &
                                                  dL, dim, buff_size_in)
 
-        type(int_bounds_info), intent(IN) :: ix, iy, iz, iv_in
-        integer :: buff_size_in, dim
-
-        real(kind(0d0)), dimension(-buff_size_in:dim + buff_size_in) :: dL
         ! arrays of cell widths
-
         type(scalar_field), &
             dimension(iv%beg:iv%end), &
-            intent(IN) :: vL_vf, vR_vf
+            intent(in) :: vL_vf, vR_vf
 
         type(scalar_field), &
             dimension(iv%beg:iv%end), &
-            intent(INOUT) :: dv_ds_vf
+            intent(inout) :: dv_ds_vf
 
-        integer, intent(IN) :: norm_dir
+        integer, intent(in) :: norm_dir
+        type(int_bounds_info), intent(in) :: ix, iy, iz, iv_in
+        integer, intent(in) :: dim, buff_size_in
+        real(kind(0d0)), dimension(-buff_size_in:dim + buff_size_in), intent(in) :: dL
 
         integer :: i, j, k, l !< Generic loop iterators
 
@@ -1322,17 +1314,15 @@ contains
     subroutine s_compute_fd_gradient(var, grad_x, grad_y, grad_z, &
                                      ix, iy, iz, buff_size_in)
 
-        type(scalar_field), intent(IN) :: var
-        type(scalar_field), intent(INOUT) :: grad_x
-        type(scalar_field), intent(INOUT) :: grad_y
-        type(scalar_field), intent(INOUT) :: grad_z
-
-        integer, intent(IN) :: buff_size_in
+        type(scalar_field), intent(in) :: var
+        type(scalar_field), intent(inout) :: grad_x
+        type(scalar_field), intent(inout) :: grad_y
+        type(scalar_field), intent(inout) :: grad_z
+        type(int_bounds_info), intent(inout) :: ix, iy, iz
+        integer, intent(in) :: buff_size_in
 
         integer :: j, k, l !< Generic loop iterators
 
-        type(int_bounds_info) :: ix, iy, iz
-
         ix%beg = -buff_size_in; ix%end = m + buff_size_in; 
         if (n > 0) then
             iy%beg = -buff_size_in; iy%end = n + buff_size_in
diff --git a/src/simulation/m_weno.fpp b/src/simulation/m_weno.fpp
index 8201685b86..6a05b9d9c3 100644
--- a/src/simulation/m_weno.fpp
+++ b/src/simulation/m_weno.fpp
@@ -267,8 +267,8 @@ contains
         !! @param is Index bounds in the s-direction
     subroutine s_compute_weno_coefficients(weno_dir, is)
 
-        integer, intent(IN) :: weno_dir
-        type(int_bounds_info), intent(IN) :: is
+        integer, intent(in) :: weno_dir
+        type(int_bounds_info), intent(in) :: is
         integer :: s
 
         real(kind(0d0)), pointer, dimension(:) :: s_cb => null() !<
@@ -506,11 +506,12 @@ contains
                       norm_dir, weno_dir, &
                       is1_weno_d, is2_weno_d, is3_weno_d)
 
-        type(scalar_field), dimension(1:), intent(IN) :: v_vf
-        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(INOUT) :: vL_rs_vf_x, vL_rs_vf_y, vL_rs_vf_z, vR_rs_vf_x, vR_rs_vf_y, vR_rs_vf_z
-        integer, intent(IN) :: norm_dir
-        integer, intent(IN) :: weno_dir
-        type(int_bounds_info), intent(IN) :: is1_weno_d, is2_weno_d, is3_weno_d
+        type(scalar_field), dimension(1:), intent(in) :: v_vf
+        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(inout) :: vL_rs_vf_x, vL_rs_vf_y, vL_rs_vf_z
+        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(inout) :: vR_rs_vf_x, vR_rs_vf_y, vR_rs_vf_z
+        integer, intent(in) :: norm_dir
+        integer, intent(in) :: weno_dir
+        type(int_bounds_info), intent(in) :: is1_weno_d, is2_weno_d, is3_weno_d
 
         real(kind(0d0)), dimension(-weno_polyn:weno_polyn - 1) :: dvd
         real(kind(0d0)), dimension(0:weno_polyn) :: poly