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
    use m_global_parameters
    use m_derived_variables
    use m_mpi_proxy
    use m_compile_specific
    use m_helper
    use m_variables_conversion
    use m_constants, only: model_eqns_gamma_law, model_eqns_5eq, model_eqns_6eq, format_silo, format_binary, precision_single
 
    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_ib_state_files, s_write_intf_data_file, &
        & s_write_energy_data_file, s_write_ib_bodies_to_formatted_database_file, s_close_formatted_database_file, &
        & s_close_intf_data_file, s_close_energy_data_file, s_finalize_data_output_module, out
 
    ! Include Silo-HDF5 interface library
    include 'silo_f9x.inc'
 
    !> Output workspace: flow variable buffers, VisIt extents/offsets, directory paths, file handles, and variable count.
    type(output_context) :: out
 
    ! Generic error flag for Silo-HDF5 and Binary I/O operations
    integer, private :: err
 
contains
 
    !> Allocate storage arrays, configure output directories, and count flow variables for formatted database output.
    impure subroutine s_initialize_data_output_module()
 
        character(LEN=len_trim(case_dir) + 2*name_len) :: file_loc
        logical                                        :: dir_check
        integer                                        :: i
 
        allocate (out%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 (out%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 == precision_single) then
            allocate (out%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 (out%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 (out%q_root_sf(0:m_root,0:0,0:0))
            if (precision == precision_single) then
                allocate (out%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 == format_silo) then
            allocate (out%data_extents(1:2,0:num_procs - 1))
 
            if (p > 0) then
                allocate (out%spatial_extents(1:6,0:num_procs - 1))
                allocate (out%lo_offset(1:3))
                allocate (out%hi_offset(1:3))
                allocate (out%dims(1:3))
            else if (n > 0) then
                allocate (out%spatial_extents(1:4,0:num_procs - 1))
                allocate (out%lo_offset(1:2))
                allocate (out%hi_offset(1:2))
                allocate (out%dims(1:2))
            else
                allocate (out%spatial_extents(1:2,0:num_procs - 1))
                allocate (out%lo_offset(1:1))
                allocate (out%hi_offset(1:1))
                allocate (out%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 == format_silo) then
            if (p > 0) then
                if (grid_geometry == 3) then
                    out%lo_offset(:) = (/offset_y%beg, offset_z%beg, offset_x%beg/)
                    out%hi_offset(:) = (/offset_y%end, offset_z%end, offset_x%end/)
                else
                    out%lo_offset(:) = (/offset_x%beg, offset_y%beg, offset_z%beg/)
                    out%hi_offset(:) = (/offset_x%end, offset_y%end, offset_z%end/)
                end if
 
                if (grid_geometry == 3) then
                    out%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
                    out%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
            else if (n > 0) then
                out%lo_offset(:) = (/offset_x%beg, offset_y%beg/)
                out%hi_offset(:) = (/offset_x%end, offset_y%end/)
 
                out%dims(:) = (/m + offset_x%beg + offset_x%end + 2, n + offset_y%beg + offset_y%end + 2/)
            else
                out%lo_offset(:) = (/offset_x%beg/)
                out%hi_offset(:) = (/offset_x%end/)
                out%dims(:) = (/m + offset_x%beg + offset_x%end + 2/)
            end if
        end if
 
        if (format == format_silo) then
            out%dbdir = trim(case_dir) // '/silo_hdf5'
 
            write (out%proc_rank_dir, '(A,I0)') '/p', proc_rank
 
            out%proc_rank_dir = trim(out%dbdir) // trim(out%proc_rank_dir)
 
            file_loc = trim(out%proc_rank_dir) // '/.'
 
            call my_inquire(file_loc, dir_check)
            if (dir_check .neqv. .true.) then
                call s_create_directory(trim(out%proc_rank_dir))
            end if
 
            if (proc_rank == 0) then
                out%rootdir = trim(out%dbdir) // '/root'
 
                file_loc = trim(out%rootdir) // '/.'
 
                call my_inquire(file_loc, dir_check)
                if (dir_check .neqv. .true.) then
                    call s_create_directory(trim(out%rootdir))
                end if
            end if
        else
            out%dbdir = trim(case_dir) // '/binary'
 
            write (out%proc_rank_dir, '(A,I0)') '/p', proc_rank
 
            out%proc_rank_dir = trim(out%dbdir) // trim(out%proc_rank_dir)
 
            file_loc = trim(out%proc_rank_dir) // '/.'
 
            call my_inquire(file_loc, dir_check)
 
            if (dir_check .neqv. .true.) then
                call s_create_directory(trim(out%proc_rank_dir))
            end if
 
            if (n == 0 .and. proc_rank == 0) then
                out%rootdir = trim(out%dbdir) // '/root'
 
                file_loc = trim(out%rootdir) // '/.'
 
                call my_inquire(file_loc, dir_check)
 
                if (dir_check .neqv. .true.) then
                    call s_create_directory(trim(out%rootdir))
                end if
            end if
        end if
 
        if (bubbles_lagrange) then  ! Lagrangian solver
            if (lag_txt_wrt) then
                out%dbdir = trim(case_dir) // '/lag_bubbles_post_process'
                file_loc = trim(out%dbdir) // '/.'
                call my_inquire(file_loc, dir_check)
 
                if (dir_check .neqv. .true.) then
                    call s_create_directory(trim(out%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 == format_binary) then
            if (n == 0 .and. proc_rank == 0) out%dbroot = 2
            out%dbfile = 1
        end if
 
        if (format == format_binary) then
            out%dbvars = 0
 
            if ((model_eqns == model_eqns_5eq) .or. (model_eqns == model_eqns_6eq)) then
                do i = 1, num_fluids
                    if (alpha_rho_wrt(i) .or. (cons_vars_wrt .or. prim_vars_wrt)) then
                        out%dbvars = out%dbvars + 1
                    end if
                end do
            end if
 
            if ((rho_wrt .or. (model_eqns == model_eqns_gamma_law .and. (cons_vars_wrt .or. prim_vars_wrt))) &
                & .and. (.not. relativity)) then
                out%dbvars = out%dbvars + 1
            end if
 
            if (relativity .and. (rho_wrt .or. prim_vars_wrt)) out%dbvars = out%dbvars + 1
            if (relativity .and. (rho_wrt .or. cons_vars_wrt)) out%dbvars = out%dbvars + 1
 
            do i = 1, eqn_idx%E - eqn_idx%mom%beg
                if (mom_wrt(i) .or. cons_vars_wrt) out%dbvars = out%dbvars + 1
            end do
 
            do i = 1, eqn_idx%E - eqn_idx%mom%beg
                if (vel_wrt(i) .or. prim_vars_wrt) out%dbvars = out%dbvars + 1
            end do
 
            do i = 1, eqn_idx%E - eqn_idx%mom%beg
                if (flux_wrt(i)) out%dbvars = out%dbvars + 1
            end do
 
            if (e_wrt .or. cons_vars_wrt) out%dbvars = out%dbvars + 1
            if (pres_wrt .or. prim_vars_wrt) out%dbvars = out%dbvars + 1
            if (hypoelasticity) out%dbvars = out%dbvars + (num_dims*(num_dims + 1))/2
            if (cont_damage) out%dbvars = out%dbvars + 1
            if (hyper_cleaning) out%dbvars = out%dbvars + 1
 
            if (mhd) then
                if (n == 0) then
                    out%dbvars = out%dbvars + 2
                else
                    out%dbvars = out%dbvars + 3
                end if
            end if
 
            if ((model_eqns == model_eqns_5eq) .or. (model_eqns == model_eqns_6eq)) then
                do i = 1, num_fluids - 1
                    if (alpha_wrt(i) .or. (cons_vars_wrt .or. prim_vars_wrt)) then
                        out%dbvars = out%dbvars + 1
                    end if
                end do
 
                if (alpha_wrt(num_fluids) .or. (cons_vars_wrt .or. prim_vars_wrt)) then
                    out%dbvars = out%dbvars + 1
                end if
            end if
 
            if (gamma_wrt .or. (model_eqns == model_eqns_gamma_law .and. (cons_vars_wrt .or. prim_vars_wrt))) then
                out%dbvars = out%dbvars + 1
            end if
 
            if (heat_ratio_wrt) out%dbvars = out%dbvars + 1
 
            if (pi_inf_wrt .or. (model_eqns == model_eqns_gamma_law .and. (cons_vars_wrt .or. prim_vars_wrt))) then
                out%dbvars = out%dbvars + 1
            end if
 
            if (pres_inf_wrt) out%dbvars = out%dbvars + 1
            if (c_wrt) out%dbvars = out%dbvars + 1
 
            if (p > 0) then
                do i = 1, num_vels
                    if (omega_wrt(i)) out%dbvars = out%dbvars + 1
                end do
            else if (n > 0) then
                do i = 1, num_vels
                    if (omega_wrt(i)) out%dbvars = out%dbvars + 1
                end do
            end if
 
            if (schlieren_wrt) out%dbvars = out%dbvars + 1
        end if
 
    end subroutine s_initialize_data_output_module
 
    !> 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
 
# 285 "/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
# 285 "/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
# 285 "/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
# 310 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
 
    end subroutine s_define_output_region
 
    !> 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)
 
        integer, intent(in)                            :: t_step
        character(LEN=len_trim(case_dir) + 3*name_len) :: file_loc
        integer                                        :: ierr
 
        if (format == format_silo) then
            write (file_loc, '(A,I0,A)') '/', t_step, '.silo'
            file_loc = trim(out%proc_rank_dir) // trim(file_loc)
 
            ierr = dbcreate(trim(file_loc), len_trim(file_loc), db_clobber, db_local, 'MFC v3.0', 8, db_hdf5, out%dbfile)
 
            if (out%dbfile == -1) then
                call s_mpi_abort('Unable to create Silo-HDF5 database ' // 'slave file ' // trim(file_loc) // '. ' // 'Exiting.')
            end if
 
            if (proc_rank == 0) then
                write (file_loc, '(A,I0,A)') '/collection_', t_step, '.silo'
                file_loc = trim(out%rootdir) // trim(file_loc)
 
                ierr = dbcreate(trim(file_loc), len_trim(file_loc), db_clobber, db_local, 'MFC v3.0', 8, db_hdf5, out%dbroot)
 
                if (out%dbroot == -1) then
                    call s_mpi_abort('Unable to create Silo-HDF5 database ' // 'master file ' // trim(file_loc) // '. ' &
                                     & // 'Exiting.')
                end if
            end if
        else
            write (file_loc, '(A,I0,A)') '/', t_step, '.dat'
            file_loc = trim(out%proc_rank_dir) // trim(file_loc)
 
            open (out%dbfile, iostat=err, file=trim(file_loc), form='unformatted', status='replace')
 
            if (err /= 0) then
                call s_mpi_abort('Unable to create Binary database slave ' // 'file ' // trim(file_loc) // '. Exiting.')
            end if
 
            if (output_partial_domain) then
                write (out%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, out%dbvars
            else
                write (out%dbfile) m, n, p, out%dbvars
            end if
 
            if (n == 0 .and. proc_rank == 0) then
                write (file_loc, '(A,I0,A)') '/', t_step, '.dat'
                file_loc = trim(out%rootdir) // trim(file_loc)
 
                open (out%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 (out%dbroot) x_output_idx%end - x_output_idx%beg, 0, 0, out%dbvars
                else
                    write (out%dbroot) m_root, 0, 0, out%dbvars
                end if
            end if
        end if
 
    end subroutine s_open_formatted_database_file
 
    !> 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
 
        write (file_path, '(A)') '/intf_data.dat'
        file_path = trim(case_dir) // trim(file_path)
 
        open (211, file=trim(file_path), form='formatted', position='append', status='unknown')
 
    end subroutine s_open_intf_data_file
 
    !> 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
 
        write (file_path, '(A)') '/eng_data.dat'
        file_path = trim(case_dir) // trim(file_path)
 
        open (251, file=trim(file_path), form='formatted', position='append', status='unknown')
 
    end subroutine s_open_energy_data_file
 
    !> 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)
 
        integer, intent(in) :: t_step
 
        ! NAG compiler requires these to be statically sized
        character(LEN=4*name_len), dimension(num_procs) :: meshnames
        integer, dimension(num_procs)                   :: meshtypes
        integer                                         :: i
        integer                                         :: ierr
 
        if (format == format_silo) 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(out%spatial_extents)
            else if (p > 0) then
                if (grid_geometry == 3) then
                    out%spatial_extents(:,0) = (/minval(y_cb), minval(z_cb), minval(x_cb), maxval(y_cb), maxval(z_cb), &
                                        & maxval(x_cb)/)
                else
                    out%spatial_extents(:,0) = (/minval(x_cb), minval(y_cb), minval(z_cb), maxval(x_cb), maxval(y_cb), &
                                        & maxval(z_cb)/)
                end if
            else if (n > 0) then
                out%spatial_extents(:,0) = (/minval(x_cb), minval(y_cb), maxval(x_cb), maxval(y_cb)/)
            else
                out%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, out%optlist)
                err = dbaddiopt(out%optlist, dbopt_extents_size, size(out%spatial_extents, 1))
                err = dbadddopt(out%optlist, dbopt_extents, out%spatial_extents)
                err = dbputmmesh(out%dbroot, 'rectilinear_grid', 16, num_procs, meshnames, len_trim(meshnames), meshtypes, &
                                 & out%optlist, ierr)
                err = dbfreeoptlist(out%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, out%optlist)
                err = dbaddiaopt(out%optlist, dbopt_lo_offset, size(out%lo_offset), out%lo_offset)
                err = dbaddiaopt(out%optlist, dbopt_hi_offset, size(out%hi_offset), out%hi_offset)
                if (grid_geometry == 3) then
                    err = dbputqm(out%dbfile, 'rectilinear_grid', 16, 'x', 1, 'y', 1, 'z', 1, y_cb, z_cb, x_cb, out%dims, 3, &
                                  & db_double, db_collinear, out%optlist, ierr)
                else
                    err = dbputqm(out%dbfile, 'rectilinear_grid', 16, 'x', 1, 'y', 1, 'z', 1, x_cb, y_cb, z_cb, out%dims, 3, &
                                  & db_double, db_collinear, out%optlist, ierr)
                end if
                err = dbfreeoptlist(out%optlist)
            else if (n > 0) then
                err = dbmkoptlist(2, out%optlist)
                err = dbaddiaopt(out%optlist, dbopt_lo_offset, size(out%lo_offset), out%lo_offset)
                err = dbaddiaopt(out%optlist, dbopt_hi_offset, size(out%hi_offset), out%hi_offset)
                err = dbputqm(out%dbfile, 'rectilinear_grid', 16, 'x', 1, 'y', 1, 'z', 1, x_cb, y_cb, db_f77null, out%dims, 2, &
                              & db_double, db_collinear, out%optlist, ierr)
                err = dbfreeoptlist(out%optlist)
            else
                err = dbmkoptlist(2, out%optlist)
                err = dbaddiaopt(out%optlist, dbopt_lo_offset, size(out%lo_offset), out%lo_offset)
                err = dbaddiaopt(out%optlist, dbopt_hi_offset, size(out%hi_offset), out%hi_offset)
                err = dbputqm(out%dbfile, 'rectilinear_grid', 16, 'x', 1, 'y', 1, 'z', 1, x_cb, db_f77null, db_f77null, out%dims, &
                              & 1, db_double, db_collinear, out%optlist, ierr)
                err = dbfreeoptlist(out%optlist)
            end if
        else if (format == format_binary) 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 == precision_single) then
                    write (out%dbfile) real(x_cb, sp), real(y_cb, sp), real(z_cb, sp)
                else
                    if (output_partial_domain) then
                        write (out%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 (out%dbfile) x_cb, y_cb, z_cb
                    end if
                end if
            else if (n > 0) then
                if (precision == precision_single) then
                    write (out%dbfile) real(x_cb, sp), real(y_cb, sp)
                else
                    if (output_partial_domain) then
                        write (out%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 (out%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 == precision_single) then
                    write (out%dbfile) real(x_cb, sp)
                else
                    if (output_partial_domain) then
                        write (out%dbfile) x_cb(x_output_idx%beg - 1:x_output_idx%end)
                    else
                        write (out%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 == precision_single) then
                        write (out%dbroot) real(x_root_cb, wp)
                    else
                        if (output_partial_domain) then
                            write (out%dbroot) x_root_cb(x_output_idx%beg - 1:x_output_idx%end)
                        else
                            write (out%dbroot) x_root_cb
                        end if
                    end if
                end if
            end if
        end if
 
    end subroutine s_write_grid_to_formatted_database_file
 
    !> 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)
 
        character(LEN=*), intent(in) :: varname
        integer, intent(in)          :: t_step
 
        ! NAG compiler requires these to be statically sized
        character(LEN=4*name_len), dimension(num_procs) :: varnames
        integer, dimension(num_procs)                   :: vartypes
        integer                                         :: i, j, k
        integer                                         :: ierr
 
        if (format == format_silo) 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(out%q_sf, out%data_extents)
            else
                out%data_extents(:,0) = (/minval(out%q_sf), maxval(out%q_sf)/)
            end if
 
            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, out%optlist)
                err = dbaddiopt(out%optlist, dbopt_extents_size, 2)
                err = dbadddopt(out%optlist, dbopt_extents, out%data_extents)
                err = dbputmvar(out%dbroot, trim(varname), len_trim(varname), num_procs, varnames, len_trim(varnames), vartypes, &
                                & out%optlist, ierr)
                err = dbfreeoptlist(out%optlist)
            end if
 
            if (wp == dp) then
                if (precision == precision_single) 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
                                out%q_sf_s(i, j, k) = real(out%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
                                    out%cyl_q_sf_s(j, k, i) = out%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
                                    out%cyl_q_sf(j, k, i) = out%q_sf(i, j, k)
                                end do
                            end do
                        end do
                    end if
                end if
            else if (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
                            out%q_sf_s(i, j, k) = out%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
                                out%cyl_q_sf_s(j, k, i) = out%q_sf_s(i, j, k)
                            end do
                        end do
                    end do
                end if
            end if
 
# 626 "/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(out%dbfile, trim(varname), len_trim(varname), 'rectilinear_grid', 16, &
                                           & out%cyl_q_sf_s, out%dims - 1, 3, db_f77null, 0, db_float, db_zonecent, &
                                           & db_f77null, ierr)
                        else
                            err = dbputqv1(out%dbfile, trim(varname), len_trim(varname), 'rectilinear_grid', 16, out%q_sf_s, &
                                           & out%dims - 1, 3, db_f77null, 0, db_float, db_zonecent, db_f77null, ierr)
                        end if
                    else if (n > 0) then
                        err = dbputqv1(out%dbfile, trim(varname), len_trim(varname), 'rectilinear_grid', 16, out%q_sf_s, &
                                       & out%dims - 1, 2, db_f77null, 0, db_float, db_zonecent, db_f77null, ierr)
                    else
                        err = dbputqv1(out%dbfile, trim(varname), len_trim(varname), 'rectilinear_grid', 16, out%q_sf_s, &
                                       & out%dims - 1, 1, db_f77null, 0, db_float, db_zonecent, db_f77null, ierr)
                    end if
                end if
# 626 "/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(out%dbfile, trim(varname), len_trim(varname), 'rectilinear_grid', 16, &
                                           & out%cyl_q_sf, out%dims - 1, 3, db_f77null, 0, db_double, db_zonecent, &
                                           & db_f77null, ierr)
                        else
                            err = dbputqv1(out%dbfile, trim(varname), len_trim(varname), 'rectilinear_grid', 16, out%q_sf, &
                                           & out%dims - 1, 3, db_f77null, 0, db_double, db_zonecent, db_f77null, ierr)
                        end if
                    else if (n > 0) then
                        err = dbputqv1(out%dbfile, trim(varname), len_trim(varname), 'rectilinear_grid', 16, out%q_sf, &
                                       & out%dims - 1, 2, db_f77null, 0, db_double, db_zonecent, db_f77null, ierr)
                    else
                        err = dbputqv1(out%dbfile, trim(varname), len_trim(varname), 'rectilinear_grid', 16, out%q_sf, &
                                       & out%dims - 1, 1, db_f77null, 0, db_double, db_zonecent, db_f77null, ierr)
                    end if
                end if
# 645 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
        else
            ! Writing the name of the flow variable and its data, associated with the local processor, to the formatted database
            ! slave file
            if (precision == precision_single) then
                write (out%dbfile) varname, real(out%q_sf, wp)
            else
                write (out%dbfile) varname, out%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(out%q_sf, out%q_root_sf)
                else
                    out%q_root_sf(:,:,:) = out%q_sf(:,:,:)
                end if
 
                if (proc_rank == 0) then
                    if (precision == precision_single) then
                        write (out%dbroot) varname, real(out%q_root_sf, wp)
                    else
                        write (out%dbroot) varname, out%q_root_sf
                    end if
                end if
            end if
        end if
 
    end subroutine s_write_variable_to_formatted_database_file
 
    !> Write the post-processed results in the folder 'lag_bubbles_data'
    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                                :: file_exist
        integer, dimension(2)                  :: gsizes, lsizes, start_idx_part
        integer                                :: ifile
        integer                                :: ierr
        real(wp)                               :: file_time, file_dt
        integer                                :: file_num_procs, file_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)
        time_real = file_time
 
        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
 
    !> 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)                                        :: time_real
        integer, dimension(MPI_STATUS_SIZE)             :: status
        integer(KIND=MPI_OFFSET_KIND)                   :: disp
        integer                                         :: view
        logical                                         :: file_exist
        integer, dimension(2)                           :: gsizes, lsizes, start_idx_part
        integer                                         :: ifile, ierr, 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
        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)
        time_real = file_time
 
        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
# 923 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (bub_id(nbub))
# 923 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (px(nbub))
# 923 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (py(nbub))
# 923 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (pz(nbub))
# 923 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (ppx(nbub))
# 923 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (ppy(nbub))
# 923 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (ppz(nbub))
# 923 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (vx(nbub))
# 923 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (vy(nbub))
# 923 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (vz(nbub))
# 923 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (radius(nbub))
# 923 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (rvel(nbub))
# 923 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (rnot(nbub))
# 923 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (rmax(nbub))
# 923 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (rmin(nbub))
# 923 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (dphidt(nbub))
# 923 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (pressure(nbub))
# 923 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (mv(nbub))
# 923 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (mg(nbub))
# 923 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (betat(nbub))
# 923 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                allocate (betac(nbub))
# 925 "/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
# 947 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                bub_id(:) = mpi_io_data_lg_bubbles(:,1)
# 947 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                px(:) = mpi_io_data_lg_bubbles(:,2)
# 947 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                py(:) = mpi_io_data_lg_bubbles(:,3)
# 947 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                pz(:) = mpi_io_data_lg_bubbles(:,4)
# 947 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                ppx(:) = mpi_io_data_lg_bubbles(:,5)
# 947 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                ppy(:) = mpi_io_data_lg_bubbles(:,6)
# 947 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                ppz(:) = mpi_io_data_lg_bubbles(:,7)
# 947 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                vx(:) = mpi_io_data_lg_bubbles(:,8)
# 947 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                vy(:) = mpi_io_data_lg_bubbles(:,9)
# 947 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                vz(:) = mpi_io_data_lg_bubbles(:,10)
# 947 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                radius(:) = mpi_io_data_lg_bubbles(:,11)
# 947 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                rvel(:) = mpi_io_data_lg_bubbles(:,12)
# 947 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                rnot(:) = mpi_io_data_lg_bubbles(:,13)
# 947 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                rmax(:) = mpi_io_data_lg_bubbles(:,14)
# 947 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                rmin(:) = mpi_io_data_lg_bubbles(:,15)
# 947 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                dphidt(:) = mpi_io_data_lg_bubbles(:,16)
# 947 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                pressure(:) = mpi_io_data_lg_bubbles(:,17)
# 947 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                mv(:) = mpi_io_data_lg_bubbles(:,18)
# 947 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                mg(:) = mpi_io_data_lg_bubbles(:,19)
# 947 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                betat(:) = mpi_io_data_lg_bubbles(:,20)
# 947 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
                betac(:) = mpi_io_data_lg_bubbles(:,21)
# 949 "/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(out%dbroot, 'lag_bubbles', 16, num_procs, meshnames, len_trim(meshnames), meshtypes, db_f77null, &
                                 & ierr)
            end if
 
            err = dbputpm(out%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(out%dbroot, 'lag_bubbles', 16, num_procs, meshnames, len_trim(meshnames), meshtypes, db_f77null, &
                                 & ierr)
            end if
 
            err = dbsetemptyok(1)
            err = dbputpm(out%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
 
    !> 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
        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(out%dbroot, trim(varname), len_trim(varname), num_procs, var_names, len_trim(var_names), &
                                & var_types, db_f77null, ierr)
            end if
 
            err = dbputpv1(out%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(out%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(out%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
 
    !> Convert the binary immersed-boundary state file to per-body formatted text files
    impure subroutine s_write_ib_state_files()
 
        character(len=len_trim(case_dir) + 4*name_len) :: in_file, out_file, file_loc
        integer                                        :: iu_in, ios, i, rec_id
        integer, allocatable, dimension(:)             :: iu_out
        real(wp)                                       :: rec_time
        real(wp), dimension(3)                         :: rec_force, rec_torque
        real(wp), dimension(3)                         :: rec_vel, rec_angular_vel
        real(wp), dimension(3)                         :: rec_angles, rec_centroid
 
        file_loc = trim(case_dir) // '/D'
 
        in_file = trim(file_loc) // '/ib_state.dat'
        open (newunit=iu_in, file=trim(in_file), form='unformatted', access='stream', status='old', action='read', iostat=ios)
        if (ios /= 0) then
            call s_mpi_abort('Cannot open IB state input file: ' // trim(in_file))
        end if
 
        allocate (iu_out(num_ibs))
        do i = 1, num_ibs
            write (out_file, '(A,I0,A)') trim(file_loc) // '/ib_', i, '.txt'
            open (newunit=iu_out(i), file=trim(out_file), form='formatted', status='replace', action='write', iostat=ios)
            if (ios /= 0) then
                call s_mpi_abort('Cannot open IB state output file: ' // trim(out_file))
            end if
            write (iu_out(i), &
                   & '(A)') 'mytime fx fy fz Tau_x Tau_y Tau_z vx vy vz omega_x omega_y omega_z angle_x angle_y angle_z x_c y_c z_c'
        end do
 
        do
            read (iu_in, iostat=ios) rec_time, rec_id, rec_force, rec_torque, rec_vel, rec_angular_vel, rec_angles, &
                  & rec_centroid(1), rec_centroid(2), rec_centroid(3)
            if (ios /= 0) exit
 
            if (rec_id >= 1 .and. rec_id <= num_ibs) then
                write (iu_out(rec_id), '(19(ES24.16E3,1X))') rec_time, rec_force(1), rec_force(2), rec_force(3), rec_torque(1), &
                       & rec_torque(2), rec_torque(3), rec_vel(1), rec_vel(2), rec_vel(3), rec_angular_vel(1), &
                       & rec_angular_vel(2), rec_angular_vel(3), rec_angles(1), rec_angles(2), rec_angles(3), rec_centroid(1), &
                       & rec_centroid(2), rec_centroid(3)
            end if
        end do
 
        close (iu_in)
        do i = 1, num_ibs
            close (iu_out(i))
        end do
        deallocate (iu_out)
 
    end subroutine s_write_ib_state_files
 
    !> 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
        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(eqn_idx%E + 2)%sf(i, j, k) > maxalph_loc) then
                        maxalph_loc = q_prim_vf(eqn_idx%E + 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(eqn_idx%E + 2)%sf(j + 1, k, cent)
                axm = q_prim_vf(eqn_idx%E + 2)%sf(j, k, cent)
                ayp = q_prim_vf(eqn_idx%E + 2)%sf(j, k + 1, cent)
                aym = q_prim_vf(eqn_idx%E + 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)
                    else
                        tgp = sqrt(dx(j)**2 + dy(k)**2)
                        do i = 1, counter
                            euc_d = sqrt((x_cc(j) - x_d1(i))**2 + (y_cc(k) - y_d1(i))**2)
                            if (euc_d < tgp) then
                                exit
                            else if (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
 
    !> 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(eqn_idx%E)%sf(i, j, k)
                    egint = egint + q_prim_vf(eqn_idx%E + 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(eqn_idx%E + 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(eqn_idx%E + 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, eqn_idx%adv%end - eqn_idx%E
                        adv(l) = q_prim_vf(eqn_idx%E + 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(eqn_idx%E)%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
 
    !> Read IB state and write a Silo point mesh with per-body scalar fields.
    impure subroutine s_write_ib_bodies_to_formatted_database_file(t_step)
 
        integer, intent(in)                            :: t_step
        character(len=len_trim(case_dir) + 3*name_len) :: file_loc
 
#ifdef MFC_MPI
        integer, parameter                              :: nfields_per_ib = 20
        real(wp)                                        :: ib_buf(nfields_per_ib)
        real(wp), dimension(:,:), allocatable           :: ib_data
        logical                                         :: file_exist
        character(LEN=4*name_len), dimension(num_procs) :: meshnames
        integer, dimension(num_procs)                   :: meshtypes
        integer                                         :: i, ios, file_unit
        integer                                         :: ierr, nbodies
        integer                                         :: r, nlocal, gbl_id
        character(len=10)                               :: t_step_string
        real(wp), dimension(:), allocatable             :: px, py, pz
        real(wp), dimension(:), allocatable             :: force_x, force_y, force_z
        real(wp), dimension(:), allocatable             :: torque_x, torque_y, torque_z
        real(wp), dimension(:), allocatable             :: vel_x, vel_y, vel_z
        real(wp), dimension(:), allocatable             :: omega_x, omega_y, omega_z
        real(wp), dimension(:), allocatable             :: angle_x, angle_y, angle_z
        real(wp), dimension(:), allocatable             :: ib_diameter
 
        if (proc_rank == 0) then
            nbodies = num_ibs
 
            if (nbodies > 0) then
                allocate (ib_data(nbodies, nfields_per_ib))
                allocate (px(nbodies), py(nbodies), pz(nbodies))
                allocate (force_x(nbodies), force_y(nbodies), force_z(nbodies))
                allocate (torque_x(nbodies), torque_y(nbodies), torque_z(nbodies))
                allocate (vel_x(nbodies), vel_y(nbodies), vel_z(nbodies))
                allocate (omega_x(nbodies), omega_y(nbodies), omega_z(nbodies))
                allocate (angle_x(nbodies), angle_y(nbodies), angle_z(nbodies))
                allocate (ib_diameter(nbodies))
 
                if (file_per_process) then
                    call s_int_to_str(t_step, t_step_string)
                    ib_data = 0._wp
                    do r = 0, num_procs - 1
                        write (file_loc, '(A,I0,A,i7.7,A)') 'ib_state_', t_step, '_', r, '.dat'
                        file_loc = trim(case_dir) // '/restart_data/lustre_' // trim(t_step_string) // '/' // trim(file_loc)
 
                        inquire (file=trim(file_loc), exist=file_exist)
                        if (.not. file_exist) call s_mpi_abort('Restart file ' // trim(file_loc) // ' does not exist!')
 
                        open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='old', iostat=ios)
                        if (ios /= 0) call s_mpi_abort('Cannot open IB state file: ' // trim(file_loc))
 
                        read (file_unit, iostat=ios) nlocal
                        if (ios /= 0) call s_mpi_abort('Error reading IB state file header: ' // trim(file_loc))
 
                        do i = 1, nlocal
                            read (file_unit, iostat=ios) gbl_id
                            if (ios /= 0) call s_mpi_abort('Error reading IB patch ID: ' // trim(file_loc))
                            read (file_unit, iostat=ios) ib_buf
                            if (ios /= 0) call s_mpi_abort('Error reading IB state data: ' // trim(file_loc))
                            ib_data(gbl_id,:) = ib_buf(:)
                        end do
 
                        close (file_unit)
                    end do
                else
                    ! Build path to per-timestep IB state file
                    write (file_loc, '(A,I0,A)') '/restart_data/ib_state_', t_step, '.dat'
                    file_loc = trim(case_dir) // trim(file_loc)
 
                    inquire (file=trim(file_loc), exist=file_exist)
                    if (.not. file_exist) call s_mpi_abort('Restart file ' // trim(file_loc) // ' does not exist!')
 
                    open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='old', iostat=ios)
                    if (ios /= 0) call s_mpi_abort('Cannot open IB state file: ' // trim(file_loc))
 
                    do i = 1, nbodies
                        read (file_unit, iostat=ios) ib_buf
                        if (ios /= 0) call s_mpi_abort('Error reading IB state file')
                        ib_data(i,:) = ib_buf(:)
                    end do
 
                    close (file_unit)
                end if
 
                do i = 1, nbodies
                    force_x(i) = ib_data(i, 2); force_y(i) = ib_data(i, 3); force_z(i) = ib_data(i, 4)
                    torque_x(i) = ib_data(i, 5); torque_y(i) = ib_data(i, 6); torque_z(i) = ib_data(i, 7)
                    vel_x(i) = ib_data(i, 8); vel_y(i) = ib_data(i, 9); vel_z(i) = ib_data(i, 10)
                    omega_x(i) = ib_data(i, 11); omega_y(i) = ib_data(i, 12); omega_z(i) = ib_data(i, 13)
                    angle_x(i) = ib_data(i, 14); angle_y(i) = ib_data(i, 15); angle_z(i) = ib_data(i, 16)
                    px(i) = ib_data(i, 17); py(i) = ib_data(i, 18); pz(i) = ib_data(i, 19)
                    ib_diameter(i) = ib_data(i, 20)*2.0_wp
                end do
 
                write (meshnames(1), '(A,I0,A)') '../p0/', t_step, '.silo:ib_bodies'
                meshtypes(1) = db_pointmesh
                err = dbset2dstrlen(len(meshnames(1)))
                err = dbputmmesh(out%dbroot, 'ib_bodies', 16, 1, meshnames, len_trim(meshnames), meshtypes, db_f77null, ierr)
 
                err = dbputpm(out%dbfile, 'ib_bodies', 9, 3, px, py, pz, nbodies, db_double, db_f77null, ierr)
 
                call s_write_ib_variable('ib_force_x', t_step, force_x, nbodies)
                call s_write_ib_variable('ib_force_y', t_step, force_y, nbodies)
                call s_write_ib_variable('ib_force_z', t_step, force_z, nbodies)
                call s_write_ib_variable('ib_torque_x', t_step, torque_x, nbodies)
                call s_write_ib_variable('ib_torque_y', t_step, torque_y, nbodies)
                call s_write_ib_variable('ib_torque_z', t_step, torque_z, nbodies)
                call s_write_ib_variable('ib_vel_x', t_step, vel_x, nbodies)
                call s_write_ib_variable('ib_vel_y', t_step, vel_y, nbodies)
                call s_write_ib_variable('ib_vel_z', t_step, vel_z, nbodies)
                call s_write_ib_variable('ib_omega_x', t_step, omega_x, nbodies)
                call s_write_ib_variable('ib_omega_y', t_step, omega_y, nbodies)
                call s_write_ib_variable('ib_omega_z', t_step, omega_z, nbodies)
                call s_write_ib_variable('ib_angle_x', t_step, angle_x, nbodies)
                call s_write_ib_variable('ib_angle_y', t_step, angle_y, nbodies)
                call s_write_ib_variable('ib_angle_z', t_step, angle_z, nbodies)
                call s_write_ib_variable('ib_diameter', t_step, ib_diameter, nbodies)
 
                deallocate (ib_data, px, py, pz, force_x, force_y, force_z)
                deallocate (torque_x, torque_y, torque_z, vel_x, vel_y, vel_z)
                deallocate (omega_x, omega_y, omega_z, angle_x, angle_y, angle_z)
                deallocate (ib_diameter)
            end if
        end if
#endif
 
    end subroutine s_write_ib_bodies_to_formatted_database_file
 
    !> Write a single IB point-variable to the Silo database slave and master files.
    subroutine s_write_ib_variable(varname, t_step, data, nBodies)
 
        character(len=*), intent(in)       :: varname
        integer, intent(in)                :: t_step
        real(wp), dimension(:), intent(in) :: data
        integer, intent(in)                :: nBodies
        character(len=4*name_len)          :: var_name_entry
        integer                            :: var_type_entry, ierr
 
        write (var_name_entry, '(A,I0,A)') '../p0/', t_step, '.silo:' // trim(varname)
        var_type_entry = db_pointvar
        err = dbset2dstrlen(len(var_name_entry))
        err = dbputmvar(out%dbroot, trim(varname), len_trim(varname), 1, var_name_entry, len_trim(var_name_entry), &
                        & var_type_entry, db_f77null, ierr)
 
        err = dbputpv1(out%dbfile, trim(varname), len_trim(varname), 'ib_bodies', 9, data, nbodies, db_double, db_f77null, ierr)
 
    end subroutine s_write_ib_variable
 
    !> Close the formatted database slave file and, for the root process, the master file.
    impure subroutine s_close_formatted_database_file()
 
        integer :: ierr
 
        if (format == format_silo) then
            ierr = dbclose(out%dbfile)
            if (proc_rank == 0) ierr = dbclose(out%dbroot)
        else
            close (out%dbfile)
            if (n == 0 .and. proc_rank == 0) close (out%dbroot)
        end if
 
    end subroutine s_close_formatted_database_file
 
    !> Close the interface data file.
    impure subroutine s_close_intf_data_file()
 
        close (211)
 
    end subroutine s_close_intf_data_file
 
    !> Close the energy data file.
    impure subroutine s_close_energy_data_file()
 
        close (251)
 
    end subroutine s_close_energy_data_file
 
    !> Deallocate module arrays and release all data-output resources.
    impure subroutine s_finalize_data_output_module()
 
        deallocate (out%q_sf)
        if (n == 0) deallocate (out%q_root_sf)
        if (grid_geometry == 3) then
            deallocate (out%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 == format_silo) then
            deallocate (out%spatial_extents)
            deallocate (out%data_extents)
            deallocate (out%lo_offset)
            deallocate (out%hi_offset)
            deallocate (out%dims)
        end if
 
    end subroutine s_finalize_data_output_module
 
end module m_data_output