diff --git a/src/common/m_helper.fpp b/src/common/m_helper.fpp
index 40a1f562d..b6d487297 100644
--- a/src/common/m_helper.fpp
+++ b/src/common/m_helper.fpp
@@ -104,7 +104,7 @@ contains
 
         end if
 
-    end subroutine s_compute_finite_difference_coefficients ! --------------
+    end subroutine s_compute_finite_difference_coefficients
 
     !> Computes the bubble number density n from the primitive variables
         !! @param vftmp is the void fraction
diff --git a/src/common/m_mpi_common.fpp b/src/common/m_mpi_common.fpp
index 951cc3ca5..ab86ce52c 100644
--- a/src/common/m_mpi_common.fpp
+++ b/src/common/m_mpi_common.fpp
@@ -30,7 +30,7 @@ contains
     !> The subroutine initializes the MPI execution environment
         !!      and queries both the number of processors which will be
         !!      available for the job and the local processor rank.
-    subroutine s_mpi_initialize ! ----------------------------------------
+    subroutine s_mpi_initialize
 
 #ifndef MFC_MPI
 
@@ -58,9 +58,11 @@ contains
 
 #endif
 
-    end subroutine s_mpi_initialize ! --------------------------------------
+    end subroutine s_mpi_initialize
 
-    subroutine s_initialize_mpi_data(q_cons_vf, ib_markers) ! --------------------------
+    !! @param q_cons_vf Conservative variables
+    !! @param ib_markers track if a cell is within the immersed boundary
+    subroutine s_initialize_mpi_data(q_cons_vf, ib_markers)
 
         type(scalar_field), &
             dimension(sys_size), &
@@ -188,7 +190,7 @@ contains
 
 #endif
 
-    end subroutine s_initialize_mpi_data ! ---------------------------------
+    end subroutine s_initialize_mpi_data
 
     subroutine mpi_bcast_time_step_values(proc_time, time_avg)
 
@@ -217,7 +219,7 @@ contains
         !!  @param icfl_max_glb Global maximum ICFL stability criterion
         !!  @param vcfl_max_glb Global maximum VCFL stability criterion
         !!  @param Rc_min_glb Global minimum Rc stability criterion
