m_data_output.fpp.f90

Go to the documentation of this file.

# 1 "/home/runner/work/MFC/MFC/src/pre_process/m_data_output.fpp"
!>
!! @file
!! @brief Contains module m_data_output
 
!> @brief Writes grid and initial condition data to serial or parallel output files
module m_data_output
 
    use m_derived_types
    use m_global_parameters
    use m_helper
    use m_mpi_proxy
#ifdef MFC_MPI
    use mpi
#endif
 
    use m_compile_specific
    use m_variables_conversion
    use m_helper
    use m_delay_file_access
    use m_boundary_common
    use m_boundary_conditions
    use m_boundary_io
    use m_thermochem, only: species_names
    use m_helper
    use m_constants, only: model_eqns_5eq, precision_single
 
    implicit none
 
    private
    public :: s_write_serial_data_files, s_write_parallel_data_files, s_write_data_files, s_initialize_data_output_module, &
        & s_finalize_data_output_module, s_write_ib_state_0_file
 
    type(scalar_field), allocatable, dimension(:) :: q_cons_temp
 
    abstract interface
 
        !> Interface for the conservative data
        impure subroutine s_write_abstract_data_files(q_cons_vf, q_prim_vf, bc_type, q_T_sf)
 
            import :: scalar_field, integer_field, sys_size, m, n, p, pres_field, num_dims
 
            type(scalar_field), dimension(sys_size), intent(inout)      :: q_cons_vf, q_prim_vf
            type(integer_field), dimension(1:num_dims,-1:1), intent(in) :: bc_type
            type(scalar_field), intent(inout), optional                 :: q_t_sf
 
        end subroutine s_write_abstract_data_files
    end interface
 
    !> Time-step folder into which grid and initial condition data will be placed
    character(LEN=path_len + 2*name_len), private :: t_step_dir
    character(LEN=path_len + 2*name_len), public  :: restart_dir  !< Restart data folder
    procedure(s_write_abstract_data_files), pointer :: s_write_data_files => null()
 
