m_data_output.fpp.f90

Go to the documentation of this file.

# 1 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
!>
!! @file
!! @brief Contains module m_data_output
 
!> @brief Writes post-processed grid and flow-variable data to Silo-HDF5 or binary database files
module m_data_output
 
    use m_derived_types         ! Definitions of the derived types
 
    use m_global_parameters     ! Global parameters
 
    use m_derived_variables     !< procedures used to compute quantities derived
 
    use m_mpi_proxy             ! Message passing interface (MPI) module proxy
 
    use m_compile_specific
 
    use m_helper
 
    use m_variables_conversion
 
    implicit none
 
    private; public :: s_initialize_data_output_module, &
 s_define_output_region, &
 s_open_formatted_database_file, &
 s_open_intf_data_file, &
 s_open_energy_data_file, &
 s_write_grid_to_formatted_database_file, &
 s_write_variable_to_formatted_database_file, &
 s_write_lag_bubbles_results_to_text, &
 s_write_lag_bubbles_to_formatted_database_file, &
 s_write_intf_data_file, &
 s_write_energy_data_file, &
 s_close_formatted_database_file, &
 s_close_intf_data_file, &
 s_close_energy_data_file, &
 s_finalize_data_output_module
 
    ! Including the Silo Fortran interface library that features the subroutines
    ! and parameters that are required to write in the Silo-HDF5 database format
    ! INCLUDE 'silo.inc'
    include 'silo_f9x.inc'
 
    ! Generic storage for flow variable(s) that are to be written to formatted
    ! database file(s). Note that for 1D simulations, q_root_sf is employed to
    ! gather the flow variable(s) from all sub-domains on to the root process.
    ! If the run is not parallel, but serial, then q_root_sf is equal to q_sf.
    real(wp), allocatable, dimension(:, :, :), public :: q_sf
    real(wp), allocatable, dimension(:, :, :) :: q_root_sf
    real(wp), allocatable, dimension(:, :, :) :: cyl_q_sf
 
    ! Single precision storage for flow variables
    real(sp), allocatable, dimension(:, :, :), public :: q_sf_s
    real(sp), allocatable, dimension(:, :, :) :: q_root_sf_s
    real(sp), allocatable, dimension(:, :, :) :: cyl_q_sf_s
 
    ! The spatial and data extents array variables contain information about the
    ! minimum and maximum values of the grid and flow variable(s), respectively.
    ! The purpose of bookkeeping this information is to boost the visualization
    ! of the Silo-HDF5 database file(s) in VisIt.
    real(wp), allocatable, dimension(:, :) :: spatial_extents
    real(wp), allocatable, dimension(:, :) :: data_extents
 
    ! The size of the ghost zone layer at beginning of each coordinate direction
    ! (lo) and at end of each coordinate direction (hi). Adding this information
    ! to Silo-HDF5 database file(s) is recommended since it supplies VisIt with
    ! connectivity information between the sub-domains of a parallel data set.
    integer, allocatable, dimension(:) :: lo_offset
    integer, allocatable, dimension(:) :: hi_offset
 
    ! For Silo-HDF5 database format, this variable is used to keep track of the
    ! number of cell-boundaries, for the grid associated with the local process,
    ! in each of the active coordinate directions.
    integer, allocatable, dimension(:) :: dims
 
    ! Locations of various folders in the case's directory tree, associated with
    ! the choice of the formatted database format. These include, in order, the
    ! location of the folder named after the selected formatted database format,
    ! and the locations of two sub-directories of the latter, the first of which
    ! is named after the local processor rank, while the second is named 'root'.
    ! The folder associated with the local processor rank contains only the data
    ! pertaining to the part of the domain taken care of by the local processor.
    ! The root directory, on the other hand, will contain either the information
    ! about the connectivity required to put the entire domain back together, or
    ! the actual data associated with the entire computational domain. This all
    ! depends on dimensionality and the choice of the formatted database format.
    character(LEN=path_len + name_len) :: dbdir
    character(LEN=path_len + 2*name_len) :: proc_rank_dir
    character(LEN=path_len + 2*name_len) :: rootdir
 
    ! Handles of the formatted database master/root file, slave/local processor
    ! file and options list. The list of options is explicitly used in the Silo-
    ! HDF5 database format to provide additional details about the contents of a
    ! formatted database file, such as the previously described spatial and data
    ! extents.
    integer :: dbroot
    integer :: dbfile
    integer :: optlist
 
    ! The total number of flow variable(s) to be stored in a formatted database
    ! file. Note that this is only needed when using the Binary format.
    integer :: dbvars
 
    ! Generic error flags utilized in the handling, checking and the reporting
    ! of the input and output operations errors with a formatted database file
    integer, private :: err
 