-    subroutine s_mpi_reduce_stability_criteria_extrema(icfl_max_loc, & ! --
+    subroutine s_mpi_reduce_stability_criteria_extrema(icfl_max_loc, &
                                                        vcfl_max_loc, &
                                                        ccfl_max_loc, &
                                                        Rc_min_loc, &
@@ -257,7 +259,7 @@ contains
 #endif
 #endif
 
-    end subroutine s_mpi_reduce_stability_criteria_extrema ! ---------------
+    end subroutine s_mpi_reduce_stability_criteria_extrema
 
     !>  The following subroutine takes the input local variable
         !!      from all processors and reduces to the sum of all
@@ -266,7 +268,7 @@ contains
         !!  @param var_loc Some variable containing the local value which should be
         !!  reduced amongst all the processors in the communicator.
         !!  @param var_glb The globally reduced value
-    subroutine s_mpi_allreduce_sum(var_loc, var_glb) ! ---------------------
+    subroutine s_mpi_allreduce_sum(var_loc, var_glb)
 
         real(kind(0d0)), intent(in) :: var_loc
         real(kind(0d0)), intent(out) :: var_glb
@@ -279,7 +281,7 @@ contains
 
 #endif
 
-    end subroutine s_mpi_allreduce_sum ! -----------------------------------
+    end subroutine s_mpi_allreduce_sum
 
     !>  The following subroutine takes the input local variable
         !!      from all processors and reduces to the minimum of all
@@ -288,7 +290,7 @@ contains
         !!  @param var_loc Some variable containing the local value which should be
         !!  reduced amongst all the processors in the communicator.
         !!  @param var_glb The globally reduced value
-    subroutine s_mpi_allreduce_min(var_loc, var_glb) ! ---------------------
+    subroutine s_mpi_allreduce_min(var_loc, var_glb)
 
         real(kind(0d0)), intent(in) :: var_loc
         real(kind(0d0)), intent(out) :: var_glb
@@ -301,7 +303,7 @@ contains
 
 #endif
 
-    end subroutine s_mpi_allreduce_min ! -----------------------------------
+    end subroutine s_mpi_allreduce_min
 
     !>  The following subroutine takes the input local variable
         !!      from all processors and reduces to the maximum of all
@@ -310,7 +312,7 @@ contains
         !!  @param var_loc Some variable containing the local value which should be
         !!  reduced amongst all the processors in the communicator.
         !!  @param var_glb The globally reduced value
-    subroutine s_mpi_allreduce_max(var_loc, var_glb) ! ---------------------
+    subroutine s_mpi_allreduce_max(var_loc, var_glb)
 
         real(kind(0d0)), intent(in) :: var_loc
         real(kind(0d0)), intent(out) :: var_glb
@@ -323,7 +325,7 @@ contains
 
 #endif
 
-    end subroutine s_mpi_allreduce_max ! -----------------------------------
+    end subroutine s_mpi_allreduce_max
 
     !>  The following subroutine takes the inputted variable and
         !!      determines its minimum value on the entire computational
@@ -331,7 +333,7 @@ contains
         !!  @param var_loc holds the local value to be reduced among
         !!      all the processors in communicator. On output, the variable holds
         !!      the minimum value, reduced amongst all of the local values.
-    subroutine s_mpi_reduce_min(var_loc) ! ---------------------------------
+    subroutine s_mpi_reduce_min(var_loc)
 
         real(kind(0d0)), intent(inout) :: var_loc
 
@@ -352,7 +354,7 @@ contains
 
 #endif
 
-    end subroutine s_mpi_reduce_min ! --------------------------------------
+    end subroutine s_mpi_reduce_min
 
     !>  The following subroutine takes the first element of the
         !!      2-element inputted variable and determines its maximum
@@ -366,7 +368,7 @@ contains
         !!  On output, this variable holds the maximum value, reduced amongst
         !!  all of the local values, and the process rank to which the value
         !!  belongs.
-    subroutine s_mpi_reduce_maxloc(var_loc) ! ------------------------------
+    subroutine s_mpi_reduce_maxloc(var_loc)
 
         real(kind(0d0)), dimension(2), intent(inout) :: var_loc
 
@@ -388,10 +390,11 @@ contains
 
 #endif
 
-    end subroutine s_mpi_reduce_maxloc ! -----------------------------------
+    end subroutine s_mpi_reduce_maxloc
 
     !> The subroutine terminates the MPI execution environment.
-    subroutine s_mpi_abort(prnt) ! ---------------------------------------------
+        !! @param prnt error message to be printed
+    subroutine s_mpi_abort(prnt)
 
         character(len=*), intent(in), optional :: prnt
 
@@ -412,10 +415,10 @@ contains
 
 #endif
 
-    end subroutine s_mpi_abort ! -------------------------------------------
+    end subroutine s_mpi_abort
 
     !>Halts all processes until all have reached barrier.
-    subroutine s_mpi_barrier ! -------------------------------------------
+    subroutine s_mpi_barrier
 
 #ifdef MFC_MPI
 
@@ -424,10 +427,10 @@ contains
 
 #endif
 
-    end subroutine s_mpi_barrier ! -----------------------------------------
+    end subroutine s_mpi_barrier
 
     !> The subroutine finalizes the MPI execution environment.
-    subroutine s_mpi_finalize ! ------------------------------------------
+    subroutine s_mpi_finalize
 
 #ifdef MFC_MPI
 
@@ -436,6 +439,6 @@ contains
 
 #endif
 
-    end subroutine s_mpi_finalize ! ----------------------------------------
+    end subroutine s_mpi_finalize
 
 end module m_mpi_common
diff --git a/src/common/m_phase_change.fpp b/src/common/m_phase_change.fpp
index 33737d107..f63723228 100644
--- a/src/common/m_phase_change.fpp
+++ b/src/common/m_phase_change.fpp
@@ -28,10 +28,11 @@ module m_phase_change
 
     implicit none
 
-    private; public :: s_initialize_phasechange_module, &
- s_relaxation_solver, &
- s_infinite_relaxation_k, &
- s_finalize_relaxation_solver_module
+    private; 
+    public :: s_initialize_phasechange_module, &
+              s_relaxation_solver, &
+              s_infinite_relaxation_k, &
+              s_finalize_relaxation_solver_module
 
     !> @name Abstract interface for creating function pointers
     !> @{
@@ -39,7 +40,7 @@ module m_phase_change
 
         !> @name Abstract subroutine for the infinite relaxation solver
         !> @{
-        subroutine s_abstract_relaxation_solver(q_cons_vf) ! -------
+        subroutine s_abstract_relaxation_solver(q_cons_vf)
             import :: scalar_field, sys_size
             type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
         end subroutine
@@ -86,14 +87,14 @@ contains
         D = ((gs_min(lp) - 1.0d0)*cvs(lp)) &
             /((gs_min(vp) - 1.0d0)*cvs(vp))
 
-    end subroutine s_initialize_phasechange_module !-------------------------------
+    end subroutine s_initialize_phasechange_module
 
     !>  This subroutine is created to activate either the pT- (N fluids) or the
         !!      pTg-equilibrium (2 fluids for g-equilibrium)
         !!      model, also considering mass depletion, depending on the incoming
         !!      state conditions.
         !!  @param q_cons_vf Cell-average conservative variables
-    subroutine s_infinite_relaxation_k(q_cons_vf) ! ----------------
+    subroutine s_infinite_relaxation_k(q_cons_vf)
 
         type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
         real(kind(0.0d0)) :: pS, pSOV, pSSL !< equilibrium pressure for mixture, overheated vapor, and subcooled liquid
@@ -283,7 +284,7 @@ contains
             end do
         end do
 
-    end subroutine s_infinite_relaxation_k ! ----------------
+    end subroutine s_infinite_relaxation_k
 
     !>  This auxiliary subroutine is created to activate the pT-equilibrium for N fluids
         !!  @param j generic loop iterator for x direction
@@ -387,7 +388,7 @@ contains
         ! common temperature
         TS = (rhoe + pS - mQ)/mCP
 
-    end subroutine s_infinite_pt_relaxation_k ! -----------------------
+    end subroutine s_infinite_pt_relaxation_k
 
     !>  This auxiliary subroutine is created to activate the pTg-equilibrium for N fluids under pT
         !!      and 2 fluids under pTg-equilibrium. There is a final common p and T during relaxation
@@ -510,7 +511,7 @@ contains
 
         ! common temperature
         TS = (rhoe + pS - mQ)/mCP
-    end subroutine s_infinite_ptg_relaxation_k ! -----------------------
+    end subroutine s_infinite_ptg_relaxation_k
 
     !>  This auxiliary subroutine corrects the partial densities of the REACTING fluids in case one of them is negative
         !!      but their sum is positive. Inert phases are not corrected at this moment
@@ -782,7 +783,7 @@ contains
     end subroutine s_TSat
 
     !>  This subroutine finalizes the phase change module
-    subroutine s_finalize_relaxation_solver_module()
+    subroutine s_finalize_relaxation_solver_module
     end subroutine
 
 #endif
diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp
index 328cf500a..6228ad50d 100644
--- a/src/common/m_variables_conversion.fpp
+++ b/src/common/m_variables_conversion.fpp
@@ -216,7 +216,7 @@ contains
         qv_sf(i, j, k) = qv
 #endif
 
-    end subroutine s_convert_mixture_to_mixture_variables ! ----------------
+    end subroutine s_convert_mixture_to_mixture_variables
 
     !>  This procedure is used alongside with the gamma/pi_inf
         !!      model to transfer the density, the specific heat ratio
@@ -341,7 +341,7 @@ contains
         qv_sf(j, k, l) = qv
 #endif
 
-    end subroutine s_convert_species_to_mixture_variables_bubbles ! ----------------
+    end subroutine s_convert_species_to_mixture_variables_bubbles
 
     !>  This subroutine is designed for the volume fraction model
         !!              and provided a set of either conservative or primitive
@@ -440,7 +440,7 @@ contains
         qv_sf(k, l, r) = qv
 #endif
 
-    end subroutine s_convert_species_to_mixture_variables ! ----------------
+    end subroutine s_convert_species_to_mixture_variables
 
     subroutine s_convert_species_to_mixture_variables_acc(rho_K, &
                                                           gamma_K, pi_inf_K, qv_K, &
@@ -523,7 +523,7 @@ contains
         end if
 #endif
 
-    end subroutine s_convert_species_to_mixture_variables_acc ! ----------------
+    end subroutine s_convert_species_to_mixture_variables_acc
 
     subroutine s_convert_species_to_mixture_variables_bubbles_acc(rho_K, &
                                                                   gamma_K, pi_inf_K, qv_K, &
@@ -596,7 +596,7 @@ contains
     !>  The computation of parameters, the allocation of memory,
         !!      the association of pointers and/or the execution of any
         !!      other procedures that are necessary to setup the module.
-    subroutine s_initialize_variables_conversion_module ! ----------------
+    subroutine s_initialize_variables_conversion_module
 
         integer :: i, j
 
@@ -752,7 +752,7 @@ contains
             s_convert_to_mixture_variables => &
                 s_convert_species_to_mixture_variables
         end if
-    end subroutine s_initialize_variables_conversion_module ! --------------
+    end subroutine s_initialize_variables_conversion_module
 
     !Initialize mv at the quadrature nodes based on the initialized moments and sigma
     subroutine s_initialize_mv(qK_cons_vf, mv)
@@ -1014,7 +1014,7 @@ contains
         end do
         !$acc end parallel loop
 
-    end subroutine s_convert_conservative_to_primitive_variables ! ---------
+    end subroutine s_convert_conservative_to_primitive_variables
 
     !>  The following procedure handles the conversion between
         !!      the primitive variables and the conservative variables.
@@ -1175,7 +1175,7 @@ contains
         end if
 #endif
 
-    end subroutine s_convert_primitive_to_conservative_variables ! ---------
+    end subroutine s_convert_primitive_to_conservative_variables
 
     !>  The following subroutine handles the conversion between
         !!      the primitive variables and the Eulerian flux variables.
@@ -1185,7 +1185,7 @@ contains
         !!  @param ix Index bounds in the first coordinate direction
         !!  @param iy Index bounds in the second coordinate direction
         !!  @param iz Index bounds in the third coordinate direction
-    subroutine s_convert_primitive_to_flux_variables(qK_prim_vf, & ! ------
+    subroutine s_convert_primitive_to_flux_variables(qK_prim_vf, &
                                                      FK_vf, &
                                                      FK_src_vf, &
                                                      is1, is2, is3, s2b, s3b)
@@ -1309,9 +1309,9 @@ contains
         end do
 #endif
 
-    end subroutine s_convert_primitive_to_flux_variables ! -----------------
+    end subroutine s_convert_primitive_to_flux_variables
 
-    subroutine s_finalize_variables_conversion_module ! ------------------
+    subroutine s_finalize_variables_conversion_module
 
         ! Deallocating the density, the specific heat ratio function and the
         ! liquid stiffness function
@@ -1335,6 +1335,6 @@ contains
         ! computing the mixture/species variables to the mixture variables
         s_convert_to_mixture_variables => null()
 
-    end subroutine s_finalize_variables_conversion_module ! ----------------
+    end subroutine s_finalize_variables_conversion_module
 
 end module m_variables_conversion
diff --git a/src/post_process/m_data_input.f90 b/src/post_process/m_data_input.f90
index 0a0982607..9b06e223c 100644
--- a/src/post_process/m_data_input.f90
+++ b/src/post_process/m_data_input.f90
@@ -37,11 +37,11 @@ module m_data_input
 
         !> Subroutine for reading data files
         !!  @param t_step Current time-step to input
-        subroutine s_read_abstract_data_files(t_step) ! ------------
+        subroutine s_read_abstract_data_files(t_step)
 
             integer, intent(IN) :: t_step
 
-        end subroutine s_read_abstract_data_files ! ----------------
+        end subroutine s_read_abstract_data_files
 
     end interface ! ========================================================
 
@@ -60,7 +60,7 @@ end subroutine s_read_abstract_data_files ! ----------------
         !!      present in the corresponding time-step directory and to
         !!      populate the associated grid and conservative variables.
         !!  @param t_step Current time-step
-    subroutine s_read_serial_data_files(t_step) ! -----------------------------
+    subroutine s_read_serial_data_files(t_step)
 
         integer, intent(IN) :: t_step
 
@@ -207,14 +207,14 @@ subroutine s_read_serial_data_files(t_step) ! -----------------------------
 
         ! ==================================================================
 
-    end subroutine s_read_serial_data_files ! ---------------------------------
+    end subroutine s_read_serial_data_files
 
     !>  This subroutine is called at each time-step that has to
         !!      be post-processed in order to parallel-read the raw data files
         !!      present in the corresponding time-step directory and to
         !!      populate the associated grid and conservative variables.
         !!  @param t_step Current time-step
-    subroutine s_read_parallel_data_files(t_step) ! ---------------------------
+    subroutine s_read_parallel_data_files(t_step)
 
         integer, intent(IN) :: t_step
 
@@ -417,7 +417,7 @@ subroutine s_read_parallel_data_files(t_step) ! ---------------------------
 
 #endif
 
-    end subroutine s_read_parallel_data_files ! -------------------------------
+    end subroutine s_read_parallel_data_files
 
     !>  The following subroutine populates the buffer regions of
         !!      the cell-width spacings, the cell-boundary locations and
@@ -427,7 +427,7 @@ end subroutine s_read_parallel_data_files ! -------------------------------
         !!      are used in aiding the multidimensional visualization of
         !!      Silo database files, in VisIt, when processor boundary
         !!      conditions are present.
-    subroutine s_populate_grid_variables_buffer_regions() ! ----------------
+    subroutine s_populate_grid_variables_buffer_regions
 
         integer :: i !< Generic loop iterator
 
@@ -667,12 +667,12 @@ subroutine s_populate_grid_variables_buffer_regions() ! ----------------
 
         ! END: Populating Buffer Regions in the z-direction ================
 
-    end subroutine s_populate_grid_variables_buffer_regions ! --------------
+    end subroutine s_populate_grid_variables_buffer_regions
 
     !>  The purpose of this procedure is to populate the buffers
         !!      of the cell-average conservative variables, depending on
         !!      the boundary conditions.
-    subroutine s_populate_conservative_variables_buffer_regions() ! --------
+    subroutine s_populate_conservative_variables_buffer_regions
 
         integer :: i, j, k !< Generic loop iterators
 
@@ -1046,11 +1046,11 @@ subroutine s_populate_conservative_variables_buffer_regions() ! --------
 
         ! END: Populating Buffer Regions in the z-direction ================
 
-    end subroutine s_populate_conservative_variables_buffer_regions ! ------
+    end subroutine s_populate_conservative_variables_buffer_regions
 
     !>  Computation of parameters, allocation procedures, and/or
         !!      any other tasks needed to properly setup the module
-    subroutine s_initialize_data_input_module() ! -----------------------------
+    subroutine s_initialize_data_input_module
 
         integer :: i !< Generic loop iterator
 
@@ -1113,10 +1113,10 @@ subroutine s_initialize_data_input_module() ! -----------------------------
             s_read_data_files => s_read_parallel_data_files
         end if
 
-    end subroutine s_initialize_data_input_module ! ---------------------------
+    end subroutine s_initialize_data_input_module
 
     !> Deallocation procedures for the module
-    subroutine s_finalize_data_input_module() ! --------------------------
+    subroutine s_finalize_data_input_module()
 
         integer :: i !< Generic loop iterator
 
@@ -1131,6 +1131,6 @@ subroutine s_finalize_data_input_module() ! --------------------------
 
         s_read_data_files => null()
 
-    end subroutine s_finalize_data_input_module ! ------------------------
+    end subroutine s_finalize_data_input_module
 
 end module m_data_input
diff --git a/src/post_process/m_data_output.fpp b/src/post_process/m_data_output.fpp
index 7f3d7d9bd..ca01e0d95 100644
--- a/src/post_process/m_data_output.fpp
+++ b/src/post_process/m_data_output.fpp
@@ -102,7 +102,7 @@ module m_data_output
 
 contains
 
-    subroutine s_initialize_data_output_module() ! ----------------------------
+    subroutine s_initialize_data_output_module()
         ! Description: Computation of parameters, allocation procedures, and/or
         !              any other tasks needed to properly setup the module
 
@@ -410,9 +410,9 @@ contains
 
         ! END: Querying Number of Flow Variable(s) in Binary Output ========
 
-    end subroutine s_initialize_data_output_module ! --------------------------
+    end subroutine s_initialize_data_output_module
 
-    subroutine s_open_formatted_database_file(t_step) ! --------------------
+    subroutine s_open_formatted_database_file(t_step)
         ! Description: This subroutine opens a new formatted database file, or
         !              replaces an old one, and readies it for the data storage
         !              of the grid and the flow variable(s) associated with the
@@ -529,9 +529,9 @@ contains
 
         ! END: Binary Database Format ======================================
 
-    end subroutine s_open_formatted_database_file ! ------------------------
+    end subroutine s_open_formatted_database_file
 
-    subroutine s_write_grid_to_formatted_database_file(t_step) ! -----------
+    subroutine s_write_grid_to_formatted_database_file(t_step)
         ! Description: The general objective of this subroutine is to write the
         !              necessary grid data to the formatted database file, for
         !              the current time-step, t_step. The local processor will
@@ -728,7 +728,7 @@ contains
 
         ! ==================================================================
 
-    end subroutine s_write_grid_to_formatted_database_file ! ---------------
+    end subroutine s_write_grid_to_formatted_database_file
 
     subroutine s_write_variable_to_formatted_database_file(varname, t_step)
         ! Description: The goal of this subroutine is to write to the formatted
@@ -937,9 +937,9 @@ contains
 
         ! ==================================================================
 
-    end subroutine s_write_variable_to_formatted_database_file ! -----------
+    end subroutine s_write_variable_to_formatted_database_file
 
-    subroutine s_close_formatted_database_file() ! -------------------------
+    subroutine s_close_formatted_database_file
         ! Description: The purpose of this subroutine is to close any formatted
         !              database file(s) that may be opened at the time-step that
         !              is currently being post-processed. The root process must
@@ -964,9 +964,9 @@ contains
 
         end if
 
-    end subroutine s_close_formatted_database_file ! -----------------------
+    end subroutine s_close_formatted_database_file
 
-    subroutine s_finalize_data_output_module() ! -------------------------
+    subroutine s_finalize_data_output_module
         ! Description: Deallocation procedures for the module
 
         ! Deallocating the generic storage employed for the flow variable(s)
@@ -990,6 +990,6 @@ contains
             deallocate (dims)
         end if
 
-    end subroutine s_finalize_data_output_module ! -----------------------
+    end subroutine s_finalize_data_output_module
 
 end module m_data_output
diff --git a/src/post_process/m_derived_variables.fpp b/src/post_process/m_derived_variables.fpp
index 1c690ad37..d64dfaf25 100644
--- a/src/post_process/m_derived_variables.fpp
+++ b/src/post_process/m_derived_variables.fpp
@@ -63,7 +63,7 @@ contains
 
     !>  Computation of parameters, allocation procedures, and/or
         !!      any other tasks needed to properly setup the module
-    subroutine s_initialize_derived_variables_module() ! ----------------------
+    subroutine s_initialize_derived_variables_module
 
         ! Allocating the gradient magnitude of the density variable provided
         ! that numerical Schlieren function is outputted during post-process
@@ -112,14 +112,14 @@ contains
             flg = 0
         end if
 
-    end subroutine s_initialize_derived_variables_module ! --------------------
+    end subroutine s_initialize_derived_variables_module
 
     !>  This subroutine receives as input the specific heat ratio
         !!      function, gamma_sf, and derives from it the specific heat
         !!      ratio. The latter is stored in the derived flow quantity
         !!      storage variable, q_sf.
         !!  @param q_sf Specific heat ratio
-    subroutine s_derive_specific_heat_ratio(q_sf) ! --------------
+    subroutine s_derive_specific_heat_ratio(q_sf)
 
         real(kind(0d0)), &
             dimension(-offset_x%beg:m + offset_x%end, &
@@ -138,7 +138,7 @@ contains
             end do
         end do
 
-    end subroutine s_derive_specific_heat_ratio ! --------------------------
+    end subroutine s_derive_specific_heat_ratio
 
     !>  This subroutine admits as inputs the specific heat ratio
         !!      function and the liquid stiffness function, gamma_sf and
@@ -146,7 +146,7 @@ contains
         !!      values of the liquid stiffness, which are stored in the
         !!      derived flow quantity storage variable, q_sf.
         !!  @param q_sf Liquid stiffness
-    subroutine s_derive_liquid_stiffness(q_sf) ! ------
+    subroutine s_derive_liquid_stiffness(q_sf)
 
         real(kind(0d0)), &
             dimension(-offset_x%beg:m + offset_x%end, &
@@ -166,7 +166,7 @@ contains
             end do
         end do
 
-    end subroutine s_derive_liquid_stiffness ! -----------------------------
+    end subroutine s_derive_liquid_stiffness
 
     !> This subroutine admits as inputs the primitive variables,
         !!      the density, the specific heat ratio function and liquid
@@ -222,7 +222,7 @@ contains
             end do
         end do
 
-    end subroutine s_derive_sound_speed ! ----------------------------------
+    end subroutine s_derive_sound_speed
 
     !>  This subroutine derives the flux_limiter at cell boundary
         !!      i+1/2. This is an approximation because the velocity used
@@ -232,7 +232,7 @@ contains
         !!  @param i Component indicator
         !!  @param q_prim_vf Primitive variables
         !!  @param q_sf Flux limiter
-    subroutine s_derive_flux_limiter(i, q_prim_vf, q_sf) ! -----------------
+    subroutine s_derive_flux_limiter(i, q_prim_vf, q_sf)
 
         integer, intent(IN) :: i
 
@@ -319,7 +319,7 @@ contains
                 end do
             end do
         end do
-    end subroutine s_derive_flux_limiter ! ---------------------------------
+    end subroutine s_derive_flux_limiter
 
     !>  Computes the solution to the linear system Ax=b w/ sol = x
         !!  @param A Input matrix
@@ -367,7 +367,7 @@ contains
             end do
         end do
 
-    end subroutine s_solve_linear_system ! -------------------------------------
+    end subroutine s_solve_linear_system
 
     !>  This subroutine receives as inputs the indicator of the
         !!      component of the vorticity that should be outputted and
@@ -378,7 +378,7 @@ contains
         !!  @param i Vorticity component indicator
         !!  @param q_prim_vf Primitive variables
         !!  @param q_sf Vorticity component
-    subroutine s_derive_vorticity_component(i, q_prim_vf, q_sf) ! ----------
+    subroutine s_derive_vorticity_component(i, q_prim_vf, q_sf)
 
         integer, intent(IN) :: i
 
@@ -472,7 +472,7 @@ contains
             end do
         end if
 
-    end subroutine s_derive_vorticity_component ! --------------------------
+    end subroutine s_derive_vorticity_component
 
     !> This subroutine gets as inputs the primitive variables. From those
         !!      inputs, it proceeds to calculate the value of the Q_M
@@ -714,10 +714,10 @@ contains
 
         ! ==================================================================
 
-    end subroutine s_derive_numerical_schlieren_function ! -----------------
+    end subroutine s_derive_numerical_schlieren_function
 
     !>  Deallocation procedures for the module
-    subroutine s_finalize_derived_variables_module() ! -------------------
+    subroutine s_finalize_derived_variables_module
 
         ! Deallocating the variable containing the gradient magnitude of the
         ! density field provided that the numerical Schlieren function was
@@ -730,6 +730,6 @@ contains
         if (allocated(fd_coeff_y)) deallocate (fd_coeff_y)
         if (allocated(fd_coeff_z)) deallocate (fd_coeff_z)
 
-    end subroutine s_finalize_derived_variables_module ! -----------------
+    end subroutine s_finalize_derived_variables_module
 
 end module m_derived_variables
diff --git a/src/post_process/m_global_parameters.fpp b/src/post_process/m_global_parameters.fpp
index 5c88ca9f0..5c5ed3ef1 100644
--- a/src/post_process/m_global_parameters.fpp
+++ b/src/post_process/m_global_parameters.fpp
@@ -259,7 +259,7 @@ contains
     !> Assigns default values to user inputs prior to reading
         !!      them in. This allows for an easier consistency check of
         !!      these parameters once they are read from the input file.
-    subroutine s_assign_default_values_to_user_inputs() ! ------------------
+    subroutine s_assign_default_values_to_user_inputs
 
         integer :: i !< Generic loop iterator
 
@@ -354,11 +354,11 @@ contains
         sigma = dflt_real
         adv_n = .false.
 
-    end subroutine s_assign_default_values_to_user_inputs ! ----------------
+    end subroutine s_assign_default_values_to_user_inputs
 
     !>  Computation of parameters, allocation procedures, and/or
         !!      any other tasks needed to properly setup the module
-    subroutine s_initialize_global_parameters_module() ! ----------------------
+    subroutine s_initialize_global_parameters_module
 
         integer :: i, j, fac
 
@@ -683,10 +683,10 @@ contains
             grid_geometry = 3
         end if
 
-    end subroutine s_initialize_global_parameters_module ! --------------------
+    end subroutine s_initialize_global_parameters_module
 
     !> Subroutine to initialize parallel infrastructure
-    subroutine s_initialize_parallel_io() ! --------------------------------
+    subroutine s_initialize_parallel_io
 
         num_dims = 1 + min(1, n) + min(1, p)
 
@@ -711,10 +711,10 @@ contains
 
 #endif
 
-    end subroutine s_initialize_parallel_io ! ------------------------------
+    end subroutine s_initialize_parallel_io
 
     !> Deallocation procedures for the module
-    subroutine s_finalize_global_parameters_module() ! -------------------
+    subroutine s_finalize_global_parameters_module
 
         integer :: i
 
@@ -754,6 +754,6 @@ contains
 
 #endif
 
-    end subroutine s_finalize_global_parameters_module ! -----------------
+    end subroutine s_finalize_global_parameters_module
 
 end module m_global_parameters
diff --git a/src/post_process/m_mpi_proxy.fpp b/src/post_process/m_mpi_proxy.fpp
index cd175ef7a..b7cee50e7 100644
--- a/src/post_process/m_mpi_proxy.fpp
+++ b/src/post_process/m_mpi_proxy.fpp
@@ -50,7 +50,7 @@ contains
 
     !>  Computation of parameters, allocation procedures, and/or
         !!      any other tasks needed to properly setup the module
-    subroutine s_initialize_mpi_proxy_module() ! ------------------------------
+    subroutine s_initialize_mpi_proxy_module
 
 #ifdef MFC_MPI
 
@@ -139,14 +139,14 @@ contains
 
 #endif
 
-    end subroutine s_initialize_mpi_proxy_module ! ----------------------------
+    end subroutine s_initialize_mpi_proxy_module
 
     !>  Since only processor with rank 0 is in charge of reading
         !!      and checking the consistency of the user provided inputs,
         !!      these are not available to the remaining processors. This
         !!      subroutine is then in charge of broadcasting the required
         !!      information.
-    subroutine s_mpi_bcast_user_inputs() ! ---------------------------------
+    subroutine s_mpi_bcast_user_inputs
 
 #ifdef MFC_MPI
         integer :: i !< Generic loop iterator
@@ -195,14 +195,14 @@ contains
         call MPI_BCAST(schlieren_alpha(1), num_fluids_max, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
 #endif
 
-    end subroutine s_mpi_bcast_user_inputs ! -------------------------------
+    end subroutine s_mpi_bcast_user_inputs
 
     !>  This subroutine takes care of efficiently distributing
         !!      the computational domain among the available processors
         !!      as well as recomputing some of the global parameters so
         !!      that they reflect the configuration of sub-domain that
         !!      is overseen by the local processor.
-    subroutine s_mpi_decompose_computational_domain() ! --------------------
+    subroutine s_mpi_decompose_computational_domain
 
 #ifdef MFC_MPI
 
@@ -638,7 +638,7 @@ contains
 
 #endif
 
-    end subroutine s_mpi_decompose_computational_domain ! ------------------
+    end subroutine s_mpi_decompose_computational_domain
 
     !>  Communicates the buffer regions associated with the grid
         !!      variables with processors in charge of the neighboring
@@ -838,7 +838,7 @@ contains
 
 #endif
 
-    end subroutine s_mpi_sendrecv_grid_vars_buffer_regions ! ---------------
+    end subroutine s_mpi_sendrecv_grid_vars_buffer_regions
 
     !>  Communicates buffer regions associated with conservative
         !!      variables with processors in charge of the neighboring
@@ -1403,7 +1403,7 @@ contains
 
 #endif
 
-    end subroutine s_mpi_sendrecv_cons_vars_buffer_regions ! ---------------
+    end subroutine s_mpi_sendrecv_cons_vars_buffer_regions
 
     !>  This subroutine gathers the Silo database metadata for
         !!      the spatial extents in order to boost the performance of
@@ -1411,7 +1411,7 @@ contains
         !!  @param spatial_extents Spatial extents for each processor's sub-domain. First dimension
         !!  corresponds to the minimum and maximum values, respectively, while
         !!  the second dimension corresponds to the processor rank.
-    subroutine s_mpi_gather_spatial_extents(spatial_extents) ! -------------
+    subroutine s_mpi_gather_spatial_extents(spatial_extents)
 
         real(kind(0d0)), dimension(1:, 0:), intent(INOUT) :: spatial_extents
 
@@ -1523,14 +1523,14 @@ contains
 
 #endif
 
-    end subroutine s_mpi_gather_spatial_extents ! --------------------------
+    end subroutine s_mpi_gather_spatial_extents
 
     !>  This subroutine collects the sub-domain cell-boundary or
         !!      cell-center locations data from all of the processors and
         !!      puts back together the grid of the entire computational
         !!      domain on the rank 0 processor. This is only done for 1D
         !!      simulations.
-    subroutine s_mpi_defragment_1d_grid_variable() ! -----------------------
+    subroutine s_mpi_defragment_1d_grid_variable
 
 #ifdef MFC_MPI
 
@@ -1556,7 +1556,7 @@ contains
 
 #endif
 
-    end subroutine s_mpi_defragment_1d_grid_variable ! ---------------------
+    end subroutine s_mpi_defragment_1d_grid_variable
 
     !>  This subroutine gathers the Silo database metadata for
         !!      the flow variable's extents as to boost performance of
@@ -1566,7 +1566,7 @@ contains
         !!   First dimension of array corresponds to the former's minimum and
         !!  maximum values, respectively, while second dimension corresponds
         !!  to each processor's rank.
-    subroutine s_mpi_gather_data_extents(q_sf, data_extents) ! -------------
+    subroutine s_mpi_gather_data_extents(q_sf, data_extents)
 
         real(kind(0d0)), dimension(:, :, :), intent(IN) :: q_sf
 
@@ -1588,7 +1588,7 @@ contains
 
 #endif
 
-    end subroutine s_mpi_gather_data_extents ! -----------------------------
+    end subroutine s_mpi_gather_data_extents
 
     !>  This subroutine gathers the sub-domain flow variable data
         !!      from all of the processors and puts it back together for
@@ -1596,7 +1596,7 @@ contains
         !!      This is only done for 1D simulations.
         !!  @param q_sf Flow variable defined on a single computational sub-domain
         !!  @param q_root_sf Flow variable defined on the entire computational domain
-    subroutine s_mpi_defragment_1d_flow_variable(q_sf, q_root_sf) ! --------
+    subroutine s_mpi_defragment_1d_flow_variable(q_sf, q_root_sf)
 
         real(kind(0d0)), &
             dimension(0:m, 0:0, 0:0), &
@@ -1617,10 +1617,10 @@ contains
 
 #endif
 
-    end subroutine s_mpi_defragment_1d_flow_variable ! ---------------------
+    end subroutine s_mpi_defragment_1d_flow_variable
 
     !> Deallocation procedures for the module
-    subroutine s_finalize_mpi_proxy_module() ! ---------------------------
+    subroutine s_finalize_mpi_proxy_module
 
 #ifdef MFC_MPI
 
@@ -1639,6 +1639,6 @@ contains
 
 #endif
 
-    end subroutine s_finalize_mpi_proxy_module ! -------------------------
+    end subroutine s_finalize_mpi_proxy_module
 
 end module m_mpi_proxy
diff --git a/src/post_process/m_start_up.f90 b/src/post_process/m_start_up.f90
index 4fb03d24b..6ba76958c 100644
--- a/src/post_process/m_start_up.f90
+++ b/src/post_process/m_start_up.f90
@@ -41,7 +41,7 @@ module m_start_up
     !>  Reads the configuration file post_process.inp, in order
         !!      to populate parameters in module m_global_parameters.f90
         !!      with the user provided inputs
-    subroutine s_read_input_file() ! ---------------------------------------
+    subroutine s_read_input_file
 
         character(LEN=name_len) :: file_loc !<
             !! Generic string used to store the address of a particular file
@@ -104,13 +104,13 @@ subroutine s_read_input_file() ! ---------------------------------------
             call s_mpi_abort('File post_process.inp is missing. Exiting ...')
         end if
 
-    end subroutine s_read_input_file ! -------------------------------------
+    end subroutine s_read_input_file
 
     !>  Checking that the user inputs make sense, i.e. that the
         !!      individual choices are compatible with the code's options
         !!      and that the combination of these choices results into a
         !!      valid configuration for the post-process
-    subroutine s_check_input_file() ! --------------------------------------
+    subroutine s_check_input_file
 
         character(LEN=len_trim(case_dir)) :: file_loc !<
             !! Generic string used to store the address of a particular file
@@ -133,7 +133,7 @@ subroutine s_check_input_file() ! --------------------------------------
 
         call s_check_inputs()
 
-    end subroutine s_check_input_file ! ------------------------------------
+    end subroutine s_check_input_file
 
     subroutine s_perform_time_step(t_step)
 
diff --git a/src/pre_process/m_assign_variables.f90 b/src/pre_process/m_assign_variables.f90
index 17e7621f3..3440d9be7 100644
--- a/src/pre_process/m_assign_variables.f90
+++ b/src/pre_process/m_assign_variables.f90
@@ -62,7 +62,7 @@ end subroutine s_assign_patch_xxxxx_primitive_variables
 
 contains
 
-    subroutine s_initialize_assign_variables_module()
+    subroutine s_initialize_assign_variables_module
 
         allocate (alf_sum%sf(0:m, 0:n, 0:p))
 
@@ -158,7 +158,7 @@ subroutine s_assign_patch_mixture_primitive_variables(patch_id, j, k, l, &
         ! Updating the patch identities bookkeeping variable
         if (1d0 - eta < 1d-16) patch_id_fp(j, k, l) = patch_id
 
-    end subroutine s_assign_patch_mixture_primitive_variables ! ------------
+    end subroutine s_assign_patch_mixture_primitive_variables
 
     !Stable perturbation in pressure (Ando)
     !! @param j the x-dir node index
@@ -595,7 +595,7 @@ subroutine s_assign_patch_species_primitive_variables(patch_id, j, k, l, &
         ! Updating the patch identities bookkeeping variable
         if (1d0 - eta < 1d-16) patch_id_fp(j, k, l) = patch_id
 
-    end subroutine s_assign_patch_species_primitive_variables! ------------
+    end subroutine s_assign_patch_species_primitive_variables
 
     subroutine s_finalize_assign_variables_module
 
diff --git a/src/pre_process/m_check_ib_patches.fpp b/src/pre_process/m_check_ib_patches.fpp
index 69f6080a7..f3cbac476 100644
--- a/src/pre_process/m_check_ib_patches.fpp
+++ b/src/pre_process/m_check_ib_patches.fpp
@@ -73,7 +73,7 @@ contains
         !!      the circle patch have consistently been inputted by the
         !!      user.
         !!  @param patch_id Patch identifier
-    subroutine s_check_circle_ib_patch_geometry(patch_id) ! -------------------
+    subroutine s_check_circle_ib_patch_geometry(patch_id)
 
         integer, intent(in) :: patch_id
 
@@ -92,13 +92,13 @@ contains
 
         end if
 
-    end subroutine s_check_circle_ib_patch_geometry ! -------------------------
+    end subroutine s_check_circle_ib_patch_geometry
 
     !>  This subroutine verifies that the geometric parameters of
         !!      the airfoil patch have consistently been inputted by the
         !!      user.
         !!  @param patch_id Patch identifier
-    subroutine s_check_airfoil_ib_patch_geometry(patch_id) ! -------------------
+    subroutine s_check_airfoil_ib_patch_geometry(patch_id)
 
         integer, intent(in) :: patch_id
 
@@ -116,13 +116,13 @@ contains
 
         end if
 
-    end subroutine s_check_airfoil_ib_patch_geometry ! -------------------------
+    end subroutine s_check_airfoil_ib_patch_geometry
 
     !>  This subroutine verifies that the geometric parameters of
         !!      the 3d airfoil patch have consistently been inputted by the
         !!      user.
         !!  @param patch_id Patch identifier
-    subroutine s_check_3d_airfoil_ib_patch_geometry(patch_id) ! -------------------
+    subroutine s_check_3d_airfoil_ib_patch_geometry(patch_id)
 
         integer, intent(in) :: patch_id
 
@@ -141,13 +141,13 @@ contains
 
         end if
 
-    end subroutine s_check_3d_airfoil_ib_patch_geometry ! -------------------------
+    end subroutine s_check_3d_airfoil_ib_patch_geometry
 
     !>  This subroutine verifies that the geometric parameters of
         !!      the rectangle patch have consistently been inputted by
         !!      the user.
         !!  @param patch_id Patch identifier
-    subroutine s_check_rectangle_ib_patch_geometry(patch_id) ! ----------------
+    subroutine s_check_rectangle_ib_patch_geometry(patch_id)
 
         integer, intent(in) :: patch_id
 
@@ -170,13 +170,13 @@ contains
 
         end if
 
-    end subroutine s_check_rectangle_ib_patch_geometry ! ----------------------
+    end subroutine s_check_rectangle_ib_patch_geometry
 
     !>  This subroutine verifies that the geometric parameters of
         !!      the sphere patch have consistently been inputted by
         !!      the user.
         !!  @param patch_id Patch identifier
-    subroutine s_check_sphere_ib_patch_geometry(patch_id) ! ----------------
+    subroutine s_check_sphere_ib_patch_geometry(patch_id)
 
         integer, intent(in) :: patch_id
 
@@ -199,13 +199,13 @@ contains
 
         end if
 
-    end subroutine s_check_sphere_ib_patch_geometry ! ----------------------
+    end subroutine s_check_sphere_ib_patch_geometry
 
     !>  This subroutine verifies that the geometric parameters of
         !!      the cylinder patch have consistently been inputted by
         !!      the user.
         !!  @param patch_id Patch identifier
-    subroutine s_check_cylinder_ib_patch_geometry(patch_id) ! -----------------
+    subroutine s_check_cylinder_ib_patch_geometry(patch_id)
 
         integer, intent(in) :: patch_id
 
@@ -244,12 +244,12 @@ contains
 
         end if
 
-    end subroutine s_check_cylinder_ib_patch_geometry ! -----------------------
+    end subroutine s_check_cylinder_ib_patch_geometry
 
     !!>  This subroutine verifies that the geometric parameters of
         !!      the inactive patch remain unaltered by the user inputs.
         !!  @param patch_id Patch identifier
-    subroutine s_check_inactive_ib_patch_geometry(patch_id) ! -----------------
+    subroutine s_check_inactive_ib_patch_geometry(patch_id)
 
         integer, intent(in) :: patch_id
 
@@ -276,6 +276,6 @@ contains
 
         end if
 
-    end subroutine s_check_inactive_ib_patch_geometry ! -----------------------
+    end subroutine s_check_inactive_ib_patch_geometry
 
 end module m_check_ib_patches
diff --git a/src/pre_process/m_check_patches.fpp b/src/pre_process/m_check_patches.fpp
index 01c356e36..535c3208f 100644
--- a/src/pre_process/m_check_patches.fpp
+++ b/src/pre_process/m_check_patches.fpp
@@ -141,7 +141,7 @@ contains
         !!      the line segment patch have consistently been inputted by
         !!      the user.
         !!  @param patch_id Patch identifier
-    subroutine s_check_line_segment_patch_geometry(patch_id) ! -------------
+    subroutine s_check_line_segment_patch_geometry(patch_id)
 
         integer, intent(in) :: patch_id
         call s_int_to_str(patch_id, iStr)
@@ -159,13 +159,13 @@ contains
 
         end if
 
-    end subroutine s_check_line_segment_patch_geometry ! -------------------
+    end subroutine s_check_line_segment_patch_geometry
 
     !>  This subroutine verifies that the geometric parameters of
         !!      the circle patch have consistently been inputted by the
         !!      user.
         !!  @param patch_id Patch identifier
-    subroutine s_check_circle_patch_geometry(patch_id) ! -------------------
+    subroutine s_check_circle_patch_geometry(patch_id)
 
         integer, intent(in) :: patch_id
         call s_int_to_str(patch_id, iStr)
@@ -183,13 +183,13 @@ contains
 
         end if
 
-    end subroutine s_check_circle_patch_geometry ! -------------------------
+    end subroutine s_check_circle_patch_geometry
 
     !>  This subroutine verifies that the geometric parameters of
         !!      the rectangle patch have consistently been inputted by
         !!      the user.
         !!  @param patch_id Patch identifier
-    subroutine s_check_rectangle_patch_geometry(patch_id) ! ----------------
+    subroutine s_check_rectangle_patch_geometry(patch_id)
 
         integer, intent(in) :: patch_id
         call s_int_to_str(patch_id, iStr)
@@ -211,13 +211,13 @@ contains
 
         end if
 
-    end subroutine s_check_rectangle_patch_geometry ! ----------------------
+    end subroutine s_check_rectangle_patch_geometry
 
     !> This subroutine verifies that the geometric parameters of
         !!      the line sweep patch have consistently been inputted by
         !!      the user.
         !!  @param patch_id Patch identifier
-    subroutine s_check_line_sweep_patch_geometry(patch_id) ! ---------------
+    subroutine s_check_line_sweep_patch_geometry(patch_id)
 
         integer, intent(in) :: patch_id
         call s_int_to_str(patch_id, iStr)
@@ -241,13 +241,13 @@ contains
 
         end if
 
-    end subroutine s_check_line_sweep_patch_geometry ! ---------------------
+    end subroutine s_check_line_sweep_patch_geometry
 
     !>  This subroutine verifies that the geometric parameters of
         !!      the ellipse patch have consistently been inputted by
         !!      the user.
         !!  @param patch_id Patch identifier
-    subroutine s_check_ellipse_patch_geometry(patch_id) ! ------------------
+    subroutine s_check_ellipse_patch_geometry(patch_id)
 
         integer, intent(in) :: patch_id
         call s_int_to_str(patch_id, iStr)
@@ -271,13 +271,13 @@ contains
 
         end if
 
-    end subroutine s_check_ellipse_patch_geometry ! ------------------------
+    end subroutine s_check_ellipse_patch_geometry
 
     !>  This subroutine verifies that the geometric parameters of
         !!      the Taylor Green vortex patch have been entered by the user
         !!      consistently.
         !!  @param patch_id Patch identifier
-    subroutine s_check_2D_TaylorGreen_vortex_patch_geometry(patch_id) ! --------
+    subroutine s_check_2D_TaylorGreen_vortex_patch_geometry(patch_id)
 
         integer, intent(in) :: patch_id
         call s_int_to_str(patch_id, iStr)
@@ -301,13 +301,13 @@ contains
 
         end if
 
-    end subroutine s_check_2D_TaylorGreen_vortex_patch_geometry! --------------
+    end subroutine s_check_2D_TaylorGreen_vortex_patch_geometry
 
     !>  This subroutine verifies that the geometric parameters of
         !!      the analytical patch have consistently been inputted by
         !!      the user.
         !!  @param patch_id Patch identifier
-    subroutine s_check_1D_analytical_patch_geometry(patch_id) ! ---------------
+    subroutine s_check_1D_analytical_patch_geometry(patch_id)
 
         integer, intent(in) :: patch_id
         call s_int_to_str(patch_id, iStr)
@@ -325,13 +325,13 @@ contains
                              'geometric parameters of 1D analytical '// &
                              'patch '//trim(iStr)//'. Exiting...')
         end if
-    end subroutine s_check_1D_analytical_patch_geometry ! ---------------------
+    end subroutine s_check_1D_analytical_patch_geometry
 
     !>  This subroutine verifies that the geometric parameters of
         !!      the analytical patch have consistently been inputted by
         !!      the user.
         !!  @param patch_id Patch identifier
-    subroutine s_check_2D_analytical_patch_geometry(patch_id) ! ---------------
+    subroutine s_check_2D_analytical_patch_geometry(patch_id)
 
         integer, intent(in) :: patch_id
         call s_int_to_str(patch_id, iStr)
@@ -351,13 +351,13 @@ contains
                              'geometric parameters of 2D analytical '// &
                              'patch '//trim(iStr)//'. Exiting...')
         end if
-    end subroutine s_check_2D_analytical_patch_geometry ! ---------------------
+    end subroutine s_check_2D_analytical_patch_geometry
 
     !>  This subroutine verifies that the geometric parameters of
         !!      the analytical patch have consistently been inputted by
         !!      the user.
         !!  @param patch_id Patch identifier
-    subroutine s_check_3D_analytical_patch_geometry(patch_id) ! ---------------
+    subroutine s_check_3D_analytical_patch_geometry(patch_id)
 
         integer, intent(in) :: patch_id
         call s_int_to_str(patch_id, iStr)
@@ -381,13 +381,13 @@ contains
                              'geometric parameters of 3D analytical '// &
                              'patch '//trim(iStr)//'. Exiting...')
         end if
-    end subroutine s_check_3D_analytical_patch_geometry ! ---------------------
+    end subroutine s_check_3D_analytical_patch_geometry
 
     !> This subroutine verifies that the geometric parameters of
         !!      the sphere patch have consistently been inputted by the
         !!      user.
         !!  @param patch_id Patch identifier
-    subroutine s_check_sphere_patch_geometry(patch_id) ! -------------------
+    subroutine s_check_sphere_patch_geometry(patch_id)
 
         integer, intent(in) :: patch_id
         call s_int_to_str(patch_id, iStr)
@@ -409,14 +409,13 @@ contains
 
         end if
 
-    end subroutine s_check_sphere_patch_geometry ! -------------------------
+    end subroutine s_check_sphere_patch_geometry
 
     !>  This subroutine verifies that the geometric parameters of
         !!      the spherical harmonic patch have consistently been
         !!      inputted by the user.
         !!  @param patch_id Patch identifier
-    subroutine s_check_spherical_harmonic_patch_geometry(patch_id) ! -------
-
+    subroutine s_check_spherical_harmonic_patch_geometry(patch_id)
         integer, intent(in) :: patch_id
 
         call s_int_to_str(patch_id, iStr)
@@ -444,13 +443,13 @@ contains
 
         end if
 
-    end subroutine s_check_spherical_harmonic_patch_geometry ! -------------
+    end subroutine s_check_spherical_harmonic_patch_geometry
 
     !>  This subroutine verifies that the geometric parameters of
         !!      the cuboid patch have consistently been inputted by the
         !!      user.
         !!  @param patch_id Patch identifier
-    subroutine s_check_cuboid_patch_geometry(patch_id) ! -------------------
+    subroutine s_check_cuboid_patch_geometry(patch_id)
 
         ! Patch identifier
         integer, intent(in) :: patch_id
@@ -477,13 +476,13 @@ contains
 
         end if
 
-    end subroutine s_check_cuboid_patch_geometry ! -------------------------
+    end subroutine s_check_cuboid_patch_geometry
 
     !>  This subroutine verifies that the geometric parameters of
         !!      the cylinder patch have consistently been inputted by the
         !!      user.
         !!  @param patch_id Patch identifier
-    subroutine s_check_cylinder_patch_geometry(patch_id) ! -----------------
+    subroutine s_check_cylinder_patch_geometry(patch_id)
 
         ! Patch identifier
         integer, intent(in) :: patch_id
@@ -522,13 +521,13 @@ contains
 
         end if
 
-    end subroutine s_check_cylinder_patch_geometry ! -----------------------
+    end subroutine s_check_cylinder_patch_geometry
 
     !>  This subroutine verifies that the geometric parameters of
         !!      the plane sweep patch have consistently been inputted by
         !!      the user.
         !!  @param patch_id Patch identifier
-    subroutine s_check_plane_sweep_patch_geometry(patch_id) ! --------------
+    subroutine s_check_plane_sweep_patch_geometry(patch_id)
 
         ! Patch identifier
         integer, intent(in) :: patch_id
@@ -555,13 +554,13 @@ contains
 
         end if
 
-    end subroutine s_check_plane_sweep_patch_geometry ! --------------------
+    end subroutine s_check_plane_sweep_patch_geometry
 
     !> This subroutine verifies that the geometric parameters of
         !!      the ellipsoid patch have consistently been inputted by
         !!      the user.
         !!  @param patch_id Patch identifier
-    subroutine s_check_ellipsoid_patch_geometry(patch_id) ! ----------------
+    subroutine s_check_ellipsoid_patch_geometry(patch_id)
 
         integer, intent(in) :: patch_id
         call s_int_to_str(patch_id, iStr)
@@ -587,12 +586,12 @@ contains
 
         end if
 
-    end subroutine s_check_ellipsoid_patch_geometry ! ----------------------
+    end subroutine s_check_ellipsoid_patch_geometry
 
     !!>  This subroutine verifies that the geometric parameters of
         !!      the inactive patch remain unaltered by the user inputs.
         !!  @param patch_id Patch identifier
-    subroutine s_check_inactive_patch_geometry(patch_id) ! -----------------
+    subroutine s_check_inactive_patch_geometry(patch_id)
 
         integer, intent(in) :: patch_id
         call s_int_to_str(patch_id, iStr)
@@ -634,13 +633,13 @@ contains
 
         end if
 
-    end subroutine s_check_inactive_patch_geometry ! -----------------------
+    end subroutine s_check_inactive_patch_geometry
 
     !>  This subroutine verifies that any rights granted to the
         !!      given active patch, to overwrite the preceding active
         !!      patches, were consistently inputted by the user.
         !!  @param patch_id Patch identifier
-    subroutine s_check_active_patch_alteration_rights(patch_id) ! ----------
+    subroutine s_check_active_patch_alteration_rights(patch_id)
 
         integer, intent(in) :: patch_id
         call s_int_to_str(patch_id, iStr)
@@ -656,13 +655,13 @@ contains
 
         end if
 
-    end subroutine s_check_active_patch_alteration_rights ! ----------------
+    end subroutine s_check_active_patch_alteration_rights
 
     !>  This subroutine verifies that the rights of the given
         !!      inactive patch, to overwrite the preceding patches,
         !!      remain unaltered by the user inputs.
         !!  @param patch_id Patch identifier
-    subroutine s_check_inactive_patch_alteration_rights(patch_id) ! --------
+    subroutine s_check_inactive_patch_alteration_rights(patch_id)
 
         ! Patch identifier
         integer, intent(in) :: patch_id
@@ -679,13 +678,13 @@ contains
 
         end if
 
-    end subroutine s_check_inactive_patch_alteration_rights ! --------------
+    end subroutine s_check_inactive_patch_alteration_rights
 
     !> This subroutine verifies that the smoothing parameters of
         !!      the given patch, which supports the smoothing out of its
         !!      boundaries, have consistently been inputted by the user.
         !!  @param patch_id Patch identifier
-    subroutine s_check_supported_patch_smoothing(patch_id) ! ---------------
+    subroutine s_check_supported_patch_smoothing(patch_id)
 
         integer, intent(in) :: patch_id
         call s_int_to_str(patch_id, iStr)
@@ -711,13 +710,13 @@ contains
 
         end if
 
-    end subroutine s_check_supported_patch_smoothing ! ---------------------
+    end subroutine s_check_supported_patch_smoothing
 
     !> This subroutine verifies that the smoothing parameters of
         !!      the given patch, which does not support the smoothing out
         !!          of its boundaries, remain unaltered by the user inputs.
         !!  @param patch_id Patch identifier
-    subroutine s_check_unsupported_patch_smoothing(patch_id) ! -------------
+    subroutine s_check_unsupported_patch_smoothing(patch_id)
 
         ! Patch identifier
         integer, intent(in) :: patch_id
@@ -736,13 +735,13 @@ contains
 
         end if
 
-    end subroutine s_check_unsupported_patch_smoothing ! -------------------
+    end subroutine s_check_unsupported_patch_smoothing
 
     !>  This subroutine verifies that the primitive variables
         !!      associated with the given active patch are physically
         !!      consistent.
         !!  @param patch_id Patch identifier
-    subroutine s_check_active_patch_primitive_variables(patch_id) ! --------
+    subroutine s_check_active_patch_primitive_variables(patch_id)
 
         integer, intent(in) :: patch_id
         call s_int_to_str(patch_id, iStr)
@@ -797,13 +796,13 @@ contains
 
         end if
 
-    end subroutine s_check_active_patch_primitive_variables ! --------------
+    end subroutine s_check_active_patch_primitive_variables
 
     !>  This subroutine verifies that the primitive variables
         !!      associated with the given inactive patch remain unaltered
         !!      by the user inputs.
         !!  @param patch_id Patch identifier
-    subroutine s_check_inactive_patch_primitive_variables(patch_id) ! ------
+    subroutine s_check_inactive_patch_primitive_variables(patch_id)
 
         integer, intent(in) :: patch_id
         call s_int_to_str(patch_id, iStr)
@@ -829,9 +828,9 @@ contains
 
         end if
 
-    end subroutine s_check_inactive_patch_primitive_variables ! ------------
+    end subroutine s_check_inactive_patch_primitive_variables
 
-    subroutine s_check_model_geometry(patch_id) ! ------------------------------
+    subroutine s_check_model_geometry(patch_id)
 
         integer, intent(in) :: patch_id
 
@@ -848,6 +847,6 @@ contains
 
         end if
 
-    end subroutine s_check_model_geometry ! -----------------------------------
+    end subroutine s_check_model_geometry
 
 end module m_check_patches
diff --git a/src/pre_process/m_data_output.fpp b/src/pre_process/m_data_output.fpp
index 6b38acfe5..94d5d8fb6 100644
--- a/src/pre_process/m_data_output.fpp
+++ b/src/pre_process/m_data_output.fpp
@@ -60,7 +60,7 @@ module m_data_output
             type(integer_field), &
                 intent(in) :: ib_markers
 
-        end subroutine s_write_abstract_data_files ! -------------------
+        end subroutine s_write_abstract_data_files
     end interface ! ========================================================
 
     character(LEN=path_len + 2*name_len), private :: t_step_dir !<
@@ -77,7 +77,7 @@ contains
         !!  time-step directory in the local processor rank folder
         !! @param q_cons_vf Conservative variables
         !! @param ib_markers track if a cell is within the immersed boundary
-    subroutine s_write_serial_data_files(q_cons_vf, ib_markers) ! -----------
+    subroutine s_write_serial_data_files(q_cons_vf, ib_markers)
         type(scalar_field), &
             dimension(sys_size), &
             intent(in) :: q_cons_vf
@@ -464,13 +464,13 @@ contains
             end do
         end if
 
-    end subroutine s_write_serial_data_files ! ------------------------------------
+    end subroutine s_write_serial_data_files
 
     !> Writes grid and initial condition data files in parallel to the "0"
         !!  time-step directory in the local processor rank folder
         !! @param q_cons_vf Conservative variables
         !! @param ib_markers track if a cell is within the immersed boundary
-    subroutine s_write_parallel_data_files(q_cons_vf, ib_markers) ! --
+    subroutine s_write_parallel_data_files(q_cons_vf, ib_markers)
 
         ! Conservative variables
         type(scalar_field), &
@@ -714,11 +714,11 @@ contains
         end if
 #endif
 
-    end subroutine s_write_parallel_data_files ! ---------------------------
+    end subroutine s_write_parallel_data_files
 
     !> Computation of parameters, allocation procedures, and/or
         !!              any other tasks needed to properly setup the module
-    subroutine s_initialize_data_output_module ! ----------------------------
+    subroutine s_initialize_data_output_module
         ! Generic string used to store the address of a particular file
         character(LEN=len_trim(case_dir) + 2*name_len) :: file_loc
 
@@ -779,13 +779,13 @@ contains
 
         close (1)
 
-    end subroutine s_initialize_data_output_module ! --------------------------
+    end subroutine s_initialize_data_output_module
 
     !> Resets s_write_data_files pointer
-    subroutine s_finalize_data_output_module ! ---------------------------
+    subroutine s_finalize_data_output_module
 
         s_write_data_files => null()
 
-    end subroutine s_finalize_data_output_module ! -------------------------
+    end subroutine s_finalize_data_output_module
 
 end module m_data_output
diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp
index c2d7b157c..ee5dfeba1 100644
--- a/src/pre_process/m_global_parameters.fpp
+++ b/src/pre_process/m_global_parameters.fpp
@@ -228,7 +228,7 @@ contains
     !>  Assigns default values to user inputs prior to reading
         !!              them in. This allows for an easier consistency check of
         !!              these parameters once they are read from the input file.
-    subroutine s_assign_default_values_to_user_inputs ! ------------------
+    subroutine s_assign_default_values_to_user_inputs
 
         integer :: i !< Generic loop operator
 
@@ -427,11 +427,11 @@ contains
             fluid_pp(i)%G = 0d0
         end do
 
-    end subroutine s_assign_default_values_to_user_inputs ! ----------------
+    end subroutine s_assign_default_values_to_user_inputs
 
     !> Computation of parameters, allocation procedures, and/or
         !! any other tasks needed to properly setup the module
-    subroutine s_initialize_global_parameters_module ! ----------------------
+    subroutine s_initialize_global_parameters_module
 
         integer :: i, j, fac
 
@@ -730,9 +730,9 @@ contains
 
         allocate (logic_grid(0:m, 0:n, 0:p))
 
-    end subroutine s_initialize_global_parameters_module ! --------------------
+    end subroutine s_initialize_global_parameters_module
 
-    subroutine s_initialize_parallel_io ! --------------------------------
+    subroutine s_initialize_parallel_io
 
         num_dims = 1 + min(1, n) + min(1, p)
 
@@ -757,9 +757,9 @@ contains
 
 #endif
 
-    end subroutine s_initialize_parallel_io ! ------------------------------
+    end subroutine s_initialize_parallel_io
 
-    subroutine s_finalize_global_parameters_module() ! ------------------------
+    subroutine s_finalize_global_parameters_module
 
         integer :: i
 
@@ -791,6 +791,6 @@ contains
 
 #endif
 
-    end subroutine s_finalize_global_parameters_module ! ----------------------
+    end subroutine s_finalize_global_parameters_module
 
 end module m_global_parameters
diff --git a/src/pre_process/m_grid.f90 b/src/pre_process/m_grid.f90
index 3b812bd2c..9ec879fa9 100644
--- a/src/pre_process/m_grid.f90
+++ b/src/pre_process/m_grid.f90
@@ -37,11 +37,11 @@ module m_grid
 
     abstract interface ! ===================================================
 
-        subroutine s_generate_abstract_grid ! ------------------------
+        subroutine s_generate_abstract_grid
 
             ! integer, intent(IN), optional :: dummy
 
-        end subroutine s_generate_abstract_grid ! ----------------------
+        end subroutine s_generate_abstract_grid
 
     end interface ! ========================================================
 
@@ -54,7 +54,7 @@ end subroutine s_generate_abstract_grid ! ----------------------
         !!              inputted by the user. The grid information is stored in
         !!              the grid variables containing coordinates of the cell-
         !!              centers and cell-boundaries.
-    subroutine s_generate_serial_grid ! -----------------------------------------
+    subroutine s_generate_serial_grid
 
         ! Generic loop iterator
         integer :: i, j             !< generic loop operatorss
@@ -189,14 +189,14 @@ subroutine s_generate_serial_grid ! -----------------------------------------
         end if
         ! ==================================================================
 
-    end subroutine s_generate_serial_grid ! ---------------------------------------
+    end subroutine s_generate_serial_grid
 
     !> The following subroutine generates either a uniform or
         !!              non-uniform rectilinear grid in parallel, defined by the parameters
         !!              inputted by the user. The grid information is stored in
         !!              the grid variables containing coordinates of the cell-
         !!              centers and cell-boundaries.
-    subroutine s_generate_parallel_grid !-------------------------
+    subroutine s_generate_parallel_grid
 
 #ifdef MFC_MPI
 
@@ -341,11 +341,11 @@ subroutine s_generate_parallel_grid !-------------------------
 
 #endif
 
-    end subroutine s_generate_parallel_grid ! ------------------------------
+    end subroutine s_generate_parallel_grid
 
     !> Computation of parameters, allocation procedures, and/or
         !!              any other tasks needed to properly setup the module
-    subroutine s_initialize_grid_module ! -----------------------------------
+    subroutine s_initialize_grid_module
 
         if (parallel_io .neqv. .true.) then
             s_generate_grid => s_generate_serial_grid
@@ -353,13 +353,13 @@ subroutine s_initialize_grid_module ! -----------------------------------
             s_generate_grid => s_generate_parallel_grid
         end if
 
-    end subroutine s_initialize_grid_module ! ---------------------------------
+    end subroutine s_initialize_grid_module
 
     !> Deallocation procedures for the module
-    subroutine s_finalize_grid_module ! --------------------------------
+    subroutine s_finalize_grid_module
 
         s_generate_grid => null()
 
-    end subroutine s_finalize_grid_module ! ------------------------------
+    end subroutine s_finalize_grid_module
 
 end module m_grid
diff --git a/src/pre_process/m_initial_condition.fpp b/src/pre_process/m_initial_condition.fpp
index 95e3a2ec0..5d58d46eb 100644
--- a/src/pre_process/m_initial_condition.fpp
+++ b/src/pre_process/m_initial_condition.fpp
@@ -62,7 +62,7 @@ contains
 
     !> Computation of parameters, allocation procedures, and/or
         !!              any other tasks needed to properly setup the module
-    subroutine s_initialize_initial_condition_module ! -------------------
+    subroutine s_initialize_initial_condition_module
 
         integer :: i !< generic loop iterator
 
@@ -107,7 +107,7 @@ contains
         patch_id_fp = 0
         ib_markers%sf = 0
 
-    end subroutine s_initialize_initial_condition_module ! -----------------
+    end subroutine s_initialize_initial_condition_module
 
     !>  This subroutine peruses the patches and depending on the
         !!              type of geometry associated with a particular patch, it
@@ -115,7 +115,7 @@ contains
         !!              on the grid using the primitive variables included with
         !!              the patch parameters. The subroutine is complete once the
         !!              primitive variables are converted to conservative ones.
-    subroutine s_generate_initial_condition ! ----------------------------
+    subroutine s_generate_initial_condition
 
         integer :: i  !< Generic loop operator
 
@@ -320,9 +320,9 @@ contains
             call s_initialize_pb(q_cons_vf, mv%sf, pb%sf)
         end if
 
-    end subroutine s_generate_initial_condition ! --------------------------
+    end subroutine s_generate_initial_condition
 
-    subroutine s_perturb_sphere ! ----------------------------------------
+    subroutine s_perturb_sphere
 
         integer :: i, j, k, l !< generic loop operators
 
@@ -352,9 +352,9 @@ contains
             end do
         end do
 
-    end subroutine s_perturb_sphere ! --------------------------------------
+    end subroutine s_perturb_sphere
 
-    subroutine s_perturb_surrounding_flow ! ------------------------------
+    subroutine s_perturb_surrounding_flow
 
         integer :: i, j, k, l !<  generic loop iterators
 
@@ -386,7 +386,7 @@ contains
             end do
         end do
 
-    end subroutine s_perturb_surrounding_flow ! ----------------------------
+    end subroutine s_perturb_surrounding_flow
 
     !>  This subroutine computes velocity perturbations for a temporal mixing
         !!              layer with hypertangent mean streamwise velocity profile
@@ -394,7 +394,7 @@ contains
         !!              instability waves with spatial wavenumbers, (4,0), (2,0),
         !!              and (1,0) are superposed. For a 3D waves, (4,4), (4,-4),
         !!              (2,2), (2,-2), (1,1), (1,-1) areadded on top of 2D waves.
-    subroutine s_superposition_instability_wave ! ------------------------
+    subroutine s_superposition_instability_wave
         real(kind(0d0)), dimension(5, 0:m, 0:n, 0:p) :: wave, wave1, wave2, wave_tmp
         real(kind(0d0)) :: tr, ti
         real(kind(0d0)) :: Lx, Lz
@@ -450,7 +450,7 @@ contains
             end do
         end do
 
-    end subroutine s_superposition_instability_wave ! ----------------------
+    end subroutine s_superposition_instability_wave
 
     !>  This subroutine computes instability waves for a given set of spatial
         !!              wavenumbers (alpha, beta) in x and z directions.
@@ -638,7 +638,7 @@ contains
     end subroutine s_generate_wave
 
     !>  Deallocation procedures for the module
-    subroutine s_finalize_initial_condition_module ! ---------------------
+    subroutine s_finalize_initial_condition_module
 
         integer :: i !< Generic loop iterator
 
@@ -655,6 +655,6 @@ contains
         deallocate (patch_id_fp)
         deallocate (ib_markers%sf)
 
-    end subroutine s_finalize_initial_condition_module ! -------------------
+    end subroutine s_finalize_initial_condition_module
 
 end module m_initial_condition
diff --git a/src/pre_process/m_mpi_proxy.fpp b/src/pre_process/m_mpi_proxy.fpp
index 529db75ea..23903aee1 100644
--- a/src/pre_process/m_mpi_proxy.fpp
+++ b/src/pre_process/m_mpi_proxy.fpp
@@ -33,7 +33,7 @@ contains
         !!       these are not available to the remaining processors. This
         !!       subroutine is then in charge of broadcasting the required
         !!       information.
-    subroutine s_mpi_bcast_user_inputs ! ---------------------------------
+    subroutine s_mpi_bcast_user_inputs
 
 #ifdef MFC_MPI
 
@@ -112,14 +112,14 @@ contains
         end do
 #endif
 
-    end subroutine s_mpi_bcast_user_inputs ! -------------------------------
+    end subroutine s_mpi_bcast_user_inputs
 
     !> Description: This subroutine takes care of efficiently distributing
         !!              the computational domain among the available processors
         !!             as well as recomputing some of the global parameters so
         !!              that they reflect the configuration of sub-domain that is
         !!              overseen by the local processor.
-    subroutine s_mpi_decompose_computational_domain ! --------------------
+    subroutine s_mpi_decompose_computational_domain
 
 #ifdef MFC_MPI
 
@@ -529,6 +529,6 @@ contains
 
 #endif
 
-    end subroutine s_mpi_decompose_computational_domain ! ------------------
+    end subroutine s_mpi_decompose_computational_domain
 
 end module m_mpi_proxy
diff --git a/src/pre_process/m_patches.fpp b/src/pre_process/m_patches.fpp
index e429446d8..9a41fc6b4 100644
--- a/src/pre_process/m_patches.fpp
+++ b/src/pre_process/m_patches.fpp
@@ -91,7 +91,7 @@ contains
     !! @param patch_id patch identifier
     !! @param patch_id_fp Array to track patch ids
     !! @param q_prim_vf Array of primitive variables
-    subroutine s_line_segment(patch_id, patch_id_fp, q_prim_vf) ! ----------------------------------
+    subroutine s_line_segment(patch_id, patch_id_fp, q_prim_vf)
 
         integer, intent(in) :: patch_id
         integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp
@@ -140,7 +140,7 @@ contains
             end if
         end do
 
-    end subroutine s_line_segment ! ----------------------------------------
+    end subroutine s_line_segment
 
     !>  The spiral patch is a 2D geometry that may be used, The geometry
         !!              of the patch is well-defined when its centroid and radius
@@ -149,7 +149,7 @@ contains
         !! @param patch_id patch identifier
         !! @param patch_id_fp Array to track patch ids
         !! @param q_prim_vf Array of primitive variables
-    subroutine s_spiral(patch_id, patch_id_fp, q_prim_vf) ! ----------------------------------------
+    subroutine s_spiral(patch_id, patch_id_fp, q_prim_vf)
 
         integer, intent(in) :: patch_id
         integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp
@@ -204,7 +204,7 @@ contains
             end do
         end do
 
-    end subroutine s_spiral ! ----------------------------------------------
+    end subroutine s_spiral
 
     !> The circular patch is a 2D geometry that may be used, for
         !!              example, in creating a bubble or a droplet. The geometry
@@ -215,7 +215,7 @@ contains
         !! @param patch_id_fp Array to track patch ids
         !! @param q_prim_vf Array of primitive variables
         !! @param ib True if this patch is an immersed boundary
-    subroutine s_circle(patch_id, patch_id_fp, q_prim_vf, ib) ! ----------------------------------------
+    subroutine s_circle(patch_id, patch_id_fp, q_prim_vf, ib)
 
         integer, intent(in) :: patch_id
         integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp
@@ -289,7 +289,7 @@ contains
             end do
         end do
 
-    end subroutine s_circle ! ----------------------------------------------
+    end subroutine s_circle
 
     !! @param patch_id is the patch identifier
     !! @param patch_id_fp Array to track patch ids
@@ -629,7 +629,7 @@ contains
         !! @param patch_id is the patch identifier
         !! @param patch_id_fp Array to track patch ids
         !! @param q_prim_vf Array of primitive variables
-    subroutine s_varcircle(patch_id, patch_id_fp, q_prim_vf) ! ----------------------------------------
+    subroutine s_varcircle(patch_id, patch_id_fp, q_prim_vf)
 
         ! Patch identifier
         integer, intent(in) :: patch_id
@@ -682,12 +682,12 @@ contains
             end do
         end do
 
-    end subroutine s_varcircle ! ----------------------------------------------
+    end subroutine s_varcircle
 
     !! @param patch_id is the patch identifier
     !! @param patch_id_fp Array to track patch ids
     !! @param q_prim_vf Array of primitive variables
-    subroutine s_3dvarcircle(patch_id, patch_id_fp, q_prim_vf) ! ----------------------------------------
+    subroutine s_3dvarcircle(patch_id, patch_id_fp, q_prim_vf)
 
         ! Patch identifier
         integer, intent(in) :: patch_id
@@ -746,7 +746,7 @@ contains
             end do
         end do
 
-    end subroutine s_3dvarcircle ! ----------------------------------------------
+    end subroutine s_3dvarcircle
 
     !> The elliptical patch is a 2D geometry. The geometry of
         !!      the patch is well-defined when its centroid and radii
@@ -755,7 +755,7 @@ contains
         !! @param patch_id is the patch identifier
         !! @param patch_id_fp Array to track patch ids
         !! @param q_prim_vf Array of primitive variables
-    subroutine s_ellipse(patch_id, patch_id_fp, q_prim_vf) ! ---------------------------------------
+    subroutine s_ellipse(patch_id, patch_id_fp, q_prim_vf)
 
         integer, intent(in) :: patch_id
         integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp
@@ -812,7 +812,7 @@ contains
             end do
         end do
 
-    end subroutine s_ellipse ! ---------------------------------------------
+    end subroutine s_ellipse
 
     !> The ellipsoidal patch is a 3D geometry. The geometry of
         !!       the patch is well-defined when its centroid and radii
@@ -821,7 +821,7 @@ contains
         !! @param patch_id is the patch identifier
         !! @param patch_id_fp Array to track patch ids
         !! @param q_prim_vf Array of primitive variables
-    subroutine s_ellipsoid(patch_id, patch_id_fp, q_prim_vf) ! -------------------------------------
+    subroutine s_ellipsoid(patch_id, patch_id_fp, q_prim_vf)
 
         ! Patch identifier
         integer, intent(in) :: patch_id
@@ -893,7 +893,7 @@ contains
             end do
         end do
 
-    end subroutine s_ellipsoid ! -------------------------------------------
+    end subroutine s_ellipsoid
 
     !> The rectangular patch is a 2D geometry that may be used,
         !!              for example, in creating a solid boundary, or pre-/post-
@@ -907,7 +907,7 @@ contains
         !! @param patch_id_fp Array to track patch ids
         !! @param q_prim_vf Array of primitive variables
         !! @param ib True if this patch is an immersed boundary
-    subroutine s_rectangle(patch_id, patch_id_fp, q_prim_vf, ib) ! -------------------------------------
+    subroutine s_rectangle(patch_id, patch_id_fp, q_prim_vf, ib)
 
         integer, intent(in) :: patch_id
         integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp
@@ -996,7 +996,7 @@ contains
             end do
         end do
 
-    end subroutine s_rectangle ! -------------------------------------------
+    end subroutine s_rectangle
 
     !> The swept line patch is a 2D geometry that may be used,
         !!      for example, in creating a solid boundary, or pre-/post-
@@ -1008,7 +1008,7 @@ contains
         !! @param patch_id is the patch identifier
         !! @param patch_id_fp Array to track patch ids
         !! @param q_prim_vf Array of primitive variables
-    subroutine s_sweep_line(patch_id, patch_id_fp, q_prim_vf) ! ------------------------------------
+    subroutine s_sweep_line(patch_id, patch_id_fp, q_prim_vf)
 
         integer, intent(in) :: patch_id
         integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp
@@ -1064,7 +1064,7 @@ contains
             end do
         end do
 
-    end subroutine s_sweep_line ! ------------------------------------------
+    end subroutine s_sweep_line
 
     !> The Taylor Green vortex is 2D decaying vortex that may be used,
         !!              for example, to verify the effects of viscous attenuation.
@@ -1073,7 +1073,7 @@ contains
         !! @param patch_id is the patch identifier
         !! @param patch_id_fp Array to track patch ids
         !! @param q_prim_vf Array of primitive variables
-    subroutine s_2D_TaylorGreen_Vortex(patch_id, patch_id_fp, q_prim_vf) ! ----------------------------
+    subroutine s_2D_TaylorGreen_Vortex(patch_id, patch_id_fp, q_prim_vf)
 
         integer, intent(in) :: patch_id
         integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp
@@ -1142,14 +1142,14 @@ contains
             end do
         end do
 
-    end subroutine s_2D_TaylorGreen_Vortex ! -----------------------------------
+    end subroutine s_2D_TaylorGreen_Vortex
 
     !>  This patch assigns the primitive variables as analytical
         !!  functions such that the code can be verified.
         !! @param patch_id is the patch identifier
         !! @param patch_id_fp Array to track patch ids
         !! @param q_prim_vf Array of primitive variables
-    subroutine s_1D_analytical(patch_id, patch_id_fp, q_prim_vf) ! ---------------------------------
+    subroutine s_1D_analytical(patch_id, patch_id_fp, q_prim_vf)
 
         ! Patch identifier
         integer, intent(in) :: patch_id
@@ -1202,12 +1202,12 @@ contains
             end if
         end do
 
-    end subroutine s_1D_analytical ! ---------------------------------------
+    end subroutine s_1D_analytical
 
         !! @param patch_id is the patch identifier
         !! @param patch_id_fp Array to track patch ids
         !! @param q_prim_vf Array of primitive variables
-    subroutine s_1d_bubble_pulse(patch_id, patch_id_fp, q_prim_vf) ! ---------------------------------
+    subroutine s_1d_bubble_pulse(patch_id, patch_id_fp, q_prim_vf)
         ! Description: This patch assigns the primitive variables as analytical
         !       functions such that the code can be verified.
 
@@ -1257,14 +1257,14 @@ contains
             end if
         end do
 
-    end subroutine s_1D_bubble_pulse ! ---------------------------------------
+    end subroutine s_1D_bubble_pulse
 
     !>  This patch assigns the primitive variables as analytical
         !!  functions such that the code can be verified.
         !! @param patch_id is the patch identifier
         !! @param patch_id_fp Array to track patch ids
         !! @param q_prim_vf Array of primitive variables
-    subroutine s_2D_analytical(patch_id, patch_id_fp, q_prim_vf) ! ---------------------------------
+    subroutine s_2D_analytical(patch_id, patch_id_fp, q_prim_vf)
 
         integer, intent(in) :: patch_id
         integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp
@@ -1326,14 +1326,14 @@ contains
             end do
         end do
 
-    end subroutine s_2D_analytical ! ---------------------------------------
+    end subroutine s_2D_analytical
 
     !> This patch assigns the primitive variables as analytical
         !!      functions such that the code can be verified.
         !! @param patch_id is the patch identifier
         !! @param patch_id_fp Array to track patch ids
         !! @param q_prim_vf Array of primitive variables
-    subroutine s_3D_analytical(patch_id, patch_id_fp, q_prim_vf) ! ---------------------------------
+    subroutine s_3D_analytical(patch_id, patch_id_fp, q_prim_vf)
 
         integer, intent(in) :: patch_id
         integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp
@@ -1410,14 +1410,14 @@ contains
             end do
         end do
 
-    end subroutine s_3D_analytical ! ---------------------------------------
+    end subroutine s_3D_analytical
 
     !> This patch generates the shape of the spherical harmonics
         !!      as a perturbation to a perfect sphere
         !! @param patch_id is the patch identifier
         !! @param patch_id_fp Array to track patch ids
         !! @param q_prim_vf Array of primitive variables
-    subroutine s_spherical_harmonic(patch_id, patch_id_fp, q_prim_vf) ! ----------------------------
+    subroutine s_spherical_harmonic(patch_id, patch_id_fp, q_prim_vf)
 
         integer, intent(in) :: patch_id
         integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp
@@ -1539,7 +1539,7 @@ contains
             end do
         end do
 
-    end subroutine s_spherical_harmonic ! ----------------------------------
+    end subroutine s_spherical_harmonic
 
     !>          The spherical patch is a 3D geometry that may be used,
         !!              for example, in creating a bubble or a droplet. The patch
@@ -1550,7 +1550,7 @@ contains
         !! @param patch_id_fp Array to track patch ids
         !! @param q_prim_vf Array of primitive variables
         !! @param ib True if this patch is an immersed boundary
-    subroutine s_sphere(patch_id, patch_id_fp, q_prim_vf, ib) ! ----------------------------------------
+    subroutine s_sphere(patch_id, patch_id_fp, q_prim_vf, ib)
 
         integer, intent(in) :: patch_id
         integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp
@@ -1641,7 +1641,7 @@ contains
             end do
         end do
 
-    end subroutine s_sphere ! ----------------------------------------------
+    end subroutine s_sphere
 
     !> The cuboidal patch is a 3D geometry that may be used, for
         !!              example, in creating a solid boundary, or pre-/post-shock
@@ -1654,7 +1654,7 @@ contains
         !! @param patch_id is the patch identifier
         !! @param patch_id_fp Array to track patch ids
         !! @param q_prim_vf Array of primitive variables
-    subroutine s_cuboid(patch_id, patch_id_fp, q_prim_vf) ! ----------------------------------------
+    subroutine s_cuboid(patch_id, patch_id_fp, q_prim_vf)
 
         integer, intent(in) :: patch_id
         integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp
@@ -1723,7 +1723,7 @@ contains
             end do
         end do
 
-    end subroutine s_cuboid ! ----------------------------------------------
+    end subroutine s_cuboid
 
     !> The cylindrical patch is a 3D geometry that may be used,
         !!              for example, in setting up a cylindrical solid boundary
@@ -1737,7 +1737,7 @@ contains
         !! @param patch_id_fp Array to track patch ids
         !! @param q_prim_vf Array of primitive variables
         !! @param ib True if this patch is an immersed boundary
-    subroutine s_cylinder(patch_id, patch_id_fp, q_prim_vf, ib) ! --------------------------------------
+    subroutine s_cylinder(patch_id, patch_id_fp, q_prim_vf, ib)
 
         integer, intent(in) :: patch_id
         integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp
@@ -1881,7 +1881,7 @@ contains
             end do
         end do
 
-    end subroutine s_cylinder ! --------------------------------------------
+    end subroutine s_cylinder
 
     !>      The swept plane patch is a 3D geometry that may be used,
         !!              for example, in creating a solid boundary, or pre-/post-
@@ -1893,7 +1893,7 @@ contains
         !! @param patch_id is the patch identifier
         !! @param patch_id_fp Array to track patch ids
         !! @param q_prim_vf Primitive variables
-    subroutine s_sweep_plane(patch_id, patch_id_fp, q_prim_vf) ! -----------------------------------
+    subroutine s_sweep_plane(patch_id, patch_id_fp, q_prim_vf)
 
         integer, intent(in) :: patch_id
         integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp
@@ -1963,13 +1963,13 @@ contains
             end do
         end do
 
-    end subroutine s_sweep_plane ! -----------------------------------------
+    end subroutine s_sweep_plane
 
     !> The STL patch is a 2/3D geometry that is imported from an STL file.
     !! @param patch_id is the patch identifier
     !! @param patch_id_fp Array to track patch ids
     !! @param q_prim_vf Primitive variables
-    subroutine s_model(patch_id, patch_id_fp, q_prim_vf) ! ---------------------
+    subroutine s_model(patch_id, patch_id_fp, q_prim_vf)
 
         integer, intent(in) :: patch_id
         integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp
@@ -2074,7 +2074,7 @@ contains
 
         call s_model_free(model)
 
-    end subroutine s_model ! ---------------------------------------------------
+    end subroutine s_model
 
     subroutine s_convert_cylindrical_to_cartesian_coord(cyl_y, cyl_z)
         !$acc routine seq
@@ -2084,7 +2084,7 @@ contains
         cart_y = cyl_y*sin(cyl_z)
         cart_z = cyl_y*cos(cyl_z)
 
-    end subroutine s_convert_cylindrical_to_cartesian_coord ! --------------
+    end subroutine s_convert_cylindrical_to_cartesian_coord
 
     function f_convert_cyl_to_cart(cyl) result(cart)
 
@@ -2106,7 +2106,7 @@ contains
 
         sph_phi = atan(cyl_y/cyl_x)
 
-    end subroutine s_convert_cylindrical_to_spherical_coord ! --------------
+    end subroutine s_convert_cylindrical_to_spherical_coord
 
     !> Archimedes spiral function
     !! @param myth Angle
diff --git a/src/pre_process/m_start_up.fpp b/src/pre_process/m_start_up.fpp
index b761b402c..1d365aa77 100644
--- a/src/pre_process/m_start_up.fpp
+++ b/src/pre_process/m_start_up.fpp
@@ -68,13 +68,13 @@ module m_start_up
 
     abstract interface ! ===================================================
 
-        subroutine s_read_abstract_grid_data_files! ----------
+        subroutine s_read_abstract_grid_data_files
 
-        end subroutine s_read_abstract_grid_data_files ! ---------------
+        end subroutine s_read_abstract_grid_data_files
 
         !! @param q_cons_vf Conservative variables
         !! @param ib_markers track if a cell is within the immersed boundary
-        subroutine s_read_abstract_ic_data_files(q_cons_vf, ib_markers) ! -----------
+        subroutine s_read_abstract_ic_data_files(q_cons_vf, ib_markers)
 
             import :: scalar_field, integer_field, sys_size, pres_field
 
@@ -85,7 +85,7 @@ module m_start_up
             type(integer_field), &
                 intent(inout) :: ib_markers
 
-        end subroutine s_read_abstract_ic_data_files ! -----------------
+        end subroutine s_read_abstract_ic_data_files
 
     end interface ! ========================================================
 
@@ -104,7 +104,7 @@ contains
     !>  Reads the configuration file pre_process.inp, in order to
         !!      populate the parameters in module m_global_parameters.f90
         !!      with the user provided inputs
-    subroutine s_read_input_file ! ---------------------------------------
+    subroutine s_read_input_file
 
         character(LEN=name_len) :: file_loc  !<
             !! Generic string used to store the address of a particular file
@@ -167,13 +167,13 @@ contains
             call s_mpi_abort('File pre_process.inp is missing. Exiting ...')
         end if
 
-    end subroutine s_read_input_file ! -------------------------------------
+    end subroutine s_read_input_file
 
     !>  Checking that the user inputs make sense, i.e. that the
     !!      individual choices are compatible with the code's options
     !!      and that the combination of these choices results into a
     !!      valid configuration for the pre-process
-    subroutine s_check_input_file ! --------------------------------------
+    subroutine s_check_input_file
 
         character(LEN=len_trim(case_dir)) :: file_loc !<
             !! Generic string used to store the address of a particular file
@@ -202,12 +202,12 @@ contains
 
         if (ib) call s_check_ib_patches()
 
-    end subroutine s_check_input_file ! ------------------------------------
+    end subroutine s_check_input_file
 
     !> The goal of this subroutine is to read in any preexisting
         !!      grid data as well as based on the imported grid, complete
         !!      the necessary global computational domain parameters.
-    subroutine s_read_serial_grid_data_files ! ---
+    subroutine s_read_serial_grid_data_files
 
         ! Generic string used to store the address of a particular file
         character(LEN=len_trim(case_dir) + 3*name_len) :: file_loc
@@ -345,13 +345,13 @@ contains
             call s_create_directory(trim(proc_rank_dir)//'/0')
         end if
 
-    end subroutine s_read_serial_grid_data_files ! --------------------------------
+    end subroutine s_read_serial_grid_data_files
 
     !> Cell-boundary data are checked for consistency by looking
         !!      at the (non-)uniform cell-width distributions for all the
         !!      active coordinate directions and making sure that all of
         !!      the cell-widths are positively valued
-    subroutine s_check_grid_data_files ! -----------------
+    subroutine s_check_grid_data_files
 
         ! Cell-boundary Data Consistency Check in x-direction ==============
 
@@ -390,7 +390,7 @@ contains
 
         ! ==================================================================
 
-    end subroutine s_check_grid_data_files ! -------------------------------
+    end subroutine s_check_grid_data_files
 
     !> The goal of this subroutine is to read in any preexisting
         !!      initial condition data files so that they may be used by
@@ -398,7 +398,7 @@ contains
         !!      all new initial condition.
         !! @param q_cons_vf Conservative variables
         !! @param ib_markers track if a cell is within the immersed boundary
-    subroutine s_read_serial_ic_data_files(q_cons_vf, ib_markers) ! ---------------------------
+    subroutine s_read_serial_ic_data_files(q_cons_vf, ib_markers)
 
         type(scalar_field), &
             dimension(sys_size), &
@@ -525,13 +525,13 @@ contains
         call s_create_directory(trim(proc_rank_dir)//'/*')
         call s_create_directory(trim(proc_rank_dir)//'/0')
 
-    end subroutine s_read_serial_ic_data_files ! ----------------------------------
+    end subroutine s_read_serial_ic_data_files
 
     !> Cell-boundary data are checked for consistency by looking
         !!      at the (non-)uniform cell-width distributions for all the
         !!      active coordinate directions and making sure that all of
         !!      the cell-widths are positively valued
-    subroutine s_read_parallel_grid_data_files ! ----------------------------------
+    subroutine s_read_parallel_grid_data_files
 
 #ifdef MFC_MPI
 
@@ -628,7 +628,7 @@ contains
 
 #endif
 
-    end subroutine s_read_parallel_grid_data_files ! -----------------------
+    end subroutine s_read_parallel_grid_data_files
 
     !> The goal of this subroutine is to read in any preexisting
         !!      initial condition data files so that they may be used by
@@ -636,7 +636,7 @@ contains
         !!      all new initial condition.
         !! @param q_cons_vf Conservative variables
         !! @param ib_markers track if a cell is within the immersed boundary
-    subroutine s_read_parallel_ic_data_files(q_cons_vf, ib_markers) ! ------------------
+    subroutine s_read_parallel_ic_data_files(q_cons_vf, ib_markers)
 
         type(scalar_field), &
             dimension(sys_size), &
@@ -749,9 +749,9 @@ contains
 
 #endif
 
-    end subroutine s_read_parallel_ic_data_files ! -------------------------
+    end subroutine s_read_parallel_ic_data_files
 
-    subroutine s_initialize_modules ! ----------------------------------
+    subroutine s_initialize_modules
         ! Computation of parameters, allocation procedures, and/or any other tasks
         ! needed to properly setup the modules
         call s_initialize_global_parameters_module()
diff --git a/src/simulation/m_cbc.fpp b/src/simulation/m_cbc.fpp
index 521bf879c..f1c1bb113 100644
--- a/src/simulation/m_cbc.fpp
+++ b/src/simulation/m_cbc.fpp
@@ -149,7 +149,7 @@ contains
     !>  The computation of parameters, the allocation of memory,
         !!      the association of pointers and/or the execution of any
         !!      other procedures that are necessary to setup the module.
-    subroutine s_initialize_cbc_module() ! ---------------------------------
+    subroutine s_initialize_cbc_module
 
         integer :: i
         logical :: is_cbc
@@ -429,12 +429,12 @@ contains
             !$acc update device(bczb, bcze)
         end if
 
-    end subroutine s_initialize_cbc_module ! -------------------------------
+    end subroutine s_initialize_cbc_module
 
     !>  Compute CBC coefficients
         !!  @param cbc_dir_in CBC coordinate direction
         !!  @param cbc_loc_in CBC coordinate location
-    subroutine s_compute_cbc_coefficients(cbc_dir_in, cbc_loc_in) ! --------------
+    subroutine s_compute_cbc_coefficients(cbc_dir_in, cbc_loc_in)
         ! Description: The purpose of this subroutine is to compute the grid
         !              dependent FD and PI coefficients, or CBC coefficients,
         !              provided the CBC coordinate direction and location.
@@ -534,7 +534,7 @@ contains
 
         ! Nullifying CBC coefficients
 
-    end subroutine s_compute_cbc_coefficients ! ----------------------------
+    end subroutine s_compute_cbc_coefficients
 
     !!  The goal of the procedure is to associate the FD and PI
     !!      coefficients, or CBC coefficients, with the appropriate
@@ -542,7 +542,7 @@ contains
     !!      of the CBC.
     !!  @param cbc_dir_in CBC coordinate direction
     !!  @param cbc_loc_in CBC coordinate location
-    subroutine s_associate_cbc_coefficients_pointers(cbc_dir_in, cbc_loc_in) ! ---
+    subroutine s_associate_cbc_coefficients_pointers(cbc_dir_in, cbc_loc_in)
 
         integer, intent(IN) :: cbc_dir_in, cbc_loc_in
 
@@ -601,7 +601,7 @@ contains
 
         ! ==================================================================
 
-    end subroutine s_associate_cbc_coefficients_pointers ! -----------------
+    end subroutine s_associate_cbc_coefficients_pointers
 
     !>  The following is the implementation of the CBC based on
         !!      the work of Thompson (1987, 1990) on hyperbolic systems.
@@ -616,7 +616,7 @@ contains
         !!  @param ix Index bound in the first coordinate direction
         !!  @param iy Index bound in the second coordinate direction
         !!  @param iz Index bound in the third coordinate direction
-    subroutine s_cbc(q_prim_vf, flux_vf, flux_src_vf, & ! -----------------
+    subroutine s_cbc(q_prim_vf, flux_vf, flux_src_vf, &
                      cbc_dir_norm, cbc_loc_norm, &
                      ix, iy, iz)
 
@@ -865,21 +865,21 @@ contains
                         lambda(3) = vel(dir_idx(1)) + c
 
                         if ((cbc_loc == -1 .and. bc${XYZ}$b == -5) .or. (cbc_loc == 1 .and. bc${XYZ}$e == -5)) then
-                            call s_compute_slip_wall_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds) ! --------------
+                            call s_compute_slip_wall_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
                         else if ((cbc_loc == -1 .and. bc${XYZ}$b == -6) .or. (cbc_loc == 1 .and. bc${XYZ}$e == -6)) then
-                            call s_compute_nonreflecting_subsonic_buffer_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds) ! --------------
+                            call s_compute_nonreflecting_subsonic_buffer_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
                         else if ((cbc_loc == -1 .and. bc${XYZ}$b == -7) .or. (cbc_loc == 1 .and. bc${XYZ}$e == -7)) then
-                            call s_compute_nonreflecting_subsonic_inflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds) ! --------------
+                            call s_compute_nonreflecting_subsonic_inflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
                         else if ((cbc_loc == -1 .and. bc${XYZ}$b == -8) .or. (cbc_loc == 1 .and. bc${XYZ}$e == -8)) then
-                            call s_compute_nonreflecting_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds) ! --------------
+                            call s_compute_nonreflecting_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
                         else if ((cbc_loc == -1 .and. bc${XYZ}$b == -9) .or. (cbc_loc == 1 .and. bc${XYZ}$e == -9)) then
-                            call s_compute_force_free_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds) ! --------------
+                            call s_compute_force_free_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
                         else if ((cbc_loc == -1 .and. bc${XYZ}$b == -10) .or. (cbc_loc == 1 .and. bc${XYZ}$e == -10)) then
-                            call s_compute_constant_pressure_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds) ! --------------
+                            call s_compute_constant_pressure_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
                         else if ((cbc_loc == -1 .and. bc${XYZ}$b == -11) .or. (cbc_loc == 1 .and. bc${XYZ}$e == -11)) then
-                            call s_compute_supersonic_inflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds) ! --------------
+                            call s_compute_supersonic_inflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
                         else
-                            call s_compute_supersonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds) ! --------------
+                            call s_compute_supersonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
                         end if
 
                         ! Be careful about the cylindrical coordinate!
@@ -1007,7 +1007,7 @@ contains
         ! CBC coordinate direction.
         call s_finalize_cbc(flux_vf, flux_src_vf, &
                             ix, iy, iz)
-    end subroutine s_cbc ! -------------------------------------------------
+    end subroutine s_cbc
 
     !>  The computation of parameters, the allocation of memory,
         !!      the association of pointers and/or the execution of any
@@ -1019,7 +1019,7 @@ contains
         !!  @param ix Index bound in the first coordinate direction
         !!  @param iy Index bound in the second coordinate direction
         !!  @param iz Index bound in the third coordinate direction
-    subroutine s_initialize_cbc(q_prim_vf, flux_vf, flux_src_vf, & ! ------
+    subroutine s_initialize_cbc(q_prim_vf, flux_vf, flux_src_vf, &
                                 ix, iy, iz)
 
         type(scalar_field), &
@@ -1291,7 +1291,7 @@ contains
 
         ! ==================================================================
 
-    end subroutine s_initialize_cbc ! --------------------------------------
+    end subroutine s_initialize_cbc
 
     !>  Deallocation and/or the disassociation procedures that
         !!      are necessary in order to finalize the CBC application
@@ -1300,7 +1300,7 @@ contains
         !!  @param ix Index bound in the first coordinate direction
         !!  @param iy Index bound in the second coordinate direction
         !!  @param iz Index bound in the third coordinate direction
-    subroutine s_finalize_cbc(flux_vf, flux_src_vf, & ! -------------------
+    subroutine s_finalize_cbc(flux_vf, flux_src_vf, &
                               ix, iy, iz)
 
         type(scalar_field), &
@@ -1479,7 +1479,7 @@ contains
 
         ! Nullifying procedural pointer used in evaluation of L for the CBC
 
-    end subroutine s_finalize_cbc ! ----------------------------------------
+    end subroutine s_finalize_cbc
 
     ! Detext if the problem has any characteristic boundary conditions
     subroutine s_any_cbc_boundaries(toggle)
@@ -1497,7 +1497,7 @@ contains
     end subroutine
 
     !> Module deallocation and/or disassociation procedures
-    subroutine s_finalize_cbc_module() ! -----------------------------------
+    subroutine s_finalize_cbc_module
 
         logical :: is_cbc
 
@@ -1561,6 +1561,6 @@ contains
         ! Disassociating the pointer to the procedure that was utilized to
         ! to convert mixture or species variables to the mixture variables
 
-    end subroutine s_finalize_cbc_module ! ---------------------------------
+    end subroutine s_finalize_cbc_module
 
 end module m_cbc
diff --git a/src/simulation/m_compute_cbc.fpp b/src/simulation/m_compute_cbc.fpp
index d323c8b89..ed045104a 100644
--- a/src/simulation/m_compute_cbc.fpp
+++ b/src/simulation/m_compute_cbc.fpp
@@ -27,7 +27,7 @@ contains
         !!      Thompson (1990). At the slip wall (frictionless wall),
         !!      the normal component of velocity is zero at all times,
         !!      while the transverse velocities may be nonzero.
-    subroutine s_compute_slip_wall_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds) ! --------------
+    subroutine s_compute_slip_wall_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
 #ifdef CRAY_ACC_WAR
         !DIR$ INLINEALWAYS s_compute_slip_wall_L
 #else
@@ -49,13 +49,13 @@ contains
 
         L(advxe) = L(1)
 
-    end subroutine s_compute_slip_wall_L ! ---------------------------------
+    end subroutine s_compute_slip_wall_L
 
     !>  The L variables for the nonreflecting subsonic buffer CBC
         !!      see pg. 13 of Thompson (1987). The nonreflecting subsonic
         !!      buffer reduces the amplitude of any reflections caused by
         !!      outgoing waves.
-    subroutine s_compute_nonreflecting_subsonic_buffer_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds) ! --------------
+    subroutine s_compute_nonreflecting_subsonic_buffer_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
 #ifdef CRAY_ACC_WAR
         !DIR$ INLINEALWAYS s_compute_nonreflecting_subsonic_buffer_L
 #else
@@ -90,13 +90,12 @@ contains
         L(advxe) = (5d-1 - 5d-1*sign(1d0, lambda(3)))*lambda(3) &
                    *(dpres_ds + rho*c*dvel_ds(dir_idx(1)))
 
-    end subroutine s_compute_nonreflecting_subsonic_buffer_L ! -------------
-
+    end subroutine s_compute_nonreflecting_subsonic_buffer_L
     !>  The L variables for the nonreflecting subsonic inflow CBC
         !!      see pg. 455, Thompson (1990). This nonreflecting subsonic
         !!      CBC assumes an incoming flow and reduces the amplitude of
         !!      any reflections caused by outgoing waves.
-    subroutine s_compute_nonreflecting_subsonic_inflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds) ! --------------
+    subroutine s_compute_nonreflecting_subsonic_inflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
 #ifdef CRAY_ACC_WAR
         !DIR$ INLINEALWAYS ss_compute_nonreflecting_subsonic_inflow_L
 #else
@@ -116,13 +115,13 @@ contains
             L(i) = 0d0
         end do
 
-    end subroutine s_compute_nonreflecting_subsonic_inflow_L ! -------------
+    end subroutine s_compute_nonreflecting_subsonic_inflow_L
 
     !>  The L variables for the nonreflecting subsonic outflow
         !!      CBC see pg. 454 of Thompson (1990). This nonreflecting
         !!      subsonic CBC presumes an outgoing flow and reduces the
         !!      amplitude of any reflections caused by outgoing waves.
-    subroutine s_compute_nonreflecting_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds) ! --------------
+    subroutine s_compute_nonreflecting_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
 #ifdef CRAY_ACC_WAR
         !DIR$ INLINEALWAYS s_compute_nonreflecting_subsonic_outflow_L
 #else
@@ -153,7 +152,7 @@ contains
         ! bubble index
         L(advxe) = 0d0
 
-    end subroutine s_compute_nonreflecting_subsonic_outflow_L ! ------------
+    end subroutine s_compute_nonreflecting_subsonic_outflow_L
 
     !>  The L variables for the force-free subsonic outflow CBC,
         !!      see pg. 454 of Thompson (1990). The force-free subsonic
@@ -162,7 +161,7 @@ contains
         !!      direction to the boundary. As a result, a fluid element
         !!      at the boundary is simply advected outward at the fluid
         !!      velocity.
-    subroutine s_compute_force_free_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds) ! --------------
+    subroutine s_compute_force_free_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
 #ifdef CRAY_ACC_WAR
         !DIR$ INLINEALWAYS s_compute_force_free_subsonic_outflow_L
 #else
@@ -192,13 +191,13 @@ contains
 
         L(advxe) = L(1) + 2d0*rho*c*lambda(2)*dvel_ds(dir_idx(1))
 
-    end subroutine s_compute_force_free_subsonic_outflow_L ! ---------------
+    end subroutine s_compute_force_free_subsonic_outflow_L
 
     !>  L variables for the constant pressure subsonic outflow
         !!      CBC see pg. 455 Thompson (1990). The constant pressure
         !!      subsonic outflow maintains a fixed pressure at the CBC
         !!      boundary in absence of any transverse effects.
-    subroutine s_compute_constant_pressure_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds) ! --------------
+    subroutine s_compute_constant_pressure_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
 #ifdef CRAY_ACC_WAR
         !DIR$ INLINEALWAYS s_compute_constant_pressure_subsonic_outflow_L
 #else
@@ -228,14 +227,14 @@ contains
 
         L(advxe) = -L(1)
 
-    end subroutine s_compute_constant_pressure_subsonic_outflow_L ! --------
+    end subroutine s_compute_constant_pressure_subsonic_outflow_L
 
     !>  L variables for the supersonic inflow CBC, see pg. 453
         !!      Thompson (1990). The supersonic inflow CBC is a steady
         !!      state, or nearly a steady state, CBC in which only the
         !!      transverse terms may generate a time dependence at the
         !!      inflow boundary.
-    subroutine s_compute_supersonic_inflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds) ! --------------
+    subroutine s_compute_supersonic_inflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
 #ifdef CRAY_ACC_WAR
         !DIR$ INLINEALWAYS s_compute_supersonic_inflow_L
 #else
@@ -253,13 +252,13 @@ contains
             L(i) = 0d0
         end do
 
-    end subroutine s_compute_supersonic_inflow_L ! -------------------------
+    end subroutine s_compute_supersonic_inflow_L
 
     !>  L variables for the supersonic outflow CBC, see pg. 453
         !!      of Thompson (1990). For the supersonic outflow CBC, the
         !!      flow evolution at the boundary is determined completely
         !!      by the interior data.
-    subroutine s_compute_supersonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds) ! --------------
+    subroutine s_compute_supersonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
 #ifdef CRAY_ACC_WAR
         !DIR$ INLINEALWAYS s_compute_supersonic_outflow_L
 #else
@@ -289,6 +288,6 @@ contains
 
         L(advxe) = lambda(3)*(dpres_ds + rho*c*dvel_ds(dir_idx(1)))
 
-    end subroutine s_compute_supersonic_outflow_L ! ------------------------
+    end subroutine s_compute_supersonic_outflow_L
 
 end module m_compute_cbc
diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp
index 04258fe64..94336b6a7 100644
--- a/src/simulation/m_data_output.fpp
+++ b/src/simulation/m_data_output.fpp
@@ -65,7 +65,7 @@ module m_data_output
 
             integer, intent(IN) :: t_step
 
-        end subroutine s_write_abstract_data_files ! -------------------
+        end subroutine s_write_abstract_data_files
     end interface ! ========================================================
 #ifdef CRAY_ACC_WAR
     @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), icfl_sf)