contains
 
    !> Writes grid and initial condition data files to the "0" time-step directory in the local processor rank folder
    impure subroutine s_write_serial_data_files(q_cons_vf, q_prim_vf, bc_type, q_T_sf)
 
        type(scalar_field), dimension(sys_size), intent(inout)      :: q_cons_vf, q_prim_vf
        type(integer_field), dimension(1:num_dims,-1:1), intent(in) :: bc_type
        logical                                                     :: file_exist
        character(LEN=15)                                           :: fmt
        character(LEN=3)                                            :: status
        character(len=int(floor(log10(real(sys_size, wp)))) + 1)    :: file_num
        character(LEN=len_trim(t_step_dir) + name_len)              :: file_loc
        integer                                                     :: i, j, k, l, r, c
        integer                                                     :: t_step
        real(wp), dimension(nb)                                     :: nrtmp
        real(wp)                                                    :: nbub
        real(wp)                                                    :: gamma, pi_inf, qv
        real(wp)                                                    :: rho
        real(wp)                                                    :: pres, t
        real(wp)                                                    :: rhoyks(1:num_species)
        real(wp)                                                    :: pres_mag
        type(scalar_field), intent(inout), optional                 :: q_t_sf
 
        pres_mag = 0._wp
 
        t = dflt_t_guess
 
        t_step = 0
 
        if (old_grid) then
            status = 'old'
        else
            status = 'new'
        end if
 
        if (bc_io) then
            if (igr) then
                call s_write_serial_boundary_condition_files(q_cons_vf, bc_type, t_step_dir, old_grid)
            else
                call s_write_serial_boundary_condition_files(q_prim_vf, bc_type, t_step_dir, old_grid, q_t_sf)
            end if
        end if
 
        file_loc = trim(t_step_dir) // '/x_cb.dat'
        open (1, file=trim(file_loc), form='unformatted', status=status)
        write (1) x_cb(-1:m)
        close (1)
 
        if (n > 0) then
            file_loc = trim(t_step_dir) // '/y_cb.dat'
            open (1, file=trim(file_loc), form='unformatted', status=status)
            write (1) y_cb(-1:n)
            close (1)
 
            if (p > 0) then
                file_loc = trim(t_step_dir) // '/z_cb.dat'
                open (1, file=trim(file_loc), form='unformatted', status=status)
                write (1) z_cb(-1:p)
                close (1)
            end if
        end if
 
        do i = 1, sys_size
            write (file_num, '(I0)') i
            file_loc = trim(t_step_dir) // '/q_cons_vf' // trim(file_num) // '.dat'
            open (1, file=trim(file_loc), form='unformatted', status=status)
            write (1) q_cons_vf(i)%sf(0:m,0:n,0:p)
            close (1)
        end do
 
        if (bubbles_lagrange) then
            block
                real(stp), allocatable                         :: beta_ones(:,:,:)
                character(LEN=len_trim(t_step_dir) + name_len) :: beta_file_loc
                integer                                        :: jj, kk, ll
                allocate (beta_ones(0:m,0:n,0:p))
                do ll = 0, p
                    do kk = 0, n
                        do jj = 0, m
                            beta_ones(jj, kk, ll) = 1.0_stp
                        end do
                    end do
                end do
                write (beta_file_loc, '(A,I0,A)') trim(t_step_dir) // '/q_cons_vf', sys_size + 1, '.dat'
                open (1, file=trim(beta_file_loc), form='unformatted', status=status)
                write (1) beta_ones
                close (1)
                deallocate (beta_ones)
            end block
        end if
 
        if (qbmm .and. .not. polytropic) then
            do i = 1, nb
                do r = 1, nnode
                    write (file_num, '(I0)') r + (i - 1)*nnode + sys_size
                    file_loc = trim(t_step_dir) // '/pb' // trim(file_num) // '.dat'
                    open (1, file=trim(file_loc), form='unformatted', status=status)
                    write (1) pb%sf(:,:,:,r, i)
                    close (1)
                end do
            end do
 
            do i = 1, nb
                do r = 1, nnode
                    write (file_num, '(I0)') r + (i - 1)*nnode + sys_size
                    file_loc = trim(t_step_dir) // '/mv' // trim(file_num) // '.dat'
                    open (1, file=trim(file_loc), form='unformatted', status=status)
                    write (1) mv%sf(:,:,:,r, i)
                    close (1)
                end do
            end do
        end if
 
        if (precision == precision_single) then
            fmt = "(2F30.3)"
        else
            fmt = "(2F40.14)"
        end if
 
        write (t_step_dir, '(A,I0,A,I0)') trim(case_dir) // '/D'
        file_loc = trim(t_step_dir) // '/.'
 
        inquire (file=trim(file_loc), exist=file_exist)
 
        if (.not. file_exist) call s_create_directory(trim(t_step_dir))
 
        if (cfl_dt) t_step = n_start
 
        if (n == 0 .and. p == 0) then
            if (model_eqns == model_eqns_5eq) then
                do i = 1, sys_size
                    write (file_loc, '(A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/prim.', i, '.', proc_rank, '.', t_step, '.dat'
 
                    open (2, file=trim(file_loc))
                    do j = 0, m
                        if (chemistry) then
                            do c = 1, num_species
                                rhoyks(c) = q_cons_vf(eqn_idx%species%beg + c - 1)%sf(j, 0, 0)
                            end do
                        end if
 
                        call s_convert_to_mixture_variables(q_cons_vf, j, 0, 0, rho, gamma, pi_inf, qv)
 
                        if ((i >= eqn_idx%species%beg) .and. (i <= eqn_idx%species%end)) then
                            write (2, fmt) x_cb(j), q_cons_vf(i)%sf(j, 0, 0)/rho
                        else if (((i >= eqn_idx%cont%beg) .and. (i <= eqn_idx%cont%end)) .or. ((i >= eqn_idx%adv%beg) &
                                 & .and. (i <= eqn_idx%adv%end)) .or. ((i >= eqn_idx%species%beg) .and. (i <= eqn_idx%species%end) &
                                 & )) then
                            write (2, fmt) x_cb(j), q_cons_vf(i)%sf(j, 0, 0)
                        else if (i == eqn_idx%mom%beg) then  ! u
                            write (2, fmt) x_cb(j), q_cons_vf(eqn_idx%mom%beg)%sf(j, 0, 0)/rho
                        else if (i == eqn_idx%stress%beg) then  ! tau_e
                            write (2, fmt) x_cb(j), q_cons_vf(eqn_idx%stress%beg)%sf(j, 0, 0)/rho
                        else if (i == eqn_idx%E) then  ! p
                            if (mhd) then
                                pres_mag = 0.5_wp*(bx0**2 + q_cons_vf(eqn_idx%B%beg)%sf(j, 0, &
                                                   & 0)**2 + q_cons_vf(eqn_idx%B%beg + 1)%sf(j, 0, 0)**2)
                            end if
 
                            call s_compute_pressure(q_cons_vf(eqn_idx%E)%sf(j, 0, 0), q_cons_vf(eqn_idx%alf)%sf(j, 0, 0), &
                                                    & 0.5_wp*(q_cons_vf(eqn_idx%mom%beg)%sf(j, 0, 0)**2._wp)/rho, pi_inf, gamma, &
                                                    & rho, qv, rhoyks, pres, t, pres_mag=pres_mag)
                            write (2, fmt) x_cb(j), pres
                        else if (mhd) then
                            if (i == eqn_idx%mom%beg + 1) then  ! v
                                write (2, fmt) x_cb(j), q_cons_vf(eqn_idx%mom%beg + 1)%sf(j, 0, 0)/rho
                            else if (i == eqn_idx%mom%beg + 2) then  ! w
                                write (2, fmt) x_cb(j), q_cons_vf(eqn_idx%mom%beg + 2)%sf(j, 0, 0)/rho
                            else if (i == eqn_idx%B%beg) then  ! By
                                write (2, fmt) x_cb(j), q_cons_vf(eqn_idx%B%beg)%sf(j, 0, 0)/rho
                            else if (i == eqn_idx%B%beg + 1) then  ! Bz
                                write (2, fmt) x_cb(j), q_cons_vf(eqn_idx%B%beg + 1)%sf(j, 0, 0)/rho
                            end if
                        else if ((i >= eqn_idx%bub%beg) .and. (i <= eqn_idx%bub%end) .and. bubbles_euler) then
                            if (qbmm) then
                                nbub = q_cons_vf(eqn_idx%bub%beg)%sf(j, 0, 0)
                            else
                                if (adv_n) then
                                    nbub = q_cons_vf(eqn_idx%n)%sf(j, 0, 0)
                                else
                                    do k = 1, nb
                                        nrtmp(k) = q_cons_vf(qbmm_idx%rs(k))%sf(j, 0, 0)
                                    end do
 
                                    call s_comp_n_from_cons(real(q_cons_vf(eqn_idx%alf)%sf(j, 0, 0), kind=wp), nrtmp, nbub, weight)
                                end if
                            end if
                            write (2, fmt) x_cb(j), q_cons_vf(i)%sf(j, 0, 0)/nbub
                        else if (i == eqn_idx%n .and. adv_n .and. bubbles_euler) then
                            write (2, fmt) x_cb(j), q_cons_vf(i)%sf(j, 0, 0)
                        else if (i == eqn_idx%damage) then
                            write (2, fmt) x_cb(j), q_cons_vf(i)%sf(j, 0, 0)
                        end if
                    end do
                    close (2)
                end do
            end if
 
            do i = 1, sys_size
                write (file_loc, '(A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/cons.', i, '.', proc_rank, '.', t_step, '.dat'
 
                open (2, file=trim(file_loc))
                do j = 0, m
                    write (2, fmt) x_cb(j), q_cons_vf(i)%sf(j, 0, 0)
                end do
                close (2)
            end do
 
            if (qbmm .and. .not. polytropic) then
                do i = 1, nb
                    do r = 1, nnode
                        write (file_loc, '(A,I0,A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/pres.', i, '.', r, '.', proc_rank, &
                               & '.', t_step, '.dat'
 
                        open (2, file=trim(file_loc))
                        do j = 0, m
                            write (2, fmt) x_cb(j), pb%sf(j, 0, 0, r, i)
                        end do
                        close (2)
                    end do
                end do
                do i = 1, nb
                    do r = 1, nnode
                        write (file_loc, '(A,I0,A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/mv.', i, '.', r, '.', proc_rank, &
                               & '.', t_step, '.dat'
 
                        open (2, file=trim(file_loc))
                        do j = 0, m
                            write (2, fmt) x_cb(j), mv%sf(j, 0, 0, r, i)
                        end do
                        close (2)
                    end do
                end do
            end if
        end if
 
        if (precision == precision_single) then
            fmt = "(3F30.7)"
        else
            fmt = "(3F40.14)"
        end if
 
        if ((n > 0) .and. (p == 0)) then
            do i = 1, sys_size
                write (file_loc, '(A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/cons.', i, '.', proc_rank, '.', t_step, '.dat'
                open (2, file=trim(file_loc))
                do j = 0, m
                    do k = 0, n
                        write (2, fmt) x_cb(j), y_cb(k), q_cons_vf(i)%sf(j, k, 0)
                    end do
                    write (2, *)
                end do
                close (2)
            end do
 
            if (qbmm .and. .not. polytropic) then
                do i = 1, nb
                    do r = 1, nnode
                        write (file_loc, '(A,I0,A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/pres.', i, '.', r, '.', proc_rank, &
                               & '.', t_step, '.dat'
 
                        open (2, file=trim(file_loc))
                        do j = 0, m
                            do k = 0, n
                                write (2, fmt) x_cb(j), y_cb(k), pb%sf(j, k, 0, r, i)
                            end do
                        end do
                        close (2)
                    end do
                end do
                do i = 1, nb
                    do r = 1, nnode
                        write (file_loc, '(A,I0,A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/mv.', i, '.', r, '.', proc_rank, &
                               & '.', t_step, '.dat'
 
                        open (2, file=trim(file_loc))
                        do j = 0, m
                            do k = 0, n
                                write (2, fmt) x_cb(j), y_cb(k), mv%sf(j, k, 0, r, i)
                            end do
                        end do
                        close (2)
                    end do
                end do
            end if
        end if
 
        if (precision == precision_single) then
            fmt = "(4F30.7)"
        else
            fmt = "(4F40.14)"
        end if
 
        if (p > 0) then
            do i = 1, sys_size
                write (file_loc, '(A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/cons.', i, '.', proc_rank, '.', t_step, '.dat'
                open (2, file=trim(file_loc))
                do j = 0, m
                    do k = 0, n
                        do l = 0, p
                            write (2, fmt) x_cb(j), y_cb(k), z_cb(l), q_cons_vf(i)%sf(j, k, l)
                        end do
                        write (2, *)
                    end do
                    write (2, *)
                end do
                close (2)
            end do
 
            if (qbmm .and. .not. polytropic) then
                do i = 1, nb
                    do r = 1, nnode
                        write (file_loc, '(A,I0,A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/pres.', i, '.', r, '.', proc_rank, &
                               & '.', t_step, '.dat'
 
                        open (2, file=trim(file_loc))
                        do j = 0, m
                            do k = 0, n
                                do l = 0, p
                                    write (2, fmt) x_cb(j), y_cb(k), z_cb(l), pb%sf(j, k, l, r, i)
                                end do
                            end do
                        end do
                        close (2)
                    end do
                end do
                do i = 1, nb
                    do r = 1, nnode
                        write (file_loc, '(A,I0,A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/mv.', i, '.', r, '.', proc_rank, &
                               & '.', t_step, '.dat'
 
                        open (2, file=trim(file_loc))
                        do j = 0, m
                            do k = 0, n
                                do l = 0, p
                                    write (2, fmt) x_cb(j), y_cb(k), z_cb(l), mv%sf(j, k, l, r, i)
                                end do
                            end do
                        end do
                        close (2)
                    end do
                end do
            end if
        end if
 
    end subroutine s_write_serial_data_files
 
    !> Writes grid and initial condition data files in parallel to the "0" time-step directory in the local processor rank folder
    impure subroutine s_write_parallel_data_files(q_cons_vf, q_prim_vf, bc_type, q_T_sf)
 
        type(scalar_field), dimension(sys_size), intent(inout)      :: q_cons_vf, q_prim_vf
        type(integer_field), dimension(1:num_dims,-1:1), intent(in) :: bc_type
        type(scalar_field), optional, intent(inout)                 :: q_t_sf
 
#ifdef MFC_MPI
        integer                              :: ifile, ierr, data_size
        integer, dimension(MPI_STATUS_SIZE)  :: status
        integer(KIND=MPI_OFFSET_KIND)        :: disp
        integer(KIND=MPI_OFFSET_KIND)        :: m_mok, n_mok, p_mok
        integer(KIND=MPI_OFFSET_KIND)        :: wp_mok, var_mok, str_mok
        integer(KIND=MPI_OFFSET_KIND)        :: nvars_mok
        integer(KIND=MPI_OFFSET_KIND)        :: mok
        character(LEN=path_len + 2*name_len) :: file_loc
        logical                              :: file_exist, dir_check
        integer                              :: i, j, k, l
        real(wp)                             :: loc_violations, glb_violations
        integer                              :: m_ds, n_ds, p_ds
        integer                              :: m_glb_ds, n_glb_ds, p_glb_ds
        integer                              :: m_glb_save, n_glb_save, p_glb_save  !< Size of array being saved
 
        loc_violations = 0._wp
 
        if (down_sample) then
            if ((mod(m + 1, 3) > 0) .or. (mod(n + 1, 3) > 0) .or. (mod(p + 1, 3) > 0)) then
                loc_violations = 1._wp
            end if
            call s_mpi_allreduce_sum(loc_violations, glb_violations)
            if (proc_rank == 0 .and. nint(glb_violations) > 0) then
                print *, &
                    & "WARNING: Attempting to downsample data but there are" &
                    & // "processors with local problem sizes that are not divisible by 3."
            end if
            call s_populate_variables_buffers(bc_type, q_cons_vf)
            call s_downsample_data(q_cons_vf, q_cons_temp, m_ds, n_ds, p_ds, m_glb_ds, n_glb_ds, p_glb_ds)
        end if
 
        if (file_per_process) then
            if (proc_rank == 0) then
                file_loc = trim(case_dir) // '/restart_data/lustre_0'
                call my_inquire(file_loc, dir_check)
                if (dir_check .neqv. .true.) then
                    call s_create_directory(trim(file_loc))
                end if
                call s_create_directory(trim(file_loc))
            end if
            call s_mpi_barrier()
            call s_delay_file_access(proc_rank)
 
            if (down_sample) then
                call s_initialize_mpi_data_ds(m_ds, n_ds, p_ds)
            else
                call s_initialize_mpi_data(q_cons_vf, qbmm_pb=pb, qbmm_mv=mv)
            end if
 
            if (cfl_dt) then
                write (file_loc, '(I0,A,i7.7,A)') n_start, '_', proc_rank, '.dat'
            else
                write (file_loc, '(I0,A,i7.7,A)') t_step_start, '_', proc_rank, '.dat'
            end if
            file_loc = trim(restart_dir) // '/lustre_0' // trim(mpiiofs) // trim(file_loc)
            inquire (file=trim(file_loc), exist=file_exist)
            if (file_exist .and. proc_rank == 0) then
                call mpi_file_delete(file_loc, mpi_info_int, ierr)
            end if
            if (file_exist) call mpi_file_delete(file_loc, mpi_info_int, ierr)
            call mpi_file_open(mpi_comm_self, file_loc, ior(mpi_mode_wronly, mpi_mode_create), mpi_info_int, ifile, ierr)
            call s_check_mpi_file_open(ierr, file_loc)
 
            if (down_sample) then
                data_size = (m_ds + 3)*(n_ds + 3)*(p_ds + 3)
                m_glb_save = m_glb_ds + 3
                n_glb_save = n_glb_ds + 3
                p_glb_save = p_glb_ds + 3
            else
                data_size = (m + 1)*(n + 1)*(p + 1)
                m_glb_save = m_glb + 1
                n_glb_save = n_glb + 1
                p_glb_save = p_glb + 1
            end if
 
            ! Resize some integers so MPI can write even the biggest files
            m_mok = int(m_glb_save, mpi_offset_kind)
            n_mok = int(n_glb_save, mpi_offset_kind)
            p_mok = int(p_glb_save, mpi_offset_kind)
            wp_mok = int(storage_size(0._stp)/8, mpi_offset_kind)
            mok = int(1._wp, mpi_offset_kind)
            str_mok = int(name_len, mpi_offset_kind)
            nvars_mok = int(sys_size, mpi_offset_kind)
 
            if (bubbles_euler) then
                do i = 1, sys_size
                    var_mok = int(i, mpi_offset_kind)
 
                    call mpi_file_write_all(ifile, mpi_io_data%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
                end do
                if (qbmm .and. .not. polytropic) then
                    do i = sys_size + 1, sys_size + 2*nb*nnode
                        var_mok = int(i, mpi_offset_kind)
 
                        call mpi_file_write_all(ifile, mpi_io_data%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
                    end do
                end if
            else
                if (down_sample) then
                    do i = 1, sys_size
                        var_mok = int(i, mpi_offset_kind)
 
                        call mpi_file_write_all(ifile, q_cons_temp(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
                    end do
                else
                    do i = 1, sys_size
                        var_mok = int(i, mpi_offset_kind)
 
                        call mpi_file_write_all(ifile, mpi_io_data%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
                    end do
                end if
            end if
 
            if (bubbles_lagrange) then
                block
                    real(stp), allocatable :: beta_ones(:,:,:)
                    integer                :: jj, kk, ll
                    allocate (beta_ones(0:m,0:n,0:p))
                    do ll = 0, p
                        do kk = 0, n
                            do jj = 0, m
                                beta_ones(jj, kk, ll) = 1.0_stp
                            end do
                        end do
                    end do
                    call mpi_file_write_all(ifile, beta_ones, data_size*mpi_io_type, mpi_io_p, status, ierr)
                    deallocate (beta_ones)
                end block
            end if
 
            call mpi_file_close(ifile, ierr)
        else
            call s_initialize_mpi_data(q_cons_vf, qbmm_pb=pb, qbmm_mv=mv)
 
            if (cfl_dt) then
                write (file_loc, '(I0,A)') n_start, '.dat'
            else
                write (file_loc, '(I0,A)') t_step_start, '.dat'
            end if
            file_loc = trim(restart_dir) // trim(mpiiofs) // trim(file_loc)
            inquire (file=trim(file_loc), exist=file_exist)
            if (file_exist .and. proc_rank == 0) then
                call mpi_file_delete(file_loc, mpi_info_int, ierr)
            end if
            call mpi_file_open(mpi_comm_world, file_loc, ior(mpi_mode_wronly, mpi_mode_create), mpi_info_int, ifile, ierr)
            call s_check_mpi_file_open(ierr, file_loc)
 
            data_size = (m + 1)*(n + 1)*(p + 1)
 
            ! Resize some integers so MPI can write even the biggest files
            m_mok = int(m_glb + 1, mpi_offset_kind)
            n_mok = int(n_glb + 1, mpi_offset_kind)
            p_mok = int(p_glb + 1, mpi_offset_kind)
            wp_mok = int(storage_size(0._stp)/8, mpi_offset_kind)
            mok = int(1._wp, mpi_offset_kind)
            str_mok = int(name_len, mpi_offset_kind)
            nvars_mok = int(sys_size, mpi_offset_kind)
 
            if (bubbles_euler) then
                do i = 1, sys_size
                    var_mok = int(i, mpi_offset_kind)
 
                    disp = m_mok*max(mok, n_mok)*max(mok, p_mok)*wp_mok*(var_mok - 1)
 
                    call mpi_file_set_view(ifile, disp, mpi_io_p, mpi_io_data%view(i), 'native', mpi_info_int, ierr)
                    call mpi_file_write_all(ifile, mpi_io_data%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
                end do
                if (qbmm .and. .not. polytropic) then
                    do i = sys_size + 1, sys_size + 2*nb*nnode
                        var_mok = int(i, mpi_offset_kind)
 
                        disp = m_mok*max(mok, n_mok)*max(mok, p_mok)*wp_mok*(var_mok - 1)
 
                        call mpi_file_set_view(ifile, disp, mpi_io_p, mpi_io_data%view(i), 'native', mpi_info_int, ierr)
                        call mpi_file_write_all(ifile, mpi_io_data%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
                    end do
                end if
            else
                do i = 1, sys_size
                    var_mok = int(i, mpi_offset_kind)
 
                    disp = m_mok*max(mok, n_mok)*max(mok, p_mok)*wp_mok*(var_mok - 1)
 
                    call mpi_file_set_view(ifile, disp, mpi_io_p, mpi_io_data%view(i), 'native', mpi_info_int, ierr)
                    call mpi_file_write_all(ifile, mpi_io_data%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
                end do
            end if
 
            if (bubbles_lagrange) then
                block
                    real(stp), allocatable :: beta_ones(:,:,:)
                    integer                :: jj, kk, ll
                    allocate (beta_ones(0:m,0:n,0:p))
                    do ll = 0, p
                        do kk = 0, n
                            do jj = 0, m
                                beta_ones(jj, kk, ll) = 1.0_stp
                            end do
                        end do
                    end do
                    var_mok = int(sys_size + 1, mpi_offset_kind)
                    disp = m_mok*max(mok, n_mok)*max(mok, p_mok)*wp_mok*(var_mok - 1)
                    call mpi_file_set_view(ifile, disp, mpi_io_p, mpi_io_data%view(1), 'native', mpi_info_int, ierr)
                    call mpi_file_write_all(ifile, beta_ones, data_size*mpi_io_type, mpi_io_p, status, ierr)
                    deallocate (beta_ones)
                end block
            end if
 
            call mpi_file_close(ifile, ierr)
        end if
#endif
 
        if (bc_io) then
            if (igr) then
                call s_write_parallel_boundary_condition_files(q_cons_vf, bc_type)
            else
                call s_write_parallel_boundary_condition_files(q_prim_vf, bc_type, q_t_sf)
            end if
        end if
 
    end subroutine s_write_parallel_data_files
 
    !> Computation of parameters, allocation procedures, and/or any other tasks needed to properly setup the module
    impure subroutine s_initialize_data_output_module
 
        character(LEN=len_trim(case_dir) + 2*name_len) :: file_loc
        character(len=15)                              :: temp
        character(LEN=1), dimension(3), parameter      :: coord = (/'x', 'y', 'z'/)
        logical                                        :: dir_check
        integer                                        :: i, iu
        integer                                        :: m_ds, n_ds, p_ds
 
        if (parallel_io .neqv. .true.) then
            write (t_step_dir, '(A,I0,A)') '/p_all/p', proc_rank, '/0'
            t_step_dir = trim(case_dir) // trim(t_step_dir)
 
            if (old_grid .neqv. .true.) then
                file_loc = trim(t_step_dir) // '/'
 
                call my_inquire(file_loc, dir_check)
 
                if (dir_check) call s_delete_directory(trim(t_step_dir))
 
                call s_create_directory(trim(t_step_dir))
            end if
 
            s_write_data_files => s_write_serial_data_files
        else
            write (restart_dir, '(A)') '/restart_data'
            restart_dir = trim(case_dir) // trim(restart_dir)
 
            if ((old_grid .neqv. .true.) .and. (proc_rank == 0)) then
                file_loc = trim(restart_dir) // '/'
                call my_inquire(file_loc, dir_check)
 
                if (dir_check) call s_delete_directory(trim(restart_dir))
                call s_create_directory(trim(restart_dir))
            end if
 
            call s_mpi_barrier()
 
            s_write_data_files => s_write_parallel_data_files
        end if
 
        open (newunit=iu, file='indices.dat', status='unknown')
 
        write (iu, '(A)') "Warning: The creation of file is currently experimental."
        write (iu, '(A)') "This file may contain errors and not support all features."
 
        write (iu, '(A3,A20,A20)') "#", "Conservative", "Primitive"
        write (iu, '(A)') "    "
        do i = eqn_idx%cont%beg, eqn_idx%cont%end
            write (temp, '(I0)') i - eqn_idx%cont%beg + 1
            write (iu, '(I3,A20,A20)') i, "\alpha_{" // trim(temp) // "} \rho_{" // trim(temp) // "}", &
                   & "\alpha_{" // trim(temp) // "} \rho"
        end do
        do i = eqn_idx%mom%beg, eqn_idx%mom%end
            write (iu, '(I3,A20,A20)') i, "\rho u_" // coord(i - eqn_idx%mom%beg + 1), "u_" // coord(i - eqn_idx%mom%beg + 1)
        end do
        if (eqn_idx%E /= 0) write (iu, '(I3,A20,A20)') eqn_idx%E, "\rho U", "p"
        do i = eqn_idx%adv%beg, eqn_idx%adv%end
            write (temp, '(I0)') i - eqn_idx%cont%beg + 1
            write (iu, '(I3,A20,A20)') i, "\alpha_{" // trim(temp) // "}", "\alpha_{" // trim(temp) // "}"
        end do
        if (chemistry) then
            do i = 1, num_species
                write (iu, '(I3,A20,A20)') eqn_idx%species%beg + i - 1, "Y_{" // trim(species_names(i)) // "} \rho", &
                       & "Y_{" // trim(species_names(i)) // "}"
            end do
        end if
 
        write (iu, '(A)') ""
        call write_range(eqn_idx%cont%beg, eqn_idx%cont%end, " Continuity")
        call write_range(eqn_idx%mom%beg, eqn_idx%mom%end, " Momentum")
        call write_range(eqn_idx%E, eqn_idx%E, " Energy/Pressure")
        call write_range(eqn_idx%adv%beg, eqn_idx%adv%end, " Advection")
        call write_range(eqn_idx%bub%beg, eqn_idx%bub%end, " Bubbles")
        call write_range(eqn_idx%stress%beg, eqn_idx%stress%end, " Stress")
        call write_range(eqn_idx%int_en%beg, eqn_idx%int_en%end, " Internal Energies")
        call write_range(eqn_idx%B%beg, eqn_idx%B%end, " Magnetic Field")
        call write_range(eqn_idx%c, eqn_idx%c, " Color Function")
        call write_range(eqn_idx%species%beg, eqn_idx%species%end, " Chemistry")
 
        close (iu)
 
        if (down_sample) then
            m_ds = int((m + 1)/3) - 1
            n_ds = int((n + 1)/3) - 1
            p_ds = int((p + 1)/3) - 1
 
            allocate (q_cons_temp(1:sys_size))
            do i = 1, sys_size
                allocate (q_cons_temp(i)%sf(-1:m_ds + 1,-1:n_ds + 1,-1:p_ds + 1))
            end do
        end if
 
    contains
 
        subroutine write_range(beg, end, label)
 
            integer, intent(in)      :: beg, end
            character(*), intent(in) :: label
 
            if (beg /= 0) write (iu, '("[",I0,",",I0,"]",A)') beg, end, label
 
        end subroutine write_range
 
    end subroutine s_initialize_data_output_module
 
    !> @brief Writes the initial IB layout (namelist patch_ib entries, then generated particle-cloud beds) that simulation reads
    !! back at startup (s_read_ib_restart_data, src/simulation/m_start_up.fpp), in the layouts simulation's own IB state writers
    !! use. Under file_per_process each rank writes only the IBs it owns to restart_data/lustre_0/ib_state_0_<rank>.dat; otherwise
    !! only rank 0 calls this and writes every IB, in global-id order, to restart_data/ib_state_0.dat.
    impure subroutine s_write_ib_state_0_file(glb_bounds, particle_cloud_ibs, num_particle_cloud_ibs)
 
        type(bounds_info), dimension(3), intent(in)         :: glb_bounds
        type(ib_patch_parameters), dimension(:), intent(in) :: particle_cloud_ibs
        integer, intent(in)                                 :: num_particle_cloud_ibs
        character(LEN=len_trim(case_dir) + 2*name_len)      :: file_loc
        integer                                             :: i, ios, file_unit
        logical, dimension(num_ibs)                         :: owned
        real(wp), dimension(3)                              :: centroid
 
        if (file_per_process) then
            do i = 1, num_ibs
                centroid = [patch_ib(i)%x_centroid, patch_ib(i)%y_centroid, patch_ib(i)%z_centroid]
                owned(i) = f_local_rank_owns_location(centroid, glb_bounds)
            end do
 
            if (proc_rank == 0) call s_create_directory(trim(case_dir) // '/restart_data/lustre_0')
            call s_mpi_barrier()
            call s_delay_file_access(proc_rank)
            write (file_loc, '(A,i7.7,A)') '/restart_data/lustre_0/ib_state_0_', proc_rank, '.dat'
        else
            owned = .true.
            call s_create_directory(trim(case_dir) // '/restart_data')
            file_loc = '/restart_data/ib_state_0.dat'
        end if
        file_loc = trim(case_dir) // trim(file_loc)
 
        open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='replace', iostat=ios)
        if (ios /= 0) call s_mpi_abort('Cannot open IB state output file: ' // trim(file_loc))
 
        if (file_per_process) write (file_unit) count(owned) + num_particle_cloud_ibs
        do i = 1, num_ibs
            if (owned(i)) call s_write_ib_state_record(patch_ib(i), i)
        end do
        do i = 1, num_particle_cloud_ibs
            call s_write_ib_state_record(particle_cloud_ibs(i), particle_cloud_ibs(i)%gbl_patch_id)
        end do
 
        close (file_unit)
 
    contains
 
        !> Writes one 20-field IB state record, prefixed by its global id under file_per_process.
        subroutine s_write_ib_state_record(ib_patch, gbl_id)
 
            type(ib_patch_parameters), intent(in) :: ib_patch
            integer, intent(in)                   :: gbl_id
            real(wp), dimension(20)               :: ib_buf
 
            ib_buf = 0._wp
            ib_buf(8:10) = ib_patch%vel
            ib_buf(11:13) = ib_patch%angular_vel
            ib_buf(14:16) = ib_patch%angles
            ib_buf(17) = ib_patch%x_centroid
            ib_buf(18) = ib_patch%y_centroid
            ib_buf(19) = ib_patch%z_centroid
            ib_buf(20) = ib_patch%radius
 
            if (file_per_process) write (file_unit) gbl_id
            write (file_unit) ib_buf
 
        end subroutine s_write_ib_state_record
 
    end subroutine s_write_ib_state_0_file
 
    !> Resets s_write_data_files pointer
    impure subroutine s_finalize_data_output_module
 
        integer :: i
 
        s_write_data_files => null()
 
        if (down_sample) then
            do i = 1, sys_size
                deallocate (q_cons_temp(i)%sf)
            end do
            deallocate (q_cons_temp)
        end if
 
    end subroutine s_finalize_data_output_module
 
end module m_data_output