contains
 
    !> @brief Allocate storage arrays, configure output directories, and count flow variables for formatted database output.
    impure subroutine s_initialize_data_output_module()
        ! Description: Computation of parameters, allocation procedures, and/or
        !              any other tasks needed to properly setup the module
 
        ! Generic string used to store the location of a particular file
        character(LEN=len_trim(case_dir) + 2*name_len) :: file_loc
 
        ! Generic logical used to test the existence of a particular folder
        logical :: dir_check
 
        integer :: i
 
        ! Allocating the generic storage for the flow variable(s) that are
        ! going to be written to the formatted database file(s). Note once
        ! more that the root variable is only required for 1D computations.
        allocate (q_sf(-offset_x%beg:m + offset_x%end, &
                       -offset_y%beg:n + offset_y%end, &
                       -offset_z%beg:p + offset_z%end))
        if (grid_geometry == 3) then
            allocate (cyl_q_sf(-offset_y%beg:n + offset_y%end, &
                               -offset_z%beg:p + offset_z%end, &
                               -offset_x%beg:m + offset_x%end))
        end if
 
        if (precision == 1) then
            allocate (q_sf_s(-offset_x%beg:m + offset_x%end, &
                             -offset_y%beg:n + offset_y%end, &
                             -offset_z%beg:p + offset_z%end))
            if (grid_geometry == 3) then
                allocate (cyl_q_sf_s(-offset_y%beg:n + offset_y%end, &
                                     -offset_z%beg:p + offset_z%end, &
                                     -offset_x%beg:m + offset_x%end))
            end if
        end if
 
        if (n == 0) then
            allocate (q_root_sf(0:m_root, 0:0, 0:0))
            if (precision == 1) then
                allocate (q_root_sf_s(0:m_root, 0:0, 0:0))
            end if
        end if
 
        ! Allocating the spatial and data extents and also the variables for
        ! the offsets and the one bookkeeping the number of cell-boundaries
        ! in each active coordinate direction. Note that all these variables
        ! are only needed by the Silo-HDF5 format for multidimensional data.
        if (format == 1) then
 
            allocate (data_extents(1:2, 0:num_procs - 1))
 
            if (p > 0) then
                allocate (spatial_extents(1:6, 0:num_procs - 1))
                allocate (lo_offset(1:3))
                allocate (hi_offset(1:3))
                allocate (dims(1:3))
            elseif (n > 0) then
                allocate (spatial_extents(1:4, 0:num_procs - 1))
                allocate (lo_offset(1:2))
                allocate (hi_offset(1:2))
                allocate (dims(1:2))
            else
                allocate (spatial_extents(1:2, 0:num_procs - 1))
                allocate (lo_offset(1:1))
                allocate (hi_offset(1:1))
                allocate (dims(1:1))
            end if
 
        end if
 
        ! The size of the ghost zone layer in each of the active coordinate
        ! directions was set in the module m_mpi_proxy.f90. The results are
        ! now transferred to the local variables of this module when they are
        ! required by the Silo-HDF5 format, for multidimensional data sets.
        ! With the same, latter, requirements, the variables bookkeeping the
        ! number of cell-boundaries in each active coordinate direction are
        ! also set here.
        if (format == 1) then
            if (p > 0) then
                if (grid_geometry == 3) then
                    lo_offset(:) = (/offset_y%beg, offset_z%beg, offset_x%beg/)
                    hi_offset(:) = (/offset_y%end, offset_z%end, offset_x%end/)
                else
                    lo_offset(:) = (/offset_x%beg, offset_y%beg, offset_z%beg/)
                    hi_offset(:) = (/offset_x%end, offset_y%end, offset_z%end/)
                end if
 
                if (grid_geometry == 3) then
                    dims(:) = (/n + offset_y%beg + offset_y%end + 2, &
                                p + offset_z%beg + offset_z%end + 2, &
                                m + offset_x%beg + offset_x%end + 2/)
                else
                    dims(:) = (/m + offset_x%beg + offset_x%end + 2, &
                                n + offset_y%beg + offset_y%end + 2, &
                                p + offset_z%beg + offset_z%end + 2/)
                end if
            elseif (n > 0) then
                lo_offset(:) = (/offset_x%beg, offset_y%beg/)
                hi_offset(:) = (/offset_x%end, offset_y%end/)
 
                dims(:) = (/m + offset_x%beg + offset_x%end + 2, &
                            n + offset_y%beg + offset_y%end + 2/)
            else
                lo_offset(:) = (/offset_x%beg/)
                hi_offset(:) = (/offset_x%end/)
                dims(:) = (/m + offset_x%beg + offset_x%end + 2/)
            end if
        end if
 
        ! Generating Silo-HDF5 Directory Tree
        if (format == 1) then
 
            ! Creating the directory associated with the local process
            dbdir = trim(case_dir)//'/silo_hdf5'
 
            write (proc_rank_dir, '(A,I0)') '/p', proc_rank
 
            proc_rank_dir = trim(dbdir)//trim(proc_rank_dir)
 
            file_loc = trim(proc_rank_dir)//'/.'
 
            call my_inquire(file_loc, dir_check)
            if (dir_check .neqv. .true.) then
                call s_create_directory(trim(proc_rank_dir))
            end if
 
            ! Creating the directory associated with the root process
            if (proc_rank == 0) then
 
                rootdir = trim(dbdir)//'/root'
 
                file_loc = trim(rootdir)//'/.'
 
                call my_inquire(file_loc, dir_check)
                if (dir_check .neqv. .true.) then
                    call s_create_directory(trim(rootdir))
                end if
 
            end if
 
            ! Generating Binary Directory Tree
 
        else
 
            ! Creating the directory associated with the local process
            dbdir = trim(case_dir)//'/binary'
 
            write (proc_rank_dir, '(A,I0)') '/p', proc_rank
 
            proc_rank_dir = trim(dbdir)//trim(proc_rank_dir)
 
            file_loc = trim(proc_rank_dir)//'/.'
 
            call my_inquire(file_loc, dir_check)
 
            if (dir_check .neqv. .true.) then
                call s_create_directory(trim(proc_rank_dir))
            end if
 
            ! Creating the directory associated with the root process
            if (n == 0 .and. proc_rank == 0) then
 
                rootdir = trim(dbdir)//'/root'
 
                file_loc = trim(rootdir)//'/.'
 
                call my_inquire(file_loc, dir_check)
 
                if (dir_check .neqv. .true.) then
                    call s_create_directory(trim(rootdir))
                end if
 
            end if
 
        end if
 
        if (bubbles_lagrange) then !Lagrangian solver
            if (lag_txt_wrt) then
                dbdir = trim(case_dir)//'/lag_bubbles_post_process'
                file_loc = trim(dbdir)//'/.'
                call my_inquire(file_loc, dir_check)
 
                if (dir_check .neqv. .true.) then
                    call s_create_directory(trim(dbdir))
                end if
            end if
        end if
 
        ! Contrary to the Silo-HDF5 database format, handles of the Binary
        ! database master/root and slave/local process files are perfectly
        ! static throughout post-process. Hence, they are set here so that
        ! they do not have to be repetitively computed in later procedures.
        if (format == 2) then
            if (n == 0 .and. proc_rank == 0) dbroot = 2
            dbfile = 1
        end if
 
        ! Querying Number of Flow Variable(s) in Binary Output
 
        if (format == 2) then
 
            ! Initializing the counter of the number of flow variable(s) to
            ! be written to the formatted database file(s)
            dbvars = 0
 
            ! Partial densities
            if ((model_eqns == 2) .or. (model_eqns == 3)) then
                do i = 1, num_fluids
                    if (alpha_rho_wrt(i) &
                        .or. &
                        (cons_vars_wrt .or. prim_vars_wrt)) then
                        dbvars = dbvars + 1
                    end if
                end do
            end if
 
            ! Density
            if ((rho_wrt .or. (model_eqns == 1 .and. (cons_vars_wrt .or. prim_vars_wrt))) &
                .and. (.not. relativity)) then
                dbvars = dbvars + 1
            end if
 
            if (relativity .and. (rho_wrt .or. prim_vars_wrt)) dbvars = dbvars + 1
            if (relativity .and. (rho_wrt .or. cons_vars_wrt)) dbvars = dbvars + 1
 
            ! Momentum
            do i = 1, e_idx - mom_idx%beg
                if (mom_wrt(i) .or. cons_vars_wrt) dbvars = dbvars + 1
            end do
 
            ! Velocity
            do i = 1, e_idx - mom_idx%beg
                if (vel_wrt(i) .or. prim_vars_wrt) dbvars = dbvars + 1
            end do
 
            ! Flux limiter function
            do i = 1, e_idx - mom_idx%beg
                if (flux_wrt(i)) dbvars = dbvars + 1
            end do
 
            ! Energy
            if (e_wrt .or. cons_vars_wrt) dbvars = dbvars + 1
 
            ! Pressure
            if (pres_wrt .or. prim_vars_wrt) dbvars = dbvars + 1
 
            ! Elastic stresses
            if (hypoelasticity) dbvars = dbvars + (num_dims*(num_dims + 1))/2
 
            ! Damage state variable
            if (cont_damage) dbvars = dbvars + 1
 
            ! Hyperbolic cleaning for MHD
            if (hyper_cleaning) dbvars = dbvars + 1
 
            ! Magnetic field
            if (mhd) then
                if (n == 0) then
                    dbvars = dbvars + 2
                else
                    dbvars = dbvars + 3
                end if
            end if
 
            ! Volume fraction(s)
            if ((model_eqns == 2) .or. (model_eqns == 3)) then
 
                do i = 1, num_fluids - 1
                    if (alpha_wrt(i) &
                        .or. &
                        (cons_vars_wrt .or. prim_vars_wrt)) then
                        dbvars = dbvars + 1
                    end if
                end do
 
                if (alpha_wrt(num_fluids) &
                    .or. &
                    (cons_vars_wrt .or. prim_vars_wrt)) &
                    then
                    dbvars = dbvars + 1
                end if
 
            end if
 
            ! Specific heat ratio function
            if (gamma_wrt &
                .or. &
                (model_eqns == 1 .and. (cons_vars_wrt .or. prim_vars_wrt))) &
                then
                dbvars = dbvars + 1
            end if
 
            ! Specific heat ratio
            if (heat_ratio_wrt) dbvars = dbvars + 1
 
            ! Liquid stiffness function
            if (pi_inf_wrt &
                .or. &
                (model_eqns == 1 .and. (cons_vars_wrt .or. prim_vars_wrt))) &
                then
                dbvars = dbvars + 1
            end if
 
            ! Liquid stiffness
            if (pres_inf_wrt) dbvars = dbvars + 1
 
            ! Speed of sound
            if (c_wrt) dbvars = dbvars + 1
 
            ! Vorticity
            if (p > 0) then
                do i = 1, num_vels
                    if (omega_wrt(i)) dbvars = dbvars + 1
                end do
            elseif (n > 0) then
                do i = 1, num_vels
                    if (omega_wrt(i)) dbvars = dbvars + 1
                end do
            end if
 
            ! Numerical Schlieren function
            if (schlieren_wrt) dbvars = dbvars + 1
 
        end if
 
        ! END: Querying Number of Flow Variable(s) in Binary Output
 
    end subroutine s_initialize_data_output_module
 
    !> @brief Compute the cell-index bounds for the user-specified partial output domain in each coordinate direction.
    impure subroutine s_define_output_region
 
        integer :: i
        integer :: lower_bound, upper_bound
 
# 447 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
 
            if (m == 0) return ! Early return for y or z if simulation is 1D or 2D
 
            lower_bound = -offset_x%beg
            upper_bound = m+offset_x%end
 
            do i = lower_bound, upper_bound
                if (x_cc(i) > x_output%beg) then
                    x_output_idx%beg = i + offset_x%beg
                    exit
                end if
            end do
 
            do i = upper_bound, lower_bound, -1
                if (x_cc(i) < x_output%end) then
                    x_output_idx%end = i + offset_x%beg
                    exit
                end if
            end do
 
            ! If no grid points are within the output region
            if ((x_cc(lower_bound) > x_output%end) .or. (x_cc(upper_bound) < x_output%beg)) then
                x_output_idx%beg = 0
                x_output_idx%end = 0
            end if
 
# 447 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
 
            if (n == 0) return ! Early return for y or z if simulation is 1D or 2D
 
            lower_bound = -offset_y%beg
            upper_bound = n+offset_y%end
 
            do i = lower_bound, upper_bound
                if (y_cc(i) > y_output%beg) then
                    y_output_idx%beg = i + offset_y%beg
                    exit
                end if
            end do
 
            do i = upper_bound, lower_bound, -1
                if (y_cc(i) < y_output%end) then
                    y_output_idx%end = i + offset_y%beg
                    exit
                end if
            end do
 
            ! If no grid points are within the output region
            if ((y_cc(lower_bound) > y_output%end) .or. (y_cc(upper_bound) < y_output%beg)) then
                y_output_idx%beg = 0
                y_output_idx%end = 0
            end if
 
# 447 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
 
            if (p == 0) return ! Early return for y or z if simulation is 1D or 2D
 
            lower_bound = -offset_z%beg
            upper_bound = p+offset_z%end
 
            do i = lower_bound, upper_bound
                if (z_cc(i) > z_output%beg) then
                    z_output_idx%beg = i + offset_z%beg
                    exit
                end if
            end do
 
            do i = upper_bound, lower_bound, -1
                if (z_cc(i) < z_output%end) then
                    z_output_idx%end = i + offset_z%beg
                    exit
                end if
            end do
 
            ! If no grid points are within the output region
            if ((z_cc(lower_bound) > z_output%end) .or. (z_cc(upper_bound) < z_output%beg)) then
                z_output_idx%beg = 0
                z_output_idx%end = 0
            end if
 