@@ -107,7 +107,7 @@ contains
         !!      In general, this requires generating a table header for
         !!      those stability criteria which will be written at every
         !!      time-step.
-    subroutine s_open_run_time_information_file() ! ------------------------
+    subroutine s_open_run_time_information_file
 
         character(LEN=name_len) :: file_name = 'run_time.inf' !<
             !! Name of the run-time information file
@@ -165,11 +165,11 @@ contains
                 '============== ICFL Max ============='
         end if
 
-    end subroutine s_open_run_time_information_file ! ----------------------
+    end subroutine s_open_run_time_information_file
 
     !>  This opens a formatted data file where the root processor
         !!      can write out flow probe information
-    subroutine s_open_probe_files() ! --------------------------------------
+    subroutine s_open_probe_files
 
         character(LEN=path_len + 3*name_len) :: file_path !<
             !! Relative path to the probe data file in the case directory
@@ -214,7 +214,7 @@ contains
             end do
         end if
 
-    end subroutine s_open_probe_files ! ------------------------------------
+    end subroutine s_open_probe_files
 
     !>  The goal of the procedure is to output to the run-time
         !!      information file the stability criteria extrema in the
@@ -223,7 +223,7 @@ contains
         !!      these stability criteria extrema over all time-steps.
         !!  @param q_prim_vf Cell-average primitive variables
         !!  @param t_step Current time step