# 474 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
 
    end subroutine s_define_output_region
 
    !> @brief Open (or create) the Silo-HDF5 or Binary formatted database slave and master files for a given time step.
    impure 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
        !              current time-step, t_step. This is performed by all the
        !              local process(es). The root processor, in addition, must
        !              also generate a master formatted database file whose job
        !              will be to link, and thus combine, the data from all of
        !              the local process(es). Note that for the Binary format,
        !              this extra task that is assigned to the root process is
        !              not performed in multidimensions.
 
        ! Time-step that is currently being post-processed
        integer, intent(IN) :: t_step
 
        ! Generic string used to store the location of a particular file
        character(LEN=len_trim(case_dir) + 3*name_len) :: file_loc
 
        integer :: ierr !< Generic flag used to identify and report database errors
 
        ! Silo-HDF5 Database Format
 
        if (format == 1) then
 
            ! Generating the relative path to the formatted database slave
            ! file, that is to be opened for the current time-step, t_step
            write (file_loc, '(A,I0,A)') '/', t_step, '.silo'
            file_loc = trim(proc_rank_dir)//trim(file_loc)
 
            ! Creating formatted database slave file at the above location
            ! and setting up the structure of the file and its header info
            ierr = dbcreate(trim(file_loc), len_trim(file_loc), &
                            db_clobber, db_local, 'MFC v3.0', 8, &
                            db_hdf5, dbfile)
 
            ! Verifying that the creation and setup process of the formatted
            ! database slave file has been performed without errors. If this
            ! is not the case, the post-process exits.
            if (dbfile == -1) then
                call s_mpi_abort('Unable to create Silo-HDF5 database '// &
                                 'slave file '//trim(file_loc)//'. '// &
                                 'Exiting.')
            end if
 
            ! Next, analogous steps to the ones above are carried out by the
            ! root process to create and setup the formatted database master
            ! file.
            if (proc_rank == 0) then
 
                write (file_loc, '(A,I0,A)') '/collection_', t_step, '.silo'
                file_loc = trim(rootdir)//trim(file_loc)
 
                ierr = dbcreate(trim(file_loc), len_trim(file_loc), &
                                db_clobber, db_local, 'MFC v3.0', 8, &
                                db_hdf5, dbroot)
 
                if (dbroot == -1) then
                    call s_mpi_abort('Unable to create Silo-HDF5 database '// &
                                     'master file '//trim(file_loc)//'. '// &
                                     'Exiting.')
                end if
 
            end if
 
            ! Binary Database Format
 
        else
 
            ! Generating the relative path to the formatted database slave
            ! file, that is to be opened for the current time-step, t_step
            write (file_loc, '(A,I0,A)') '/', t_step, '.dat'
            file_loc = trim(proc_rank_dir)//trim(file_loc)
 
            ! Creating the formatted database slave file, at the previously
            ! precised relative path location, and setting up its structure
            open (dbfile, iostat=err, file=trim(file_loc), &
                  form='unformatted', status='replace')
 
            ! Verifying that the creation and setup process of the formatted
            ! database slave file has been performed without errors. If this
            ! is not the case, the post-process exits.
            if (err /= 0) then
                call s_mpi_abort('Unable to create Binary database slave '// &
                                 'file '//trim(file_loc)//'. Exiting.')
            end if
 
            ! Further defining the structure of the formatted database slave
            ! file by describing in it the dimensionality of post-processed
            ! data as well as the total number of flow variable(s) that will
            ! eventually be stored in it
            if (output_partial_domain) then
                write (dbfile) x_output_idx%end - x_output_idx%beg, &
                    y_output_idx%end - y_output_idx%beg, &
                    z_output_idx%end - z_output_idx%beg, &
                    dbvars
            else
                write (dbfile) m, n, p, dbvars
            end if
 
            ! Next, analogous steps to the ones above are carried out by the
            ! root process to create and setup the formatted database master
            ! file. Note that this is only done in multidimensional cases.
            if (n == 0 .and. proc_rank == 0) then
 
                write (file_loc, '(A,I0,A)') '/', t_step, '.dat'
                file_loc = trim(rootdir)//trim(file_loc)
 
                open (dbroot, iostat=err, file=trim(file_loc), &
                      form='unformatted', status='replace')
 
                if (err /= 0) then
                    call s_mpi_abort('Unable to create Binary database '// &
                                     'master file '//trim(file_loc)// &
                                     '. Exiting.')
                end if
 
                if (output_partial_domain) then
                    write (dbroot) x_output_idx%end - x_output_idx%beg, 0, 0, dbvars
                else
                    write (dbroot) m_root, 0, 0, dbvars
                end if
 
            end if
 
        end if
 
    end subroutine s_open_formatted_database_file
 
    !> @brief Open the interface data file for appending extracted interface coordinates.
    impure subroutine s_open_intf_data_file()
 
        character(LEN=path_len + 3*name_len) :: file_path !<
              !! Relative path to a file in the case directory
 
        write (file_path, '(A)') '/intf_data.dat'
        file_path = trim(case_dir)//trim(file_path)
 
        ! Opening the simulation data file
        open (211, file=trim(file_path), &
              form='formatted', &
              position='append', &
              status='unknown')
 
    end subroutine s_open_intf_data_file
 
    !> @brief Open the energy data file for appending volume-integrated energy budget quantities.
    impure subroutine s_open_energy_data_file()
 
        character(LEN=path_len + 3*name_len) :: file_path !<
              !! Relative path to a file in the case directory
 
        write (file_path, '(A)') '/eng_data.dat'
        file_path = trim(case_dir)//trim(file_path)
 
        ! Opening the simulation data file
        open (251, file=trim(file_path), &
              form='formatted', &
              position='append', &
              status='unknown')
 
    end subroutine s_open_energy_data_file
 
    !> @brief Write the computational grid (cell-boundary coordinates) to the formatted database slave and master files.
    impure 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
        !              write the grid data of the domain segment that it is in
        !              charge of to the formatted database slave file. The root
        !              process will additionally take care of linking that grid
        !              data in the formatted database master file. In the Silo-
        !              HDF5 database format, the spatial extents of each local
        !              process grid are also written to the master file. In the
        !              Binary format, note that no master file is maintained in
        !              multidimensions. Finally, in 1D, no grid data is written
        !              within this subroutine for the Silo-HDF5 format because
        !              curve objects rather than quadrilateral meshes are used.
        !              For curve objects, in contrast to the quadrilateral mesh
        !              objects, the grid data is included side by side with the
        !              flow variable data. Then, in this case, we take care of
        !              writing both the grid and the flow variable data in the
        !              subroutine s_write_variable_to_formatted_database_file.
 
        ! Time-step that is currently being post-processed
        integer, intent(IN) :: t_step
 
        ! Bookkeeping variables storing the name and type of mesh that is
        ! handled by the local processor(s). Note that due to an internal
        ! NAG Fortran compiler problem, these two variables could not be
        ! allocated dynamically.
        character(LEN=4*name_len), dimension(num_procs) :: meshnames
        integer, dimension(num_procs) :: meshtypes
 
        ! Generic loop iterator
        integer :: i
 
        integer :: ierr !< Generic flag used to identify and report database errors
 
        ! Silo-HDF5 Database Format
 
        if (format == 1) then
 
            ! For multidimensional data sets, the spatial extents of all of
            ! the grid(s) handled by the local processor(s) are recorded so
            ! that they may be written, by root processor, to the formatted
            ! database master file.
            if (num_procs > 1) then
                call s_mpi_gather_spatial_extents(spatial_extents)
            elseif (p > 0) then
                if (grid_geometry == 3) then
                    spatial_extents(:, 0) = (/minval(y_cb), minval(z_cb), &
                                              minval(x_cb), maxval(y_cb), &
                                              maxval(z_cb), maxval(x_cb)/)
                else
                    spatial_extents(:, 0) = (/minval(x_cb), minval(y_cb), &
                                              minval(z_cb), maxval(x_cb), &
                                              maxval(y_cb), maxval(z_cb)/)
                end if
            elseif (n > 0) then
                spatial_extents(:, 0) = (/minval(x_cb), minval(y_cb), &
                                          maxval(x_cb), maxval(y_cb)/)
            else
                spatial_extents(:, 0) = (/minval(x_cb), maxval(x_cb)/)
            end if
 
            ! Next, the root processor proceeds to record all of the spatial
            ! extents in the formatted database master file. In addition, it
            ! also records a sub-domain connectivity map so that the entire
            ! grid may be reassembled by looking at the master file.
            if (proc_rank == 0) then
 
                do i = 1, num_procs
                    write (meshnames(i), '(A,I0,A,I0,A)') '../p', i - 1, &
                        '/', t_step, '.silo:rectilinear_grid'
                end do
 
                meshtypes = db_quad_rect
 
                err = dbset2dstrlen(len(meshnames(1)))
                err = dbmkoptlist(2, optlist)
                err = dbaddiopt(optlist, dbopt_extents_size, &
                                size(spatial_extents, 1))
                err = dbadddopt(optlist, dbopt_extents, spatial_extents)
                err = dbputmmesh(dbroot, 'rectilinear_grid', 16, &
                                 num_procs, meshnames, &
                                 len_trim(meshnames), &
                                 meshtypes, optlist, ierr)
                err = dbfreeoptlist(optlist)
 
            end if
 
            ! Finally, the local quadrilateral mesh, either 2D or 3D, along
            ! with its offsets that indicate the presence and size of ghost
            ! zone layer(s), are put in the formatted database slave file.
 
            if (p > 0) then
                err = dbmkoptlist(2, optlist)
                err = dbaddiopt(optlist, dbopt_lo_offset, lo_offset)
                err = dbaddiopt(optlist, dbopt_hi_offset, hi_offset)
                if (grid_geometry == 3) then
                    err = dbputqm(dbfile, 'rectilinear_grid', 16, &
                                  'x', 1, 'y', 1, 'z', 1, &
                                  y_cb, z_cb, x_cb, dims, 3, &
                                  db_double, db_collinear, &
                                  optlist, ierr)
                else
                    err = dbputqm(dbfile, 'rectilinear_grid', 16, &
                                  'x', 1, 'y', 1, 'z', 1, &
                                  x_cb, y_cb, z_cb, dims, 3, &
                                  db_double, db_collinear, &
                                  optlist, ierr)
                end if
                err = dbfreeoptlist(optlist)
            elseif (n > 0) then
                err = dbmkoptlist(2, optlist)
                err = dbaddiopt(optlist, dbopt_lo_offset, lo_offset)
                err = dbaddiopt(optlist, dbopt_hi_offset, hi_offset)
                err = dbputqm(dbfile, 'rectilinear_grid', 16, &
                              'x', 1, 'y', 1, 'z', 1, &
                              x_cb, y_cb, db_f77null, dims, 2, &
                              db_double, db_collinear, &
                              optlist, ierr)
                err = dbfreeoptlist(optlist)
            else
                err = dbmkoptlist(2, optlist)
                err = dbaddiopt(optlist, dbopt_lo_offset, lo_offset)
                err = dbaddiopt(optlist, dbopt_hi_offset, hi_offset)
                err = dbputqm(dbfile, 'rectilinear_grid', 16, &
                              'x', 1, 'y', 1, 'z', 1, &
                              x_cb, db_f77null, db_f77null, dims, 1, &
                              db_double, db_collinear, &
                              optlist, ierr)
                err = dbfreeoptlist(optlist)
            end if
            ! END: Silo-HDF5 Database Format
 
            ! Binary Database Format
 
        elseif (format == 2) then
 
            ! Multidimensional local grid data is written to the formatted
            ! database slave file. Recall that no master file to maintained
            ! in multidimensions.
            if (p > 0) then
                if (precision == 1) then
                    write (dbfile) real(x_cb, sp), &
                        real(y_cb, sp), &
                        real(z_cb, sp)
                else
                    if (output_partial_domain) then
                        write (dbfile) x_cb(x_output_idx%beg - 1:x_output_idx%end), &
                            y_cb(y_output_idx%beg - 1:y_output_idx%end), &
                            z_cb(z_output_idx%beg - 1:z_output_idx%end)
                    else
                        write (dbfile) x_cb, y_cb, z_cb
                    end if
                end if
 
            elseif (n > 0) then
                if (precision == 1) then
                    write (dbfile) real(x_cb, sp), &
                        real(y_cb, sp)
                else
                    if (output_partial_domain) then
                        write (dbfile) x_cb(x_output_idx%beg - 1:x_output_idx%end), &
                            y_cb(y_output_idx%beg - 1:y_output_idx%end)
                    else
                        write (dbfile) x_cb, y_cb
                    end if
                end if
 
                ! One-dimensional local grid data is written to the formatted
                ! database slave file. In addition, the local grid data is put
                ! together by the root process and written to the master file.
            else
 
                if (precision == 1) then
                    write (dbfile) real(x_cb, sp)
                else
                    if (output_partial_domain) then
                        write (dbfile) x_cb(x_output_idx%beg - 1:x_output_idx%end)
                    else
                        write (dbfile) x_cb
                    end if
                end if
 
                if (num_procs > 1) then
                    call s_mpi_defragment_1d_grid_variable()
                else
                    x_root_cb(:) = x_cb(:)
                end if
 
                if (proc_rank == 0) then
                    if (precision == 1) then
                        write (dbroot) real(x_root_cb, wp)
                    else
                        if (output_partial_domain) then
                            write (dbroot) x_root_cb(x_output_idx%beg - 1:x_output_idx%end)
                        else
                            write (dbroot) x_root_cb
                        end if
                    end if
                end if
 
            end if
 
        end if
 
    end subroutine s_write_grid_to_formatted_database_file
 
    !> @brief Write a single flow variable field to the formatted database slave and master files for a given time step.
    impure subroutine s_write_variable_to_formatted_database_file(varname, t_step)
        ! Description: The goal of this subroutine is to write to the formatted
        !              database file the flow variable at the current time-step,
        !              t_step. The local process(es) write the part of the flow
        !              variable that they handle to the formatted database slave
        !              file. The root process, on the other hand, will also take
        !              care of connecting all of the flow variable data in the
        !              formatted database master file. In the Silo-HDF5 database
        !              format, the extents of each local process flow variable
        !              are also written to the master file. Note that in Binary
        !              format, no master file is maintained in multidimensions.
        !              Finally note that in 1D, grid data is also written within
        !              this subroutine for Silo-HDF5 database format since curve
        !              and not the quadrilateral variable objects are used, see
        !              description of s_write_grid_to_formatted_database_file
        !              for more details on this topic.
 
        ! Name of the flow variable, which will be written to the formatted
        ! database file at the current time-step, t_step
        character(LEN=*), intent(IN) :: varname
 
        ! Time-step that is currently being post-processed
        integer, intent(IN) :: t_step
 
        ! Bookkeeping variables storing the name and type of flow variable
        ! that is about to be handled by the local processor(s). Note that
        ! due to an internal NAG Fortran compiler problem, these variables
        ! could not be allocated dynamically.
        character(LEN=4*name_len), dimension(num_procs) :: varnames
        integer, dimension(num_procs) :: vartypes
 
        ! Generic loop iterator
        integer :: i, j, k
 
        integer :: ierr !< Generic flag used to identify and report database errors
 
        ! Silo-HDF5 Database Format
 
        if (format == 1) then
 
            ! Determining the extents of the flow variable on each local
            ! process and gathering all this information on root process
            if (num_procs > 1) then
                call s_mpi_gather_data_extents(q_sf, data_extents)
            else
                data_extents(:, 0) = (/minval(q_sf), maxval(q_sf)/)
            end if
 
            ! Next, the root process proceeds to write the gathered flow
            ! variable data extents to formatted database master file.
            if (proc_rank == 0) then
 
                do i = 1, num_procs
                    write (varnames(i), '(A,I0,A,I0,A)') '../p', i - 1, &
                        '/', t_step, '.silo:'//trim(varname)
                end do
 
                vartypes = db_quadvar
 
                err = dbset2dstrlen(len(varnames(1)))
                err = dbmkoptlist(2, optlist)
                err = dbaddiopt(optlist, dbopt_extents_size, 2)
                err = dbadddopt(optlist, dbopt_extents, data_extents)
                err = dbputmvar(dbroot, trim(varname), &
                                len_trim(varname), num_procs, &
                                varnames, len_trim(varnames), &
                                vartypes, optlist, ierr)
                err = dbfreeoptlist(optlist)
 
            end if
 
            ! Finally, each of the local processor(s) proceeds to write
            ! the flow variable data that it is responsible for to the
            ! formatted database slave file.
            if (wp == dp) then
                if (precision == 1) then
                    do i = -offset_x%beg, m + offset_x%end
                        do j = -offset_y%beg, n + offset_y%end
                            do k = -offset_z%beg, p + offset_z%end
                                q_sf_s(i, j, k) = real(q_sf(i, j, k), sp)
                            end do
                        end do
                    end do
                    if (grid_geometry == 3) then
                        do i = -offset_x%beg, m + offset_x%end
                            do j = -offset_y%beg, n + offset_y%end
                                do k = -offset_z%beg, p + offset_z%end
                                    cyl_q_sf_s(j, k, i) = q_sf_s(i, j, k)
                                end do
                            end do
                        end do
                    end if
                else
                    if (grid_geometry == 3) then
                        do i = -offset_x%beg, m + offset_x%end
                            do j = -offset_y%beg, n + offset_y%end
                                do k = -offset_z%beg, p + offset_z%end
                                    cyl_q_sf(j, k, i) = q_sf(i, j, k)
                                end do
                            end do
                        end do
                    end if
                end if
            elseif (wp == sp) then
                do i = -offset_x%beg, m + offset_x%end
                    do j = -offset_y%beg, n + offset_y%end
                        do k = -offset_z%beg, p + offset_z%end
                            q_sf_s(i, j, k) = q_sf(i, j, k)
                        end do
                    end do
                end do
                if (grid_geometry == 3) then
                    do i = -offset_x%beg, m + offset_x%end
                        do j = -offset_y%beg, n + offset_y%end
                            do k = -offset_z%beg, p + offset_z%end
                                cyl_q_sf_s(j, k, i) = q_sf_s(i, j, k)
                            end do
                        end do
                    end do
                end if
            end if
 
# 972 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                if (precision == 1) then
                    if (p > 0) then
                        if (grid_geometry == 3) then
                            err = dbputqv1(dbfile, trim(varname), &
                                           len_trim(varname), &
                                           'rectilinear_grid', 16, &
                                           cyl_q_sf_s, dims - 1, 3, db_f77null, &
                                           0, db_float, db_zonecent, &
                                           db_f77null, ierr)
                        else
                            err = dbputqv1(dbfile, trim(varname), &
                                           len_trim(varname), &
                                           'rectilinear_grid', 16, &
                                           q_sf_s, dims - 1, 3, db_f77null, &
                                           0, db_float, db_zonecent, &
                                           db_f77null, ierr)
                        end if
                    elseif (n > 0) then
                        err = dbputqv1(dbfile, trim(varname), &
                                       len_trim(varname), &
                                       'rectilinear_grid', 16, &
                                       q_sf_s, dims - 1, 2, db_f77null, &
                                       0, db_float, db_zonecent, &
                                       db_f77null, ierr)
                    else
                        err = dbputqv1(dbfile, trim(varname), &
                                       len_trim(varname), &
                                       'rectilinear_grid', 16, &
                                       q_sf_s, dims - 1, 1, db_f77null, &
                                       0, db_float, db_zonecent, &
                                       db_f77null, ierr)
 
                    end if
                end if
# 972 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                if (precision == 2) then
                    if (p > 0) then
                        if (grid_geometry == 3) then
                            err = dbputqv1(dbfile, trim(varname), &
                                           len_trim(varname), &
                                           'rectilinear_grid', 16, &
                                           cyl_q_sf, dims - 1, 3, db_f77null, &
                                           0, db_double, db_zonecent, &
                                           db_f77null, ierr)
                        else
                            err = dbputqv1(dbfile, trim(varname), &
                                           len_trim(varname), &
                                           'rectilinear_grid', 16, &
                                           q_sf, dims - 1, 3, db_f77null, &
                                           0, db_double, db_zonecent, &
                                           db_f77null, ierr)
                        end if
                    elseif (n > 0) then
                        err = dbputqv1(dbfile, trim(varname), &
                                       len_trim(varname), &
                                       'rectilinear_grid', 16, &
                                       q_sf, dims - 1, 2, db_f77null, &
                                       0, db_double, db_zonecent, &
                                       db_f77null, ierr)
                    else
                        err = dbputqv1(dbfile, trim(varname), &
                                       len_trim(varname), &
                                       'rectilinear_grid', 16, &
                                       q_sf, dims - 1, 1, db_f77null, &
                                       0, db_double, db_zonecent, &
                                       db_f77null, ierr)
 
                    end if
                end if
# 1007 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
 
            ! END: Silo-HDF5 Database Format
 
            ! Binary Database Format
 
        else
 
            ! Writing the name of the flow variable and its data, associated
            ! with the local processor, to the formatted database slave file
            if (precision == 1) then
                write (dbfile) varname, real(q_sf, wp)
            else
                write (dbfile) varname, q_sf
            end if
 
            ! In 1D, the root process also takes care of gathering the flow
            ! variable data from all of the local processor(s) and writes it
            ! to the formatted database master file.
            if (n == 0) then
 
                if (num_procs > 1) then
                    call s_mpi_defragment_1d_flow_variable(q_sf, q_root_sf)
                else
                    q_root_sf(:, :, :) = q_sf(:, :, :)
                end if
 
                if (proc_rank == 0) then
                    if (precision == 1) then
                        write (dbroot) varname, real(q_root_sf, wp)
                    else
                        write (dbroot) varname, q_root_sf
                    end if
                end if
 
            end if
 
        end if
 
    end subroutine s_write_variable_to_formatted_database_file
 
    !>  Subroutine that writes the post processed results in the folder 'lag_bubbles_data'
            !!  @param t_step Current time step
    impure subroutine s_write_lag_bubbles_results_to_text(t_step)
 
        integer, intent(in) :: t_step
 
        character(len=len_trim(case_dir) + 3*name_len) :: file_loc
 
        integer :: id
 
#ifdef MFC_MPI
        real(wp), dimension(20) :: inputvals
        real(wp) :: time_real
        integer, dimension(MPI_STATUS_SIZE) :: status
        integer(KIND=MPI_OFFSET_KIND) :: disp
        integer :: view
 
        logical :: lg_bub_file, file_exist
 
        integer, dimension(2) :: gsizes, lsizes, start_idx_part
        integer :: ifile
        integer :: ierr !< Generic flag used to identify and report MPI errors
        real(wp) :: file_time, file_dt
        integer :: file_num_procs, file_tot_part, tot_part
        integer :: i
 
        integer, dimension(:), allocatable :: proc_bubble_counts
        real(wp), dimension(1:1, 1:lag_io_vars) :: lag_io_null
        lag_io_null = 0._wp
 
        ! Construct file path
        write (file_loc, '(A,I0,A)') 'lag_bubbles_', t_step, '.dat'
        file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc)
 
        ! Check if file exists
        inquire (file=trim(file_loc), exist=file_exist)
        if (.not. file_exist) then
            call s_mpi_abort('Restart file '//trim(file_loc)//' does not exist!')
        end if
 
        if (.not. parallel_io) return
 
        if (proc_rank == 0) then
            call mpi_file_open(mpi_comm_self, file_loc, mpi_mode_rdonly, &
                               mpi_info_int, ifile, ierr)
 
            call mpi_file_read(ifile, file_tot_part, 1, mpi_integer, status, ierr)
            call mpi_file_read(ifile, file_time, 1, mpi_p, status, ierr)
            call mpi_file_read(ifile, file_dt, 1, mpi_p, status, ierr)
            call mpi_file_read(ifile, file_num_procs, 1, mpi_integer, status, ierr)
 
            call mpi_file_close(ifile, ierr)
        end if
 
        call mpi_bcast(file_tot_part, 1, mpi_integer, 0, mpi_comm_world, ierr)
        call mpi_bcast(file_time, 1, mpi_p, 0, mpi_comm_world, ierr)
        call mpi_bcast(file_dt, 1, mpi_p, 0, mpi_comm_world, ierr)
        call mpi_bcast(file_num_procs, 1, mpi_integer, 0, mpi_comm_world, ierr)
 
        allocate (proc_bubble_counts(file_num_procs))
 
        if (proc_rank == 0) then
            call mpi_file_open(mpi_comm_self, file_loc, mpi_mode_rdonly, &
                               mpi_info_int, ifile, ierr)
 
            ! Skip to processor counts position
            disp = int(sizeof(file_tot_part) + 2*sizeof(file_time) + sizeof(file_num_procs), &
                       mpi_offset_kind)
            call mpi_file_seek(ifile, disp, mpi_seek_set, ierr)
            call mpi_file_read(ifile, proc_bubble_counts, file_num_procs, mpi_integer, status, ierr)
 
            call mpi_file_close(ifile, ierr)
        end if
 
        call mpi_bcast(proc_bubble_counts, file_num_procs, mpi_integer, 0, mpi_comm_world, ierr)
 
        gsizes(1) = file_tot_part
        gsizes(2) = lag_io_vars
        lsizes(1) = file_tot_part
        lsizes(2) = lag_io_vars
        start_idx_part(1) = 0
        start_idx_part(2) = 0
 
        call mpi_type_create_subarray(2, gsizes, lsizes, start_idx_part, &
                                      mpi_order_fortran, mpi_p, view, ierr)
        call mpi_type_commit(view, ierr)
 
        call mpi_file_open(mpi_comm_world, file_loc, mpi_mode_rdonly, &
                           mpi_info_int, ifile, ierr)
 
        disp = int(sizeof(file_tot_part) + 2*sizeof(file_time) + sizeof(file_num_procs) + &
                   file_num_procs*sizeof(proc_bubble_counts(1)), mpi_offset_kind)
        call mpi_file_set_view(ifile, disp, mpi_p, view, &
                               'native', mpi_info_null, ierr)
 
        allocate (mpi_io_data_lg_bubbles(file_tot_part, 1:lag_io_vars))
 
        call mpi_file_read_all(ifile, mpi_io_data_lg_bubbles, lag_io_vars*file_tot_part, &
                               mpi_p, status, ierr)
 
        write (file_loc, '(A,I0,A)') 'lag_bubbles_post_process_', t_step, '.dat'
        file_loc = trim(case_dir)//'/lag_bubbles_post_process/'//trim(file_loc)
 
        if (proc_rank == 0) then
            open (unit=29, file=file_loc, form='formatted', position='rewind')
 
            if (lag_header) then
                write (29, '(A)', advance='no')
                if (lag_id_wrt) write (29, '(A8)', advance='no') 'id, '
                if (lag_pos_wrt) write (29, '(3(A17))', advance='no') 'px, ', 'py, ', 'pz, '
                if (lag_pos_prev_wrt) write (29, '(3(A17))', advance='no') 'pvx, ', 'pvy, ', 'pvz, '
                if (lag_vel_wrt) write (29, '(3(A17))', advance='no') 'vx, ', 'vy, ', 'vz, '
                if (lag_rad_wrt) write (29, '(A17)', advance='no') 'radius, '
                if (lag_rvel_wrt) write (29, '(A17)', advance='no') 'rvel, '
                if (lag_r0_wrt) write (29, '(A17)', advance='no') 'r0, '
                if (lag_rmax_wrt) write (29, '(A17)', advance='no') 'rmax, '
                if (lag_rmin_wrt) write (29, '(A17)', advance='no') 'rmin, '
                if (lag_dphidt_wrt) write (29, '(A17)', advance='no') 'dphidt, '
                if (lag_pres_wrt) write (29, '(A17)', advance='no') 'pressure, '
                if (lag_mv_wrt) write (29, '(A17)', advance='no') 'mv, '
                if (lag_mg_wrt) write (29, '(A17)', advance='no') 'mg, '
                if (lag_betat_wrt) write (29, '(A17)', advance='no') 'betaT, '
                if (lag_betac_wrt) write (29, '(A17)', advance='no') 'betaC, '
                write (29, '(A15)') 'time'
            end if
 
            do i = 1, file_tot_part
                id = int(mpi_io_data_lg_bubbles(i, 1))
                inputvals(1:20) = mpi_io_data_lg_bubbles(i, 2:21)
                if (id > 0) then
                    write (29, '(100(A))', advance='no') ''
                    if (lag_id_wrt) write (29, '(I6, A)', advance='no') id, ', '
                    if (lag_pos_wrt) write (29, '(3(E15.7, A))', advance='no') inputvals(1), ', ', inputvals(2), ', ', inputvals(3), ', '
                    if (lag_pos_prev_wrt) write (29, '(3(E15.7, A))', advance='no') inputvals(4), ', ', inputvals(5), ', ', inputvals(6), ', '
                    if (lag_vel_wrt) write (29, '(3(E15.7, A))', advance='no') inputvals(7), ', ', inputvals(8), ', ', inputvals(9), ', '
                    if (lag_rad_wrt) write (29, '(E15.7, A)', advance='no') inputvals(10), ', '
                    if (lag_rvel_wrt) write (29, '(E15.7, A)', advance='no') inputvals(11), ', '
                    if (lag_r0_wrt) write (29, '(E15.7, A)', advance='no') inputvals(12), ', '
                    if (lag_rmax_wrt) write (29, '(E15.7, A)', advance='no') inputvals(13), ', '
                    if (lag_rmin_wrt) write (29, '(E15.7, A)', advance='no') inputvals(14), ', '
                    if (lag_dphidt_wrt) write (29, '(E15.7, A)', advance='no') inputvals(15), ', '
                    if (lag_pres_wrt) write (29, '(E15.7, A)', advance='no') inputvals(16), ', '
                    if (lag_mv_wrt) write (29, '(E15.7, A)', advance='no') inputvals(17), ', '
                    if (lag_mg_wrt) write (29, '(E15.7, A)', advance='no') inputvals(18), ', '
                    if (lag_betat_wrt) write (29, '(E15.7, A)', advance='no') inputvals(19), ', '
                    if (lag_betac_wrt) write (29, '(E15.7, A)', advance='no') inputvals(20), ', '
                    write (29, '(E15.7)') time_real
                end if
            end do
            close (29)
        end if
 
        deallocate (mpi_io_data_lg_bubbles)
 
        call s_mpi_barrier()
 
        call mpi_file_close(ifile, ierr)
#endif
 
    end subroutine s_write_lag_bubbles_results_to_text
 
    !> @brief Read Lagrangian bubble restart data and write bubble positions and scalar fields to the Silo database.
    impure subroutine s_write_lag_bubbles_to_formatted_database_file(t_step)
 
        integer, intent(in) :: t_step
 
        character(len=len_trim(case_dir) + 3*name_len) :: file_loc
 
        integer :: id
 
#ifdef MFC_MPI
        real(wp), dimension(20) :: inputvals
        real(wp) :: time_real
        integer, dimension(MPI_STATUS_SIZE) :: status
        integer(KIND=MPI_OFFSET_KIND) :: disp
        integer :: view
 
        logical :: lg_bub_file, file_exist
 
        integer, dimension(2) :: gsizes, lsizes, start_idx_part
        integer :: ifile, ierr, tot_data, valid_data, nbub
        real(wp) :: file_time, file_dt
        integer :: file_num_procs, file_tot_part
        integer, dimension(:), allocatable :: proc_bubble_counts
        real(wp), dimension(1:1, 1:lag_io_vars) :: dummy
        character(LEN=4*name_len), dimension(num_procs) :: meshnames
        integer, dimension(num_procs) :: meshtypes
        real(wp) :: dummy_data
 
        integer :: i, j
 
        real(wp), dimension(:), allocatable :: bub_id
        real(wp), dimension(:), allocatable :: px, py, pz, ppx, ppy, ppz, vx, vy, vz
        real(wp), dimension(:), allocatable :: radius, rvel, rnot, rmax, rmin, dphidt
        real(wp), dimension(:), allocatable :: pressure, mv, mg, betat, betac
 
        dummy = 0._wp
        dummy_data = 0._wp
 
        ! Construct file path
        write (file_loc, '(A,I0,A)') 'lag_bubbles_', t_step, '.dat'
        file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc)
 
        ! Check if file exists
        inquire (file=trim(file_loc), exist=file_exist)
        if (.not. file_exist) then
            call s_mpi_abort('Restart file '//trim(file_loc)//' does not exist!')
        end if
 
        if (.not. parallel_io) return
 
        if (proc_rank == 0) then
            call mpi_file_open(mpi_comm_self, file_loc, mpi_mode_rdonly, &
                               mpi_info_int, ifile, ierr)
 
            call mpi_file_read(ifile, file_tot_part, 1, mpi_integer, status, ierr)
            call mpi_file_read(ifile, file_time, 1, mpi_p, status, ierr)
            call mpi_file_read(ifile, file_dt, 1, mpi_p, status, ierr)
            call mpi_file_read(ifile, file_num_procs, 1, mpi_integer, status, ierr)
 
            call mpi_file_close(ifile, ierr)
        end if
 
        call mpi_bcast(file_tot_part, 1, mpi_integer, 0, mpi_comm_world, ierr)
        call mpi_bcast(file_time, 1, mpi_p, 0, mpi_comm_world, ierr)
        call mpi_bcast(file_dt, 1, mpi_p, 0, mpi_comm_world, ierr)
        call mpi_bcast(file_num_procs, 1, mpi_integer, 0, mpi_comm_world, ierr)
 
        allocate (proc_bubble_counts(file_num_procs))
 
        if (proc_rank == 0) then
            call mpi_file_open(mpi_comm_self, file_loc, mpi_mode_rdonly, &
                               mpi_info_int, ifile, ierr)
 
            ! Skip to processor counts position
            disp = int(sizeof(file_tot_part) + 2*sizeof(file_time) + sizeof(file_num_procs), &
                       mpi_offset_kind)
            call mpi_file_seek(ifile, disp, mpi_seek_set, ierr)
            call mpi_file_read(ifile, proc_bubble_counts, file_num_procs, mpi_integer, status, ierr)
 
            call mpi_file_close(ifile, ierr)
        end if
 
        call mpi_bcast(proc_bubble_counts, file_num_procs, mpi_integer, 0, mpi_comm_world, ierr)
 
        ! Set time variables from file
 
        nbub = proc_bubble_counts(proc_rank + 1)
 
        start_idx_part(1) = 0
        do i = 1, proc_rank
            start_idx_part(1) = start_idx_part(1) + proc_bubble_counts(i)
        end do
 
        start_idx_part(2) = 0
        lsizes(1) = nbub
        lsizes(2) = lag_io_vars
 
        gsizes(1) = file_tot_part
        gsizes(2) = lag_io_vars
 
        if (nbub > 0) then
 
# 1313 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (bub_id(nbub))
# 1313 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (px(nbub))
# 1313 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (py(nbub))
# 1313 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (pz(nbub))
# 1313 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (ppx(nbub))
# 1313 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (ppy(nbub))
# 1313 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (ppz(nbub))
# 1313 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (vx(nbub))
# 1313 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (vy(nbub))
# 1313 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (vz(nbub))
# 1313 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (radius(nbub))
# 1313 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (rvel(nbub))
# 1313 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (rnot(nbub))
# 1313 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (rmax(nbub))
# 1313 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (rmin(nbub))
# 1313 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (dphidt(nbub))
# 1313 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (pressure(nbub))
# 1313 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (mv(nbub))
# 1313 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (mg(nbub))
# 1313 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (betat(nbub))
# 1313 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (betac(nbub))
# 1315 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
            allocate (mpi_io_data_lg_bubbles(nbub, 1:lag_io_vars))
 
            call mpi_type_create_subarray(2, gsizes, lsizes, start_idx_part, &
                                          mpi_order_fortran, mpi_p, view, ierr)
            call mpi_type_commit(view, ierr)
 
            call mpi_file_open(mpi_comm_world, file_loc, mpi_mode_rdonly, &
                               mpi_info_int, ifile, ierr)
 
            ! Skip extended header
            disp = int(sizeof(file_tot_part) + 2*sizeof(file_time) + sizeof(file_num_procs) + &
                       file_num_procs*sizeof(proc_bubble_counts(1)), mpi_offset_kind)
            call mpi_file_set_view(ifile, disp, mpi_p, view, 'native', mpi_info_int, ierr)
 
            call mpi_file_read_all(ifile, mpi_io_data_lg_bubbles, &
                                   lag_io_vars*nbub, mpi_p, status, ierr)
 
            call mpi_file_close(ifile, ierr)
            call mpi_type_free(view, ierr)
 
            ! Extract data from MPI_IO_DATA_lg_bubbles array
            ! Adjust these indices based on your actual data layout
# 1341 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                bub_id(:) = mpi_io_data_lg_bubbles(:, 1)
# 1341 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                px(:) = mpi_io_data_lg_bubbles(:, 2)
# 1341 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                py(:) = mpi_io_data_lg_bubbles(:, 3)
# 1341 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                pz(:) = mpi_io_data_lg_bubbles(:, 4)
# 1341 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                ppx(:) = mpi_io_data_lg_bubbles(:, 5)
# 1341 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                ppy(:) = mpi_io_data_lg_bubbles(:, 6)
# 1341 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                ppz(:) = mpi_io_data_lg_bubbles(:, 7)
# 1341 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                vx(:) = mpi_io_data_lg_bubbles(:, 8)
# 1341 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                vy(:) = mpi_io_data_lg_bubbles(:, 9)
# 1341 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                vz(:) = mpi_io_data_lg_bubbles(:, 10)
# 1341 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                radius(:) = mpi_io_data_lg_bubbles(:, 11)
# 1341 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                rvel(:) = mpi_io_data_lg_bubbles(:, 12)
# 1341 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                rnot(:) = mpi_io_data_lg_bubbles(:, 13)
# 1341 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                rmax(:) = mpi_io_data_lg_bubbles(:, 14)
# 1341 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                rmin(:) = mpi_io_data_lg_bubbles(:, 15)
# 1341 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                dphidt(:) = mpi_io_data_lg_bubbles(:, 16)
# 1341 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                pressure(:) = mpi_io_data_lg_bubbles(:, 17)
# 1341 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                mv(:) = mpi_io_data_lg_bubbles(:, 18)
# 1341 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                mg(:) = mpi_io_data_lg_bubbles(:, 19)
# 1341 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                betat(:) = mpi_io_data_lg_bubbles(:, 20)
# 1341 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                betac(:) = mpi_io_data_lg_bubbles(:, 21)
# 1343 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
 
            ! Next, the root processor proceeds to record all of the spatial
            ! extents in the formatted database master file. In addition, it
            ! also records a sub-domain connectivity map so that the entire
            ! grid may be reassembled by looking at the master file.
            if (proc_rank == 0) then
 
                do i = 1, num_procs
                    write (meshnames(i), '(A,I0,A,I0,A)') '../p', i - 1, &
                        '/', t_step, '.silo:lag_bubbles'
                    meshtypes(i) = db_pointmesh
                end do
                err = dbset2dstrlen(len(meshnames(1)))
                err = dbputmmesh(dbroot, 'lag_bubbles', 16, &
                                 num_procs, meshnames, &
                                 len_trim(meshnames), &
                                 meshtypes, db_f77null, ierr)
            end if
 
            err = dbputpm(dbfile, 'lag_bubbles', 11, 3, &
                          px, py, pz, nbub, &
                          db_double, db_f77null, ierr)
 
            if (lag_id_wrt) call s_write_lag_variable_to_formatted_database_file('part_id', t_step, bub_id, nbub)
            if (lag_vel_wrt) then
                call s_write_lag_variable_to_formatted_database_file('part_vel1', t_step, vx, nbub)
                call s_write_lag_variable_to_formatted_database_file('part_vel2', t_step, vy, nbub)
                if (p > 0) call s_write_lag_variable_to_formatted_database_file('part_vel3', t_step, vz, nbub)
            end if
            if (lag_rad_wrt) call s_write_lag_variable_to_formatted_database_file('part_radius', t_step, radius, nbub)
            if (lag_rvel_wrt) call s_write_lag_variable_to_formatted_database_file('part_rdot', t_step, rvel, nbub)
            if (lag_r0_wrt) call s_write_lag_variable_to_formatted_database_file('part_r0', t_step, rnot, nbub)
            if (lag_rmax_wrt) call s_write_lag_variable_to_formatted_database_file('part_rmax', t_step, rmax, nbub)
            if (lag_rmin_wrt) call s_write_lag_variable_to_formatted_database_file('part_rmin', t_step, rmin, nbub)
            if (lag_dphidt_wrt) call s_write_lag_variable_to_formatted_database_file('part_dphidt', t_step, dphidt, nbub)
            if (lag_pres_wrt) call s_write_lag_variable_to_formatted_database_file('part_pressure', t_step, pressure, nbub)
            if (lag_mv_wrt) call s_write_lag_variable_to_formatted_database_file('part_mv', t_step, mv, nbub)
            if (lag_mg_wrt) call s_write_lag_variable_to_formatted_database_file('part_mg', t_step, mg, nbub)
            if (lag_betat_wrt) call s_write_lag_variable_to_formatted_database_file('part_betaT', t_step, betat, nbub)
            if (lag_betac_wrt) call s_write_lag_variable_to_formatted_database_file('part_betaC', t_step, betac, nbub)
 
            deallocate (bub_id, px, py, pz, ppx, ppy, ppz, vx, vy, vz, radius, &
                        rvel, rnot, rmax, rmin, dphidt, pressure, mv, mg, &
                        betat, betac)
            deallocate (mpi_io_data_lg_bubbles)
        else
            call mpi_type_contiguous(0, mpi_p, view, ierr)
            call mpi_type_commit(view, ierr)
 
            call mpi_file_open(mpi_comm_world, file_loc, mpi_mode_rdonly, &
                               mpi_info_int, ifile, ierr)
 
            ! Skip extended header
            disp = int(sizeof(file_tot_part) + 2*sizeof(file_time) + sizeof(file_num_procs) + &
                       file_num_procs*sizeof(proc_bubble_counts(1)), mpi_offset_kind)
            call mpi_file_set_view(ifile, disp, mpi_p, view, 'native', mpi_info_int, ierr)
 
            call mpi_file_read_all(ifile, dummy, 0, mpi_p, status, ierr)
 
            call mpi_file_close(ifile, ierr)
            call mpi_type_free(view, ierr)
 
            if (proc_rank == 0) then
 
                do i = 1, num_procs
                    write (meshnames(i), '(A,I0,A,I0,A)') '../p', i - 1, &
                        '/', t_step, '.silo:lag_bubbles'
                    meshtypes(i) = db_pointmesh
                end do
                err = dbset2dstrlen(len(meshnames(1)))
                err = dbputmmesh(dbroot, 'lag_bubbles', 16, &
                                 num_procs, meshnames, &
                                 len_trim(meshnames), &
                                 meshtypes, db_f77null, ierr)
            end if
 
            err = dbsetemptyok(1)
            err = dbputpm(dbfile, 'lag_bubbles', 11, 3, &
                          dummy_data, dummy_data, dummy_data, 0, &
                          db_double, db_f77null, ierr)
 
            if (lag_id_wrt) call s_write_lag_variable_to_formatted_database_file('part_id', t_step)
            if (lag_vel_wrt) then
                call s_write_lag_variable_to_formatted_database_file('part_vel1', t_step)
                call s_write_lag_variable_to_formatted_database_file('part_vel2', t_step)
                if (p > 0) call s_write_lag_variable_to_formatted_database_file('part_vel3', t_step)
            end if
            if (lag_rad_wrt) call s_write_lag_variable_to_formatted_database_file('part_radius', t_step)
            if (lag_rvel_wrt) call s_write_lag_variable_to_formatted_database_file('part_rdot', t_step)
            if (lag_r0_wrt) call s_write_lag_variable_to_formatted_database_file('part_r0', t_step)
            if (lag_rmax_wrt) call s_write_lag_variable_to_formatted_database_file('part_rmax', t_step)
            if (lag_rmin_wrt) call s_write_lag_variable_to_formatted_database_file('part_rmin', t_step)
            if (lag_dphidt_wrt) call s_write_lag_variable_to_formatted_database_file('part_dphidt', t_step)
            if (lag_pres_wrt) call s_write_lag_variable_to_formatted_database_file('part_pressure', t_step)
            if (lag_mv_wrt) call s_write_lag_variable_to_formatted_database_file('part_mv', t_step)
            if (lag_mg_wrt) call s_write_lag_variable_to_formatted_database_file('part_mg', t_step)
            if (lag_betat_wrt) call s_write_lag_variable_to_formatted_database_file('part_betaT', t_step)
            if (lag_betac_wrt) call s_write_lag_variable_to_formatted_database_file('part_betaC', t_step)
        end if
 
#endif
 
    end subroutine s_write_lag_bubbles_to_formatted_database_file
 
    !> @brief Write a single Lagrangian bubble point-variable to the Silo database slave and master files.
    subroutine s_write_lag_variable_to_formatted_database_file(varname, t_step, data, nBubs)
 
        character(len=*), intent(in) :: varname
        integer, intent(in) :: t_step
        real(wp), dimension(1:), intent(in), optional :: data
        integer, intent(in), optional :: nBubs
 
        character(len=64), dimension(num_procs) :: var_names
        integer, dimension(num_procs) :: var_types
        real(wp) :: dummy_data
 
        integer :: ierr !< Generic flag used to identify and report database errors
        integer :: i
 
        dummy_data = 0._wp
 
        if (present(nbubs) .and. present(data)) then
            if (proc_rank == 0) then
                do i = 1, num_procs
                    write (var_names(i), '(A,I0,A,I0,A)') '../p', i - 1, &
                        '/', t_step, '.silo:'//trim(varname)
                    var_types(i) = db_pointvar
                end do
                err = dbset2dstrlen(len(var_names(1)))
                err = dbputmvar(dbroot, trim(varname), len_trim(varname), &
                                num_procs, var_names, &
                                len_trim(var_names), &
                                var_types, db_f77null, ierr)
            end if
 
            err = dbputpv1(dbfile, trim(varname), len_trim(varname), &
                           'lag_bubbles', 11, data, nbubs, db_double, db_f77null, ierr)
        else
            if (proc_rank == 0) then
                do i = 1, num_procs
                    write (var_names(i), '(A,I0,A,I0,A)') '../p', i - 1, &
                        '/', t_step, '.silo:'//trim(varname)
                    var_types(i) = db_pointvar
                end do
                err = dbset2dstrlen(len(var_names(1)))
                err = dbsetemptyok(1)
                err = dbputmvar(dbroot, trim(varname), len_trim(varname), &
                                num_procs, var_names, &
                                len_trim(var_names), &
                                var_types, db_f77null, ierr)
            end if
 
            err = dbsetemptyok(1)
            err = dbputpv1(dbfile, trim(varname), len_trim(varname), &
                           'lag_bubbles', 11, dummy_data, 0, db_double, db_f77null, ierr)
        end if
 
    end subroutine s_write_lag_variable_to_formatted_database_file
 
    !> @brief Extract the volume-fraction interface contour from primitive fields and write the coordinates to the interface data file.
    impure subroutine s_write_intf_data_file(q_prim_vf)
 
        type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf
        integer :: i, j, k, l, cent !< Generic loop iterators
        integer :: counter, root !< number of data points extracted to fit shape to SH perturbations
        real(wp), allocatable :: x_td(:), y_td(:), x_d1(:), y_d1(:), y_d(:), x_d(:)
        real(wp) :: axp, axm, ayp, aym, tgp, euc_d, thres, maxalph_loc, maxalph_glb
 
        allocate (x_d1(m*n))
        allocate (y_d1(m*n))
        counter = 0
        maxalph_loc = 0._wp
        do k = 0, p
            do j = 0, n
                do i = 0, m
                    if (q_prim_vf(e_idx + 2)%sf(i, j, k) > maxalph_loc) then
                        maxalph_loc = q_prim_vf(e_idx + 2)%sf(i, j, k)
                    end if
                end do
            end do
        end do
 
        call s_mpi_allreduce_max(maxalph_loc, maxalph_glb)
        if (p > 0) then
            do l = 0, p
                if (z_cc(l) < dz(l) .and. z_cc(l) > 0) then
                    cent = l
                end if
            end do
        else
            cent = 0
        end if
 
        thres = 0.9_wp*maxalph_glb
        do k = 0, n
            do j = 0, m
                axp = q_prim_vf(e_idx + 2)%sf(j + 1, k, cent)
                axm = q_prim_vf(e_idx + 2)%sf(j, k, cent)
                ayp = q_prim_vf(e_idx + 2)%sf(j, k + 1, cent)
                aym = q_prim_vf(e_idx + 2)%sf(j, k, cent)
                if ((axp > thres .and. axm < thres) .or. (axp < thres .and. axm > thres) &
                    .or. (ayp > thres .and. aym < thres) .or. (ayp < thres .and. aym > thres)) then
                    if (counter == 0) then
                        counter = counter + 1
                        x_d1(counter) = x_cc(j)
                        y_d1(counter) = y_cc(k)
                        euc_d = sqrt((x_cc(j) - x_d1(i))**2 + (y_cc(k) - y_d1(i))**2)
                        tgp = sqrt(dx(j)**2 + dy(k)**2)
                    else
                        euc_d = sqrt((x_cc(j) - x_d1(i))**2 + (y_cc(k) - y_d1(i))**2)
                        tgp = sqrt(dx(j)**2 + dy(k)**2)
                        do i = 1, counter
                            if (euc_d < tgp) then
                                cycle
                            elseif (euc_d > tgp .and. i == counter) then
                                counter = counter + 1
                                x_d1(counter) = x_cc(j)
                                y_d1(counter) = y_cc(k)
 
                            end if
                        end do
                    end if
                end if
            end do
        end do
 
        allocate (x_d(counter), y_d(counter))
 
        do i = 1, counter
            y_d(i) = y_d1(i)
            x_d(i) = x_d1(i)
        end do
        root = 0
 
        call s_mpi_gather_data(x_d, counter, x_td, root)
        call s_mpi_gather_data(y_d, counter, y_td, root)
        if (proc_rank == 0) then
            do i = 1, size(x_td)
                if (i == size(x_td)) then
                    write (211, '(F12.9,1X,F12.9,1X,I4)') &
                        x_td(i), y_td(i), size(x_td)
                else
                    write (211, '(F12.9,1X,F12.9,1X,F3.1)') &
                        x_td(i), y_td(i), 0._wp
                end if
            end do
        end if
 
    end subroutine s_write_intf_data_file
 
    !> @brief Compute volume-integrated kinetic, potential, and internal energies and write the energy budget to the energy data file.
    impure subroutine s_write_energy_data_file(q_prim_vf, q_cons_vf)
        type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf, q_cons_vf
        real(wp) :: elk, egk, elp, egint, vb, vl, pres_av, et
        real(wp) :: rho, pres, dv, tmp, gamma, pi_inf, maxma, maxma_glb, maxvel, c, ma, h, qv
        real(wp), dimension(num_vels) :: vel
        real(wp), dimension(num_fluids) :: adv
        integer :: i, j, k, l, s !looping indices
 
        egk = 0._wp
        elp = 0._wp
        egint = 0._wp
        vb = 0._wp
        maxvel = 0._wp
        maxma = 0._wp
        vl = 0._wp
        elk = 0._wp
        et = 0._wp
        vb = 0._wp
        dv = 0._wp
        pres_av = 0._wp
        pres = 0._wp
        c = 0._wp
 
        do k = 0, p
            do j = 0, n
                do i = 0, m
                    pres = 0._wp
                    dv = dx(i)*dy(j)*dz(k)
                    rho = 0._wp
                    gamma = 0._wp
                    pi_inf = 0._wp
                    qv = 0._wp
                    pres = q_prim_vf(e_idx)%sf(i, j, k)
                    egint = egint + q_prim_vf(e_idx + 2)%sf(i, j, k)*(gammas(2)*pres)*dv
                    do s = 1, num_vels
                        vel(s) = q_prim_vf(num_fluids + s)%sf(i, j, k)
                        egk = egk + 0.5_wp*q_prim_vf(e_idx + 2)%sf(i, j, k)*q_prim_vf(2)%sf(i, j, k)*vel(s)*vel(s)*dv
                        elk = elk + 0.5_wp*q_prim_vf(e_idx + 1)%sf(i, j, k)*q_prim_vf(1)%sf(i, j, k)*vel(s)*vel(s)*dv
                        if (abs(vel(s)) > maxvel) then
                            maxvel = abs(vel(s))
                        end if
                    end do
                    do l = 1, adv_idx%end - e_idx
                        adv(l) = q_prim_vf(e_idx + l)%sf(i, j, k)
                        gamma = gamma + adv(l)*gammas(l)
                        pi_inf = pi_inf + adv(l)*pi_infs(l)
                        rho = rho + adv(l)*q_prim_vf(l)%sf(i, j, k)
                        qv = qv + adv(l)*q_prim_vf(l)%sf(i, j, k)*qvs(l)
                    end do
 
                    h = ((gamma + 1._wp)*pres + pi_inf + qv)/rho
 
                    call s_compute_speed_of_sound(pres, rho, &
                                                  gamma, pi_inf, &
                                                  h, adv, 0._wp, 0._wp, c, qv)
 
                    ma = maxvel/c
                    if (ma > maxma .and. (adv(1) > (1.0_wp - 1.0e-10_wp))) then
                        maxma = ma
                    end if
                    vl = vl + adv(1)*dv
                    vb = vb + adv(2)*dv
                    pres_av = pres_av + adv(1)*pres*dv
                    et = et + q_cons_vf(e_idx)%sf(i, j, k)*dv
                end do
            end do
        end do
 
        tmp = pres_av
        call s_mpi_allreduce_sum(tmp, pres_av)
        tmp = vl
        call s_mpi_allreduce_sum(tmp, vl)
 
        call s_mpi_allreduce_max(maxma, maxma_glb)
        tmp = elk
        call s_mpi_allreduce_sum(tmp, elk)
        tmp = egint
        call s_mpi_allreduce_sum(tmp, egint)
        tmp = egk
        call s_mpi_allreduce_sum(tmp, egk)
        tmp = vb
        call s_mpi_allreduce_sum(tmp, vb)
        tmp = et
        call s_mpi_allreduce_sum(tmp, et)
 
        elp = pres_av/vl*vb
        if (proc_rank == 0) then
            write (251, '(10X, 8F24.8)') &
                elp, &
                egint, &
                elk, &
                egk, &
                et, &
                vb, &
                vl, &
                maxma_glb
        end if
 
    end subroutine s_write_energy_data_file
 
    !> @brief Close the formatted database slave file and, for the root process, the master file.
    impure 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
        !              typically close two files, one associated with the local
        !              sub-domain and the other with the entire domain. The non-
        !              root process(es) must close one file, which is associated
        !              with the local sub-domain. Note that for the Binary data-
        !              base format and multidimensional data, the root process
        !              only has to close the file associated with the local sub-
        !              domain, because one associated with the entire domain is
        !              not generated.
 
        integer :: ierr !< Generic flag used to identify and report database errors
 
        ! Silo-HDF5 database format
        if (format == 1) then
            ierr = dbclose(dbfile)
            if (proc_rank == 0) ierr = dbclose(dbroot)
 
            ! Binary database format
        else
            close (dbfile)
            if (n == 0 .and. proc_rank == 0) close (dbroot)
 
        end if
 
    end subroutine s_close_formatted_database_file
 
    !> @brief Close the interface data file.
    impure subroutine s_close_intf_data_file()
 
        close (211)
 
    end subroutine s_close_intf_data_file
 
    !> @brief Close the energy data file.
    impure subroutine s_close_energy_data_file()
 
        close (251)
 
    end subroutine s_close_energy_data_file
 
    !> @brief Deallocate module arrays and release all data-output resources.
    impure subroutine s_finalize_data_output_module()
        ! Description: Deallocation procedures for the module
 
        ! Deallocating the generic storage employed for the flow variable(s)
        ! that were written to the formatted database file(s). Note that the
        ! root variable is only deallocated in the case of a 1D computation.
        deallocate (q_sf)
        if (n == 0) deallocate (q_root_sf)
        if (grid_geometry == 3) then
            deallocate (cyl_q_sf)
        end if
 
        ! Deallocating spatial and data extents and also the variables for
        ! the offsets and the one bookkeeping the number of cell-boundaries
        ! in each active coordinate direction. Note that all these variables
        ! were only needed by Silo-HDF5 format for multidimensional data.
        if (format == 1) then
            deallocate (spatial_extents)
            deallocate (data_extents)
            deallocate (lo_offset)
            deallocate (hi_offset)
            deallocate (dims)
        end if
 
    end subroutine s_finalize_data_output_module
 
end module m_data_output