-    subroutine s_write_run_time_information(q_prim_vf, t_step) ! -----------
+    subroutine s_write_run_time_information(q_prim_vf, t_step)
 
         type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf
         integer, intent(IN) :: t_step
@@ -450,13 +450,13 @@ contains
 
         call s_mpi_barrier()
 
-    end subroutine s_write_run_time_information ! --------------------------
+    end subroutine s_write_run_time_information
 
     !>  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 t_step Current time-step
-    subroutine s_write_serial_data_files(q_cons_vf, q_prim_vf, t_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
@@ -821,13 +821,13 @@ contains
             end if
         end if
 
-    end subroutine s_write_serial_data_files ! ------------------------------------
+    end subroutine s_write_serial_data_files
 
     !>  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 t_step Current time-step
-    subroutine s_write_parallel_data_files(q_cons_vf, q_prim_vf, t_step) ! --
+    subroutine s_write_parallel_data_files(q_cons_vf, q_prim_vf, t_step)
 
         type(scalar_field), &
             dimension(sys_size), &
@@ -1006,13 +1006,13 @@ contains
 
 #endif
 
-    end subroutine s_write_parallel_data_files ! ---------------------------
+    end subroutine s_write_parallel_data_files
 
     !>  This writes a formatted data file for the flow probe information
         !!  @param t_step Current time-step
         !!  @param q_cons_vf Conservative variables
         !!  @param accel_mag Acceleration magnitude information
-    subroutine s_write_probe_files(t_step, q_cons_vf, accel_mag) ! -----------
+    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
@@ -1628,14 +1628,14 @@ contains
             end if
         end if
 
-    end subroutine s_write_probe_files ! -----------------------------------
+    end subroutine s_write_probe_files
 
     !>  The goal of this subroutine is to write to the run-time
         !!      information file basic footer information applicable to
         !!      the current computation and to close the file when done.
         !!      The footer contains the stability criteria extrema over
         !!      all of the time-steps and the simulation run-time.
-    subroutine s_close_run_time_information_file() ! -----------------------
+    subroutine s_close_run_time_information_file
 
         real(kind(0d0)) :: run_time !< Run-time of the simulation
 
@@ -1656,10 +1656,10 @@ contains
             '========================================'
         close (1)
 
-    end subroutine s_close_run_time_information_file ! ---------------------
+    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
 
@@ -1667,12 +1667,12 @@ contains
             close (i + 30)
         end do
 
-    end subroutine s_close_probe_files ! -----------------------------------
+    end subroutine s_close_probe_files
 
     !>  The computation of parameters, the allocation of memory,
         !!      the association of pointers and/or the execution of any
         !!      other procedures that are necessary to setup the module.
-    subroutine s_initialize_data_output_module() ! -------------------------
+    subroutine s_initialize_data_output_module
 
         type(int_bounds_info) :: ix, iy, iz
 
@@ -1710,10 +1710,10 @@ contains
             s_write_data_files => s_write_parallel_data_files
         end if
 
-    end subroutine s_initialize_data_output_module ! -----------------------
+    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
 
@@ -1728,6 +1728,6 @@ contains
         s_convert_to_mixture_variables => null()
         s_write_data_files => null()
 
-    end subroutine s_finalize_data_output_module ! -------------------------
+    end subroutine s_finalize_data_output_module
 
 end module m_data_output
diff --git a/src/simulation/m_derived_variables.f90 b/src/simulation/m_derived_variables.f90
index 53021f666..f27f68284 100644
--- a/src/simulation/m_derived_variables.f90
+++ b/src/simulation/m_derived_variables.f90
@@ -51,7 +51,7 @@ module m_derived_variables
 
     !>  Computation of parameters, allocation procedures, and/or
         !!      any other tasks needed to properly setup the module
-    subroutine s_initialize_derived_variables_module() ! ----------------------
+    subroutine s_initialize_derived_variables_module
 
         ! Allocating the variables which will store the coefficients of the
         ! centered family of finite-difference schemes. Note that sufficient
@@ -81,10 +81,10 @@ subroutine s_initialize_derived_variables_module() ! ----------------------
             end if
         end if
 
-    end subroutine s_initialize_derived_variables_module ! --------------------
+    end subroutine s_initialize_derived_variables_module
 
     !> Allocate and open derived variables. Computing FD coefficients.
-    subroutine s_initialize_derived_variables() ! -----------------------------
+    subroutine s_initialize_derived_variables
 
         if (probe_wrt) then
             ! Opening and writing header of flow probe files
@@ -108,11 +108,11 @@ subroutine s_initialize_derived_variables() ! -----------------------------
 
         end if
 
-    end subroutine s_initialize_derived_variables ! -----------------------------
+    end subroutine s_initialize_derived_variables
 
     !> Writes coherent body information, communication files, and probes.
         !!  @param t_step Current time-step
-    subroutine s_compute_derived_variables(t_step) ! -----------------------
+    subroutine s_compute_derived_variables(t_step)
 
         integer, intent(IN) :: t_step
 
@@ -161,7 +161,7 @@ subroutine s_compute_derived_variables(t_step) ! -----------------------
             call s_write_probe_files(t_step, q_cons_ts(1)%vf, accel_mag)
         end if
 
-    end subroutine s_compute_derived_variables ! ---------------------------
+    end subroutine s_compute_derived_variables
 
     !> This subroutine receives as inputs the indicator of the
         !!      component of the acceleration that should be outputted and
@@ -176,7 +176,7 @@ end subroutine s_compute_derived_variables ! ---------------------------
         !!  @param q_prim_vf3 Primitive variables
         !!  @param q_sf Acceleration component
     subroutine s_derive_acceleration_component(i, q_prim_vf0, q_prim_vf1, &
-                                               q_prim_vf2, q_prim_vf3, q_sf) ! ----------
+                                               q_prim_vf2, q_prim_vf3, q_sf)
         !DIR$ INLINEALWAYS s_derive_acceleration_component
         integer, intent(IN) :: i
 
@@ -314,10 +314,10 @@ subroutine s_derive_acceleration_component(i, q_prim_vf0, q_prim_vf1, &
             end do
         end if
 
-    end subroutine s_derive_acceleration_component ! --------------------------
+    end subroutine s_derive_acceleration_component
 
     !> Deallocation procedures for the module
-    subroutine s_finalize_derived_variables_module() ! -------------------
+    subroutine s_finalize_derived_variables_module
 
         ! Closing CoM and flow probe files
         if (proc_rank == 0) then
@@ -342,6 +342,6 @@ subroutine s_finalize_derived_variables_module() ! -------------------
         if (allocated(fd_coeff_y)) deallocate (fd_coeff_y)
         if (allocated(fd_coeff_z)) deallocate (fd_coeff_z)
 
-    end subroutine s_finalize_derived_variables_module ! -----------------
+    end subroutine s_finalize_derived_variables_module
 
 end module m_derived_variables
diff --git a/src/simulation/m_fftw.fpp b/src/simulation/m_fftw.fpp
index 603d65ce1..937ab1dab 100644
--- a/src/simulation/m_fftw.fpp
+++ b/src/simulation/m_fftw.fpp
@@ -80,7 +80,7 @@ contains
     !>  The purpose of this subroutine is to create the fftw plan
         !!      that will be used in the forward and backward DFTs when
         !!      applying the Fourier filter in the azimuthal direction.
-    subroutine s_initialize_fftw_module() ! ----------------------------------
+    subroutine s_initialize_fftw_module
 
         ! Size of input array going into DFT
         real_size = p + 1
@@ -131,13 +131,14 @@ contains
 
 #endif
 
-    end subroutine s_initialize_fftw_module ! ------------------------------
+    end subroutine s_initialize_fftw_module
 
     !>  The purpose of this subroutine is to apply a Fourier low-
         !!      pass filter to the flow variables in the azimuthal direction
         !!      to remove the high-frequency content. This alleviates the
         !!      restrictive CFL condition arising from cells near the axis.
-    subroutine s_apply_fourier_filter(q_cons_vf) ! --------------------------
+        !! @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
 
@@ -297,12 +298,12 @@ contains
         end do
 #endif
 
-    end subroutine s_apply_fourier_filter ! --------------------------------
+    end subroutine s_apply_fourier_filter
 
     !>  The purpose of this subroutine is to destroy the fftw plan
         !!      that will be used in the forward and backward DFTs when
         !!      applying the Fourier filter in the azimuthal direction.
-    subroutine s_finalize_fftw_module() ! ------------------------------------
+    subroutine s_finalize_fftw_module
 
 #if defined(MFC_OpenACC)
         @:DEALLOCATE_GLOBAL(data_real_gpu, data_fltr_cmplx_gpu, data_cmplx_gpu)
@@ -323,6 +324,6 @@ contains
         call fftw_destroy_plan(bwd_plan)
 #endif
 
-    end subroutine s_finalize_fftw_module ! --------------------------------
+    end subroutine s_finalize_fftw_module
 
 end module
diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp
index 39955308d..31fffd3f8 100644
--- a/src/simulation/m_global_parameters.fpp
+++ b/src/simulation/m_global_parameters.fpp
@@ -468,7 +468,7 @@ contains
     !> Assigns default values to the user inputs before reading
         !!  them in. This enables for an easier consistency check of
         !!  these parameters once they are read from the input file.
-    subroutine s_assign_default_values_to_user_inputs() ! ------------------
+    subroutine s_assign_default_values_to_user_inputs
 
         integer :: i, j !< Generic loop iterator
 
@@ -652,12 +652,12 @@ contains
             integral(i)%ymax = dflt_real
         end do
 
-    end subroutine s_assign_default_values_to_user_inputs ! ----------------
+    end subroutine s_assign_default_values_to_user_inputs
 
     !>  The computation of parameters, the allocation of memory,
         !!      the association of pointers and/or the execution of any
         !!      other procedures that are necessary to setup the module.
-    subroutine s_initialize_global_parameters_module() ! -------------------
+    subroutine s_initialize_global_parameters_module
 
         integer :: i, j, k
         integer :: fac
@@ -1070,10 +1070,10 @@ contains
         @:ALLOCATE_GLOBAL(z_cc(-buff_size:p + buff_size))
         @:ALLOCATE_GLOBAL(dz(-buff_size:p + buff_size))
 
-    end subroutine s_initialize_global_parameters_module ! -----------------
+    end subroutine s_initialize_global_parameters_module
 
     !> Initializes parallel infrastructure
-    subroutine s_initialize_parallel_io() ! --------------------------------
+    subroutine s_initialize_parallel_io
 
         #:if not MFC_CASE_OPTIMIZATION
             num_dims = 1 + min(1, n) + min(1, p)
@@ -1101,10 +1101,10 @@ contains
 
 #endif
 
-    end subroutine s_initialize_parallel_io ! ------------------------------
+    end subroutine s_initialize_parallel_io
 
     !> Module deallocation and/or disassociation procedures
-    subroutine s_finalize_global_parameters_module() ! ---------------------
+    subroutine s_finalize_global_parameters_module
 
         integer :: i
 
@@ -1137,6 +1137,6 @@ contains
         if (p == 0) return; 
         @:DEALLOCATE_GLOBAL(z_cb, z_cc, dz)
 
-    end subroutine s_finalize_global_parameters_module ! -------------------
+    end subroutine s_finalize_global_parameters_module
 
 end module m_global_parameters
diff --git a/src/simulation/m_hypoelastic.fpp b/src/simulation/m_hypoelastic.fpp
index 48951541e..d01ce4ffa 100644
--- a/src/simulation/m_hypoelastic.fpp
+++ b/src/simulation/m_hypoelastic.fpp
@@ -49,7 +49,7 @@ module m_hypoelastic
 
 contains
 
-    subroutine s_initialize_hypoelastic_module() ! --------------------
+    subroutine s_initialize_hypoelastic_module
 
         integer :: i
 
diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp
index 847743b1a..94272f15e 100644
--- a/src/simulation/m_mpi_proxy.fpp
+++ b/src/simulation/m_mpi_proxy.fpp
@@ -93,7 +93,7 @@ contains
     !> The computation of parameters, the allocation of memory,
         !!      the association of pointers and/or the execution of any
         !!      other procedures that are necessary to setup the module.
-    subroutine s_initialize_mpi_proxy_module() ! ---------------------------
+    subroutine s_initialize_mpi_proxy_module
 
 #ifdef MFC_MPI
 
@@ -164,14 +164,14 @@ contains
 
 #endif
 
-    end subroutine s_initialize_mpi_proxy_module ! -------------------------
+    end subroutine s_initialize_mpi_proxy_module
 
     !>  Since only the processor with rank 0 reads and verifies
         !!      the consistency of user inputs, these are initially not
         !!      available to the other processors. Then, the purpose of
         !!      this subroutine is to distribute the user inputs to the
         !!      remaining processors in the communicator.
-    subroutine s_mpi_bcast_user_inputs() ! ---------------------------------
+    subroutine s_mpi_bcast_user_inputs
 
 #ifdef MFC_MPI
 
@@ -260,7 +260,7 @@ contains
 
 #endif
 
-    end subroutine s_mpi_bcast_user_inputs ! -------------------------------
+    end subroutine s_mpi_bcast_user_inputs
 
     !>  The purpose of this procedure is to optimally decompose
         !!      the computational domain among the available processors.
@@ -268,7 +268,7 @@ contains
         !!      in each of the coordinate directions, approximately the
         !!      same number of cells, and then recomputing the affected
         !!      global parameters.
-    subroutine s_mpi_decompose_computational_domain() ! --------------------
+    subroutine s_mpi_decompose_computational_domain
 
 #ifdef MFC_MPI
 
@@ -636,7 +636,7 @@ contains
 
 #endif
 
-    end subroutine s_mpi_decompose_computational_domain ! ------------------
+    end subroutine s_mpi_decompose_computational_domain
 
     !>  The goal of this procedure is to populate the buffers of
         !!      the grid variables by communicating with the neighboring
@@ -646,7 +646,7 @@ contains
         !!      directly from those of the cell-width distributions.
         !!  @param mpi_dir MPI communication coordinate direction
         !!  @param pbc_loc Processor boundary condition (PBC) location
-    subroutine s_mpi_sendrecv_grid_variables_buffers(mpi_dir, pbc_loc) ! ---
+    subroutine s_mpi_sendrecv_grid_variables_buffers(mpi_dir, pbc_loc)
 
         integer, intent(IN) :: mpi_dir
         integer, intent(IN) :: pbc_loc
@@ -821,7 +821,7 @@ contains
 
 #endif
 
-    end subroutine s_mpi_sendrecv_grid_variables_buffers ! -----------------
+    end subroutine s_mpi_sendrecv_grid_variables_buffers
 
     !>  The goal of this procedure is to populate the buffers of
         !!      the cell-average conservative variables by communicating
@@ -1242,7 +1242,7 @@ contains
 
 #endif
 
-    end subroutine s_mpi_sendrecv_variables_buffers ! ---------
+    end subroutine s_mpi_sendrecv_variables_buffers
 
     !>  The goal of this procedure is to populate the buffers of
         !!      the cell-average conservative variables by communicating
@@ -2034,7 +2034,7 @@ contains
 
 #endif
 
-    end subroutine s_mpi_sendrecv_ib_buffers ! ---------
+    end subroutine s_mpi_sendrecv_ib_buffers
 
     subroutine s_mpi_sendrecv_capilary_variables_buffers(c_divs_vf, mpi_dir, pbc_loc)
 
@@ -2244,7 +2244,7 @@ contains
     end subroutine s_mpi_sendrecv_capilary_variables_buffers
 
     !> Module deallocation and/or disassociation procedures
-    subroutine s_finalize_mpi_proxy_module() ! -----------------------------
+    subroutine s_finalize_mpi_proxy_module
 
 #ifdef MFC_MPI
 
@@ -2260,6 +2260,6 @@ contains
 
 #endif
 
-    end subroutine s_finalize_mpi_proxy_module ! ---------------------------
+    end subroutine s_finalize_mpi_proxy_module
 
 end module m_mpi_proxy
diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp
index b59adfaf2..8f451c3cf 100644
--- a/src/simulation/m_rhs.fpp
+++ b/src/simulation/m_rhs.fpp
@@ -228,7 +228,7 @@ contains
     !> The computation of parameters, the allocation of memory,
         !!      the association of pointers and/or the execution of any
         !!      other procedures that are necessary to setup the module.
-    subroutine s_initialize_rhs_module() ! ---------------------------------
+    subroutine s_initialize_rhs_module
 
         integer :: i, j, k, l, id !< Generic loop iterators
 
@@ -694,9 +694,9 @@ contains
             @:ALLOCATE_GLOBAL(nbub(0:m, 0:n, 0:p))
         end if
 
-    end subroutine s_initialize_rhs_module ! -------------------------------
+    end subroutine s_initialize_rhs_module
 
-    subroutine s_compute_rhs(q_cons_vf, q_prim_vf, rhs_vf, pb, rhs_pb, mv, rhs_mv, t_step, time_avg) ! -------
+    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
@@ -1037,7 +1037,7 @@ contains
         end if
         ! ==================================================================
 
-    end subroutine s_compute_rhs ! -----------------------------------------
+    end subroutine s_compute_rhs
 
     subroutine s_compute_advection_source_term(idir, rhs_vf, q_cons_vf, q_prim_vf, flux_src_n_vf)
 
@@ -1835,7 +1835,7 @@ contains
         !!      purpose, this pressure is finally corrected using the
         !!      mixture-total-energy equation.
         !!  @param q_cons_vf Cell-average conservative variables
-    subroutine s_pressure_relaxation_procedure(q_cons_vf) ! ----------------
+    subroutine s_pressure_relaxation_procedure(q_cons_vf)
 
         type(scalar_field), dimension(sys_size), intent(INOUT) :: q_cons_vf
 
@@ -2079,7 +2079,7 @@ contains
             end do
         end do
 
-    end subroutine s_pressure_relaxation_procedure ! -----------------------
+    end subroutine s_pressure_relaxation_procedure
 
     !>  The purpose of this subroutine is to WENO-reconstruct the
         !!      left and the right cell-boundary values, including values
@@ -2091,7 +2091,7 @@ contains
         !!  @param vR_qp Right WENO-reconstructed, cell-boundary values including
         !!          the values at the quadrature points, of the cell-average variables
         !!  @param norm_dir Splitting coordinate direction
-    subroutine s_reconstruct_cell_boundary_values(v_vf, vL_x, vL_y, vL_z, vR_x, vR_y, vR_z, & ! -
+    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
@@ -2144,9 +2144,9 @@ contains
         end if
 
         ! ==================================================================
-    end subroutine s_reconstruct_cell_boundary_values ! --------------------
+    end subroutine s_reconstruct_cell_boundary_values
 
-    subroutine s_reconstruct_cell_boundary_values_first_order(v_vf, vL_x, vL_y, vL_z, vR_x, vR_y, vR_z, & ! -
+    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
@@ -2225,7 +2225,7 @@ contains
     end subroutine s_reconstruct_cell_boundary_values_first_order
 
     !> Module deallocation and/or disassociation procedures
-    subroutine s_finalize_rhs_module() ! -----------------------------------
+    subroutine s_finalize_rhs_module
 
         integer :: i, j, k, l !< Generic loop iterators
 
@@ -2381,6 +2381,6 @@ contains
         s_riemann_solver => null()
         s_convert_to_mixture_variables => null()
 
-    end subroutine s_finalize_rhs_module ! ---------------------------------
+    end subroutine s_finalize_rhs_module
 
 end module m_rhs
diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp
index 2d0cb2bf8..58489cd47 100644
--- a/src/simulation/m_riemann_solvers.fpp
+++ b/src/simulation/m_riemann_solvers.fpp
@@ -120,7 +120,7 @@ module m_riemann_solvers
         !! For more information please refer to:
         !!      1) s_compute_cartesian_viscous_source_flux
         !!      2) s_compute_cylindrical_viscous_source_flux
-        subroutine s_compute_abstract_viscous_source_flux(velL_vf, & ! -------------
+        subroutine s_compute_abstract_viscous_source_flux(velL_vf, &
                                                           dvelL_dx_vf, &
                                                           dvelL_dy_vf, &
                                                           dvelL_dz_vf, &
@@ -266,7 +266,7 @@ contains
 
     @:s_compute_speed_of_sound()
 
-    subroutine s_hll_riemann_solver(qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, dqL_prim_dx_vf, & ! -------
+    subroutine s_hll_riemann_solver(qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, dqL_prim_dx_vf, &
                                     dqL_prim_dy_vf, &
                                     dqL_prim_dz_vf, &
                                     qL_prim_vf, &
@@ -807,7 +807,7 @@ contains
         !!  @param iy Index bounds in the y-dir
         !!  @param iz Index bounds in the z-dir
         !!  @param q_prim_vf Cell-averaged primitive variables
-    subroutine s_hllc_riemann_solver(qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, dqL_prim_dx_vf, & ! ------
+    subroutine s_hllc_riemann_solver(qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, dqL_prim_dx_vf, &
                                      dqL_prim_dy_vf, &
                                      dqL_prim_dz_vf, &
                                      qL_prim_vf, &
@@ -2328,7 +2328,7 @@ contains
     !>  The computation of parameters, the allocation of memory,
         !!      the association of pointers and/or the execution of any
         !!      other procedures that are necessary to setup the module.
-    subroutine s_initialize_riemann_solvers_module() ! ---------------------
+    subroutine s_initialize_riemann_solvers_module
 
         ! Allocating the variables that will be utilized to formulate the
         ! left, right, and average states of the Riemann problem, as well
@@ -2454,7 +2454,7 @@ contains
                 is3%beg:is3%end, 1:2))
         end if
 
-    end subroutine s_initialize_riemann_solvers_module ! -------------------
+    end subroutine s_initialize_riemann_solvers_module
 
     !>  The purpose of this subroutine is to populate the buffers
         !!      of the left and right Riemann states variables, depending
@@ -2481,7 +2481,7 @@ contains
         !!  @param ix Index bounds in the x-dir
         !!  @param iy Index bounds in the y-dir
         !!  @param iz Index bounds in the z-dir
-    subroutine s_populate_riemann_states_variables_buffers( & ! ------------
+    subroutine s_populate_riemann_states_variables_buffers( &
         qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, dqL_prim_dx_vf, &
         dqL_prim_dy_vf, &
         dqL_prim_dz_vf, &
@@ -2846,7 +2846,7 @@ contains
         end if
         ! END: Population of Buffers in z-direction ========================
 
-    end subroutine s_populate_riemann_states_variables_buffers ! -----------
+    end subroutine s_populate_riemann_states_variables_buffers
 
     !>  The computation of parameters, the allocation of memory,
         !!      the association of pointers and/or the execution of any
@@ -2980,7 +2980,7 @@ contains
 
         ! ==================================================================
 
-    end subroutine s_initialize_riemann_solver ! ---------------------------
+    end subroutine s_initialize_riemann_solver
 
     !>  The goal of this subroutine is to evaluate and account
         !!      for the contribution of viscous stresses in the source
@@ -2998,7 +2998,7 @@ contains
         !!  @param ix Index bounds in  first coordinate direction
         !!  @param iy Index bounds in second coordinate direction
         !!  @param iz Index bounds in  third coordinate direction
-    subroutine s_compute_cylindrical_viscous_source_flux(velL_vf, & ! -------------
+    subroutine s_compute_cylindrical_viscous_source_flux(velL_vf, &
                                                          dvelL_dx_vf, &
                                                          dvelL_dy_vf, &
                                                          dvelL_dz_vf, &
@@ -3507,7 +3507,7 @@ contains
         end if
         ! END: Viscous Stresses in theta-direction =============================
 
-    end subroutine s_compute_cylindrical_viscous_source_flux ! -------------------------
+    end subroutine s_compute_cylindrical_viscous_source_flux
 
     !>  The goal of this subroutine is to evaluate and account
         !!      for the contribution of viscous stresses in the source
@@ -3525,7 +3525,7 @@ contains
         !!  @param ix Index bounds in  first coordinate direction
         !!  @param iy Index bounds in second coordinate direction
         !!  @param iz Index bounds in  third coordinate direction
-    subroutine s_compute_cartesian_viscous_source_flux(velL_vf, & ! -------------
+    subroutine s_compute_cartesian_viscous_source_flux(velL_vf, &
                                                        dvelL_dx_vf, &
                                                        dvelL_dy_vf, &
                                                        dvelL_dz_vf, &
@@ -4000,7 +4000,7 @@ contains
         end if
         ! END: Viscous Stresses in z-direction =============================
 
-    end subroutine s_compute_cartesian_viscous_source_flux ! -------------------------
+    end subroutine s_compute_cartesian_viscous_source_flux
 
     !>  Deallocation and/or disassociation procedures that are
         !!      needed to finalize the selected Riemann problem solver
@@ -4011,7 +4011,7 @@ contains
         !!  @param ix   Index bounds in  first coordinate direction
         !!  @param iy   Index bounds in second coordinate direction
         !!  @param iz   Index bounds in  third coordinate direction
-    subroutine s_finalize_riemann_solver(flux_vf, flux_src_vf, & ! --------
+    subroutine s_finalize_riemann_solver(flux_vf, flux_src_vf, &
                                          flux_gsrc_vf, &
                                          norm_dir, ix, iy, iz)
 
@@ -4173,10 +4173,10 @@ contains
 
         ! ==================================================================
 
-    end subroutine s_finalize_riemann_solver ! -----------------------------
+    end subroutine s_finalize_riemann_solver
 
     !> Module deallocation and/or disassociation procedures
-    subroutine s_finalize_riemann_solvers_module() ! -----------------------
+    subroutine s_finalize_riemann_solvers_module
 
         ! Disassociating procedural pointer to the subroutine which was
         ! utilized to calculate the solution of a given Riemann problem
@@ -4227,6 +4227,6 @@ contains
             @:DEALLOCATE_GLOBAL(mom_sp_rsz_vf)
         end if
 
-    end subroutine s_finalize_riemann_solvers_module ! ---------------------
+    end subroutine s_finalize_riemann_solvers_module
 
 end module m_riemann_solvers
diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp
index 582b7231a..e8911ad29 100644
--- a/src/simulation/m_start_up.fpp
+++ b/src/simulation/m_start_up.fpp
@@ -91,7 +91,7 @@ module m_start_up
     abstract interface ! ===================================================
 
         !! @param q_cons_vf  Conservative variables
-        subroutine s_read_abstract_data_files(q_cons_vf) ! -----------
+        subroutine s_read_abstract_data_files(q_cons_vf)
 
             import :: scalar_field, sys_size, pres_field
 
@@ -99,7 +99,7 @@ module m_start_up
                 dimension(sys_size), &
                 intent(INOUT) :: q_cons_vf
 
-        end subroutine s_read_abstract_data_files ! -----------------
+        end subroutine s_read_abstract_data_files
 
     end interface ! ========================================================
 
@@ -112,7 +112,7 @@ contains
     !>  The purpose of this procedure is to first verify that an
         !!      input file has been made available by the user. Provided
         !!      that this is so, the input file is then read in.
-    subroutine s_read_input_file() ! ---------------------------------------
+    subroutine s_read_input_file
 
         ! Relative path to the input file provided by the user
         character(LEN=name_len) :: file_path = './simulation.inp'
@@ -204,12 +204,12 @@ contains
 #endif
 #endif
 
-    end subroutine s_read_input_file ! -------------------------------------
+    end subroutine s_read_input_file
 
     !> The goal of this procedure is to verify that each of the
     !!      user provided inputs is valid and that their combination
     !!      constitutes a meaningful configuration for the simulation.
-    subroutine s_check_input_file() ! --------------------------------------
+    subroutine s_check_input_file
 
         ! Relative path to the current directory file in the case directory
         character(LEN=path_len) :: file_path
@@ -232,7 +232,7 @@ contains
 
         call s_check_inputs()
 
-    end subroutine s_check_input_file ! ------------------------------------
+    end subroutine s_check_input_file
 
         !!              initial condition and grid data files. The cell-average
         !!              conservative variables constitute the former, while the
@@ -240,7 +240,7 @@ contains
         !!              up the latter. This procedure also calculates the cell-
         !!              width distributions from the cell-boundary locations.
         !! @param q_cons_vf Cell-averaged conservative variables
-    subroutine s_read_serial_data_files(q_cons_vf) ! ------------------------------
+    subroutine s_read_serial_data_files(q_cons_vf)
 
         type(scalar_field), dimension(sys_size), intent(INOUT) :: q_cons_vf
 
@@ -447,10 +447,10 @@ contains
 
         end if
 
-    end subroutine s_read_serial_data_files ! -------------------------------------
+    end subroutine s_read_serial_data_files
 
         !! @param q_cons_vf Conservative variables
-    subroutine s_read_parallel_data_files(q_cons_vf) ! ---------------------------
+    subroutine s_read_parallel_data_files(q_cons_vf)
 
         type(scalar_field), &
             dimension(sys_size), &
@@ -809,13 +809,13 @@ contains
 
 #endif
 
-    end subroutine s_read_parallel_data_files ! -------------------------------
+    end subroutine s_read_parallel_data_files
 
     !> The purpose of this subroutine is to populate the buffers
         !!          of the grid variables, which are constituted of the cell-
         !!          boundary locations and cell-width distributions, based on
         !!          the boundary conditions.
-    subroutine s_populate_grid_variables_buffers() ! -----------------------
+    subroutine s_populate_grid_variables_buffers
 
         integer :: i !< Generic loop iterator
 
@@ -1021,14 +1021,14 @@ contains
 
         ! END: Population of Buffers in z-direction ========================
 
-    end subroutine s_populate_grid_variables_buffers ! ---------------------
+    end subroutine s_populate_grid_variables_buffers
 
     !> The purpose of this procedure is to initialize the
         !!      values of the internal-energy equations of each phase
         !!      from the mass of each phase, the mixture momentum and
         !!      mixture-total-energy equations.
         !! @param v_vf conservative variables
-    subroutine s_initialize_internal_energy_equations(v_vf) !---------------
+    subroutine s_initialize_internal_energy_equations(v_vf)
 
         type(scalar_field), dimension(sys_size), intent(INOUT) :: v_vf
         real(kind(0d0)) :: rho
diff --git a/src/simulation/m_surface_tension.fpp b/src/simulation/m_surface_tension.fpp
index ebe0cd70a..de1d56794 100644
--- a/src/simulation/m_surface_tension.fpp
+++ b/src/simulation/m_surface_tension.fpp
@@ -318,7 +318,7 @@ contains
 
     end subroutine s_get_capilary
 
-    subroutine s_reconstruct_cell_boundary_values_capillary(v_vf, vL_x, vL_y, vL_z, vR_x, vR_y, vR_z, & ! -
+    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
diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp
index 470d91477..f600093c1 100644
--- a/src/simulation/m_time_steppers.fpp
+++ b/src/simulation/m_time_steppers.fpp
@@ -89,7 +89,7 @@ contains
     !> The computation of parameters, the allocation of memory,
         !!      the association of pointers and/or the execution of any
         !!      other procedures that are necessary to setup the module.
-    subroutine s_initialize_time_steppers_module() ! -----------------------
+    subroutine s_initialize_time_steppers_module
 
         type(int_bounds_info) :: ix_t, iy_t, iz_t !<
             !! Indical bounds in the x-, y- and z-directions
@@ -285,11 +285,11 @@ contains
             call s_open_run_time_information_file()
         end if
 
-    end subroutine s_initialize_time_steppers_module ! ---------------------
+    end subroutine s_initialize_time_steppers_module
 
     !> 1st order TVD RK time-stepping algorithm
         !! @param t_step Current time step
-    subroutine s_1st_order_tvd_rk(t_step, time_avg) ! --------------------------------
+    subroutine s_1st_order_tvd_rk(t_step, time_avg)
 
         integer, intent(IN) :: t_step
         real(kind(0d0)), intent(INOUT) :: time_avg
@@ -399,11 +399,11 @@ contains
 
         ! ==================================================================
 
-    end subroutine s_1st_order_tvd_rk ! ------------------------------------
+    end subroutine s_1st_order_tvd_rk
 
     !> 2nd order TVD RK time-stepping algorithm
         !! @param t_step Current time-step
-    subroutine s_2nd_order_tvd_rk(t_step, time_avg) ! --------------------------------
+    subroutine s_2nd_order_tvd_rk(t_step, time_avg)
 
         integer, intent(IN) :: t_step
         real(kind(0d0)), intent(INOUT) :: time_avg
@@ -586,11 +586,11 @@ contains
         call cpu_time(finish)
         ! ==================================================================
 
-    end subroutine s_2nd_order_tvd_rk ! ------------------------------------
+    end subroutine s_2nd_order_tvd_rk
 
     !> 3rd order TVD RK time-stepping algorithm
         !! @param t_step Current time-step
-    subroutine s_3rd_order_tvd_rk(t_step, time_avg, dt_in) ! --------------------------------
+    subroutine s_3rd_order_tvd_rk(t_step, time_avg, dt_in)
 
         integer, intent(IN) :: t_step
         real(kind(0d0)), intent(INOUT) :: time_avg
@@ -855,13 +855,13 @@ contains
         end if
         ! ==================================================================
 
-    end subroutine s_3rd_order_tvd_rk ! ------------------------------------
+    end subroutine s_3rd_order_tvd_rk
 
     !> Strang splitting scheme with 3rd order TVD RK time-stepping algorithm for
         !!      the flux term and adaptive time stepping algorithm for
         !!      the source term
         !! @param t_step Current time-step
-    subroutine s_strang_splitting(t_step, time_avg) ! --------------------------------
+    subroutine s_strang_splitting(t_step, time_avg)
 
         integer, intent(IN) :: t_step
         real(kind(0d0)), intent(INOUT) :: time_avg
@@ -890,11 +890,11 @@ contains
 
         ! ==================================================================
 
-    end subroutine s_strang_splitting ! ------------------------------------
+    end subroutine s_strang_splitting
 
     !> Bubble source part in Strang operator splitting scheme
         !! @param q_cons_vf conservative variables
-    subroutine s_adaptive_dt_bubble(t_step) ! ------------------------
+    subroutine s_adaptive_dt_bubble(t_step)
 
         integer, intent(IN) :: t_step
 
@@ -913,7 +913,7 @@ contains
 
         call s_compute_bubble_source(q_cons_ts(1)%vf, q_prim_vf, t_step, rhs_vf)
 
-    end subroutine s_adaptive_dt_bubble ! ------------------------------
+    end subroutine s_adaptive_dt_bubble
 
     !> This subroutine applies the body forces source term at each
         !! Runge-Kutta stage
@@ -943,7 +943,7 @@ contains
     !> This subroutine saves the temporary q_prim_vf vector
         !!      into the q_prim_ts vector that is then used in p_main
         !! @param t_step current time-step
-    subroutine s_time_step_cycling(t_step) ! ----------------------------
+    subroutine s_time_step_cycling(t_step)
 
         integer, intent(IN) :: t_step
 
@@ -978,10 +978,9 @@ contains
             end do
         end if
 
-    end subroutine s_time_step_cycling ! -----------------------------------
-
+    end subroutine s_time_step_cycling
     !> Module deallocation and/or disassociation procedures
-    subroutine s_finalize_time_steppers_module() ! -------------------------
+    subroutine s_finalize_time_steppers_module
 
         integer :: i, j !< Generic loop iterators
 
@@ -1046,6 +1045,6 @@ contains
             call s_close_run_time_information_file()
         end if
 
-    end subroutine s_finalize_time_steppers_module ! -----------------------
+    end subroutine s_finalize_time_steppers_module
 
 end module m_time_steppers
diff --git a/src/simulation/m_viscous.fpp b/src/simulation/m_viscous.fpp
index a3ab33325..a6f1892b5 100644
--- a/src/simulation/m_viscous.fpp
+++ b/src/simulation/m_viscous.fpp
@@ -72,7 +72,7 @@ contains
     !  @param grad_z_vf Cell-average primitive variable derivatives, z-dir
     subroutine s_compute_viscous_stress_tensor(q_prim_vf, grad_x_vf, grad_y_vf, grad_z_vf, &
                                                tau_Re_vf, &
-                                               ix, iy, iz) ! ---
+                                               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
@@ -522,7 +522,7 @@ contains
                 end do
             end do
         end if
-    end subroutine s_compute_viscous_stress_tensor ! ----------------------------------------
+    end subroutine s_compute_viscous_stress_tensor
 
 !>  Computes viscous terms
     !!  @param q_cons_vf Cell-averaged conservative variables
@@ -981,7 +981,7 @@ contains
 
     end subroutine s_get_viscous
 
-    subroutine s_reconstruct_cell_boundary_values_visc(v_vf, vL_x, vL_y, vL_z, vR_x, vR_y, vR_z, & ! -
+    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
@@ -1082,9 +1082,9 @@ contains
 
         ! ==================================================================
 
-    end subroutine s_reconstruct_cell_boundary_values_visc ! --------------------
+    end subroutine s_reconstruct_cell_boundary_values_visc
 
-    subroutine s_reconstruct_cell_boundary_values_visc_deriv(v_vf, vL_x, vL_y, vL_z, vR_x, vR_y, vR_z, & ! -
+    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
@@ -1184,7 +1184,7 @@ contains
         end if
         ! ==================================================================
 
-    end subroutine s_reconstruct_cell_boundary_values_visc_deriv ! --------------------
+    end subroutine s_reconstruct_cell_boundary_values_visc_deriv
 
     !>  The purpose of this subroutine is to employ the inputted
         !!      left and right cell-boundary integral-averaged variables
@@ -1195,7 +1195,7 @@ contains
         !!  @param vR_vf Right cell-boundary integral averages
         !!  @param dv_ds_vf Cell-average first-order spatial derivatives
         !!  @param norm_dir Splitting coordinate direction
-    subroutine s_apply_scalar_divergence_theorem(vL_vf, vR_vf, & ! --------
+    subroutine s_apply_scalar_divergence_theorem(vL_vf, vR_vf, &
                                                  dv_ds_vf, &
                                                  norm_dir, &
                                                  ix, iy, iz, iv_in, &
@@ -1311,7 +1311,7 @@ contains
         end if
         ! END: First-Order Spatial Derivatives in z-direction ==============
 
-    end subroutine s_apply_scalar_divergence_theorem ! ---------------------
+    end subroutine s_apply_scalar_divergence_theorem
 
     !>  Computes the scalar gradient fields via finite differences
         !!  @param var Variable to compute derivative of
@@ -1501,9 +1501,9 @@ contains
             end if
         end if
 
-    end subroutine s_compute_fd_gradient ! --------------------------------------
+    end subroutine s_compute_fd_gradient
 
-    subroutine s_finalize_viscous_module()
+    subroutine s_finalize_viscous_module
 
         integer :: i
 
diff --git a/src/simulation/m_weno.fpp b/src/simulation/m_weno.fpp
index 9e2d38f2b..8201685b8 100644
--- a/src/simulation/m_weno.fpp
+++ b/src/simulation/m_weno.fpp
@@ -158,7 +158,7 @@ contains
     !>  The computation of parameters, the allocation of memory,
         !!      the association of pointers and/or the execution of any
         !!      other procedures that are necessary to setup the module.
-    subroutine s_initialize_weno_module() ! --------------------------------
+    subroutine s_initialize_weno_module
 
         integer :: i, j
         if (weno_order == 1) return
@@ -256,7 +256,7 @@ contains
 
         ! ==================================================================
 
-    end subroutine s_initialize_weno_module ! ------------------------------
+    end subroutine s_initialize_weno_module
 
     !>  The purpose of this subroutine is to compute the grid
         !!      dependent coefficients of the WENO polynomials, ideal
@@ -265,7 +265,7 @@ contains
         !!      reconstruction.
         !! @param weno_dir Coordinate direction of the WENO reconstruction
         !! @param is Index bounds in the s-direction
-    subroutine s_compute_weno_coefficients(weno_dir, is) ! -------
+    subroutine s_compute_weno_coefficients(weno_dir, is)
 
         integer, intent(IN) :: weno_dir
         type(int_bounds_info), intent(IN) :: is
@@ -500,9 +500,9 @@ contains
 
         nullify (s_cb)
 
-    end subroutine s_compute_weno_coefficients ! ---------------------------
+    end subroutine s_compute_weno_coefficients
 
-    subroutine s_weno(v_vf, vL_rs_vf_x, vL_rs_vf_y, vL_rs_vf_z, vR_rs_vf_x, vR_rs_vf_y, vR_rs_vf_z, & ! -------------------
+    subroutine s_weno(v_vf, vL_rs_vf_x, vL_rs_vf_y, vL_rs_vf_z, vR_rs_vf_x, vR_rs_vf_y, vR_rs_vf_z, &
                       norm_dir, weno_dir, &
                       is1_weno_d, is2_weno_d, is3_weno_d)
 
@@ -789,7 +789,7 @@ contains
         !! @param is1_weno Index bounds in first coordinate direction
         !! @param is2_weno Index bounds in second coordinate direction
         !! @param is3_weno Index bounds in third coordinate direction
-    subroutine s_initialize_weno(v_vf, & ! ---------
+    subroutine s_initialize_weno(v_vf, &
                                  norm_dir, weno_dir)
 
         type(scalar_field), dimension(:), intent(IN) :: v_vf
@@ -897,7 +897,7 @@ contains
 
         ! ==================================================================
 
-    end subroutine s_initialize_weno ! -------------------------------------
+    end subroutine s_initialize_weno
 
     !>  The goal of this subroutine is to ensure that the WENO
         !!      reconstruction is monotonic. The latter is achieved by
@@ -910,7 +910,7 @@ contains
         !!  @param j First-coordinate cell index
         !!  @param k Second-coordinate cell index
         !!  @param l Third-coordinate cell index
-    subroutine s_preserve_monotonicity(v_rs_ws, vL_rs_vf, vR_rs_vf) ! --------------------------
+    subroutine s_preserve_monotonicity(v_rs_ws, vL_rs_vf, vR_rs_vf)
 
         real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(IN) :: v_rs_ws
         real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(INOUT) :: vL_rs_vf, vR_rs_vf
@@ -1071,10 +1071,10 @@ contains
         end do
         !$acc end parallel loop
 
-    end subroutine s_preserve_monotonicity ! -------------------------------
+    end subroutine s_preserve_monotonicity
 
     !>  Module deallocation and/or disassociation procedures
-    subroutine s_finalize_weno_module() ! ----------------------------------
+    subroutine s_finalize_weno_module()
 
         integer :: i, j
 
@@ -1113,6 +1113,6 @@ contains
         @:DEALLOCATE_GLOBAL(beta_coef_z)
         ! ==================================================================
 
-    end subroutine s_finalize_weno_module ! --------------------------------
+    end subroutine s_finalize_weno_module
 
 end module m_weno