m_start_up.fpp.f90

Go to the documentation of this file.

# 1 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
# 1 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 1
# 1 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 1
# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 12 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 1 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp" 1
! Default for Fypp runs outside the toolchain (e.g. documentation). Searched last, so the
! build/include/<target>/thermochem.fpp the toolchain writes from m_thermochem takes precedence.
# 4 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
# 5 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
# 13 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 2
 
# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
! New line at end of file is required for FYPP
# 2 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
# 1 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 1
# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 12 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 1 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp" 1
! Default for Fypp runs outside the toolchain (e.g. documentation). Searched last, so the
! build/include/<target>/thermochem.fpp the toolchain writes from m_thermochem takes precedence.
# 4 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
# 5 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
# 13 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 2
 
# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
! New line at end of file is required for FYPP
# 2 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 2
 
# 4 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 5 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 6 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 7 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 8 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 20 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 43 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 48 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 53 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 58 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 63 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 68 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 81 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 95 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 100 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 105 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 110 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 115 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 120 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 125 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 130 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 135 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 165 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 219 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 233 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 258 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 269 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 271 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 282 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 332 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 342 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 352 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 361 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 378 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 388 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 395 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 401 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 407 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 413 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 419 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 425 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
! New line at end of file is required for FYPP
# 3 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
# 1 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 1
# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 12 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 1 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp" 1
! Default for Fypp runs outside the toolchain (e.g. documentation). Searched last, so the
! build/include/<target>/thermochem.fpp the toolchain writes from m_thermochem takes precedence.
# 4 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
# 5 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
# 13 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 2
 
# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
! New line at end of file is required for FYPP
# 2 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 2
 
# 7 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 17 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 22 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 27 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 32 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 37 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 42 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 47 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 52 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 57 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 62 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 73 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 78 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 83 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 88 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 103 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 131 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 160 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 175 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 193 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 215 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 244 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 259 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 269 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 278 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 294 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 304 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 311 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
! New line at end of file is required for FYPP
# 4 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
 
! GPU parallel region (scalar reductions, maxval/minval)
# 23 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! GPU parallel loop over threads (most common GPU macro)
# 43 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Required closing for GPU_PARALLEL_LOOP
# 55 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Mark routine for device compilation
# 112 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Declare device-resident data
# 130 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Inner loop within a GPU parallel region
# 145 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Scoped GPU data region
# 164 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Host code with device pointers (for MPI with GPU buffers)
# 193 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Allocate device memory (unscoped)
# 207 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Free device memory
# 219 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Atomic operation on device
# 231 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! End atomic capture block
# 242 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Copy data between host and device
# 254 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Synchronization barrier
# 266 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Import GPU library module (openacc or omp_lib)
# 275 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Emit code only for AMD compiler
# 282 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Emit code for non-Cray compilers
# 289 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Emit code only for Cray compiler
# 296 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Emit code for non-NVIDIA compilers
# 303 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 305 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
# 306 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
! New line at end of file is required for FYPP
# 2 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 2
 
# 14 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
! Caution: This macro requires the use of a binding script to set CUDA_VISIBLE_DEVICES, such that we have one GPU device per MPI
! rank. That's because for both cudaMemAdvise (preferred location) and cudaMemPrefetchAsync we use location = device_id = 0. For an
! example see misc/nvidia_uvm/bind.sh.
# 52 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
! Allocate and create GPU device memory
# 72 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
! Free GPU device memory and deallocate
# 80 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
! Cray-specific GPU pointer setup for vector fields
# 104 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
! Cray-specific GPU pointer setup for scalar fields
# 120 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
! Cray-specific GPU pointer setup for acoustic source spatials
# 145 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 151 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 158 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
! New line at end of file is required for FYPP
# 2 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp" 2
 
!>
!! @file
!! @brief  Contains module m_start_up
 
!> @brief Reads and validates user inputs, allocates variables, and configures MPI decomposition and I/O for post-processing
 
module m_start_up
 
    use, intrinsic :: iso_c_binding
 
    use m_derived_types
    use m_global_parameters
    use m_mpi_proxy
    use m_mpi_common
    use m_boundary_common
    use m_boundary_io
    use m_variables_conversion
    use m_eos
    use m_data_input
    use m_data_output
    use m_derived_variables
    use m_helper
    use m_compile_specific
    use m_checker_common
    use m_checker
    use m_thermochem, only: num_species, species_names
    use m_finite_differences
    use m_constants, only: model_eqns_gamma_law, model_eqns_5eq, model_eqns_6eq, format_silo
    use m_chemistry
 
#ifdef MFC_MPI
    use mpi
#endif
 
    implicit none
 
    include 'fftw3.f03'
 
    type(c_ptr)                             :: fwd_plan_x, fwd_plan_y, fwd_plan_z
    complex(c_double_complex), allocatable  :: data_in(:), data_out(:)
    complex(c_double_complex), allocatable  :: data_cmplx(:,:,:), data_cmplx_y(:,:,:), data_cmplx_z(:,:,:)
    real(wp), allocatable, dimension(:,:,:) :: en_real
    real(wp), allocatable, dimension(:)     :: en
    integer                                 :: nx, ny, nz, nxloc, nyloc, nyloc2, nzloc, nf
    integer                                 :: ierr
    integer                                 :: mpi_comm_cart, mpi_comm_cart12, mpi_comm_cart13
    integer, dimension(3)                   :: cart3d_coords
    integer, dimension(2)                   :: cart2d12_coords, cart2d13_coords
    integer                                 :: proc_rank12, proc_rank13
 
contains
 
    !> Reads the configuration file post_process.inp, in order to populate parameters in module m_global_parameters.f90 with the
    !! user provided inputs
    impure subroutine s_read_input_file
 
        character(LEN=name_len) :: file_loc
        logical                 :: file_check
        integer                 :: iostatus
        character(len=1000)     :: line
 
# 1 "/home/runner/work/MFC/MFC/build/include/post_process/generated_namelist.fpp" 1
! AUTO-GENERATED - do not edit directly. Regenerate: cmake reconfigure
!
namelist /user_inputs/ bx0, ca, e_wrt, g, r0ref, re_inv, t_wrt, web, adv_n, alpha_rho_e_wrt, alpha_rho_wrt, alpha_wrt, &
    & alt_soundspeed, avg_state, bc_x, bc_y, bc_z, bub_pp, bubbles_euler, bubbles_lagrange, c_wrt, case_dir, cf_wrt, &
    & cfl_adap_dt, cfl_const_dt, cfl_target, chem_wrt_t, chem_wrt_y, cons_vars_wrt, cont_damage, cyl_coord, down_sample, &
    & fd_order, fft_wrt, file_per_process, fluid_pp, flux_lim, flux_wrt, format, gamma_wrt, heat_ratio_wrt, hyper_cleaning, &
    & hypoelasticity, ib, ib_force_stride, ib_force_wrt, ib_state_wrt, igr, igr_order, lag_betac_wrt, lag_betat_wrt, lag_db_wrt, &
    & lag_dphidt_wrt, lag_header, lag_id_wrt, lag_mg_wrt, lag_mv_wrt, lag_pos_prev_wrt, lag_pos_wrt, lag_pres_wrt, lag_r0_wrt, &
    & lag_rad_wrt, lag_rmax_wrt, lag_rmin_wrt, lag_rvel_wrt, lag_txt_wrt, lag_vel_wrt, liutex_wrt, m, mhd, mixture_err, &
    & model_eqns, mom_wrt, mpp_lim, muscl_order, n, n_start, nb, num_bc_patches, num_fluids, num_ibs, num_particle_clouds, &
    & omega_wrt, output_partial_domain, p, parallel_io, particle_cloud, pi_inf_wrt, poly_sigma, polydisperse, polytropic, &
    & precision, pres_inf_wrt, pres_wrt, prim_vars_wrt, qbmm, qm_wrt, rburn, reactive_burn, recon_type, relativity, relax, &
    & relax_model, rho_wrt, schlieren_alpha, schlieren_wrt, sigr, sigma, sim_data, surface_tension, t_save, t_step_save, &
    & t_step_start, t_step_stop, t_stop, thermal, vel_wrt, weno_order, x_output, y_output, z_output
# 65 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp" 2
 
        file_loc = 'post_process.inp'
        inquire (file=trim(file_loc), exist=file_check)
 
        if (file_check) then
            open (1, file=trim(file_loc), form='formatted', status='old', action='read')
            read (1, nml=user_inputs, iostat=iostatus)
 
            if (iostatus /= 0) then
                backspace(1)
                read (1, fmt='(A)') line
                print *, 'Invalid line in namelist: ' // trim(line)
                call s_mpi_abort('Invalid line in post_process.inp. It is ' // 'likely due to a datatype mismatch. Exiting.')
            end if
 
            close (1)
 
            call s_update_cell_bounds(cells_bounds, m, n, p)
 
            if (down_sample) then
                m = int((m + 1)/3) - 1
                n = int((n + 1)/3) - 1
                p = int((p + 1)/3) - 1
            end if
 
            m_glb = m
            n_glb = n
            p_glb = p
 
            nglobal = int(m_glb + 1, kind=8)*int(n_glb + 1, kind=8)*int(p_glb + 1, kind=8)
 
            if (cfl_adap_dt .or. cfl_const_dt) cfl_dt = .true.
 
            if (any((/bc_x%beg, bc_x%end, bc_y%beg, bc_y%end, bc_z%beg, bc_z%end/) == -17) .or. num_bc_patches > 0) then
                bc_io = .true.
            end if
        else
            call s_mpi_abort('File post_process.inp is missing. Exiting.')
        end if
 
    end subroutine s_read_input_file
 
    !> Checking that the user inputs make sense, i.e. that the individual choices are compatible with the code's options and that
    !! the combination of these choices results into a valid configuration for the post-process
    impure subroutine s_check_input_file
 
        character(LEN=len_trim(case_dir)) :: file_loc
        logical                           :: dir_check
 
        case_dir = adjustl(case_dir)
 
        file_loc = trim(case_dir) // '/.'
 
        call my_inquire(file_loc, dir_check)
 
        if (dir_check .neqv. .true.) then
            call s_mpi_abort('Unsupported choice for the value of ' // 'case_dir. Exiting.')
        end if
 
        call s_check_inputs_common(check_total_cells=.true., n_global=nglobal)
        call s_check_inputs()
 
    end subroutine s_check_input_file
 
    !> Load grid and conservative data for a time step, fill ghost-cell buffers, and convert to primitive variables.
    impure subroutine s_perform_time_step(t_step)
 
        integer, intent(inout) :: t_step
        integer                :: eta_hh, eta_mm, eta_ss
        real(wp)               :: eta_sec
 
        if (proc_rank == 0) then
            if (cfl_dt) then
                eta_sec = wall_time_avg*real(n_save - 1 - t_step, wp)
                eta_hh = int(eta_sec)/3600
                eta_mm = mod(int(eta_sec), 3600)/60
                eta_ss = mod(int(eta_sec), 60)
                print '(" [", I3, "%] Saving ", I0, " of ", I0, " t/step ", ES9.2, "s (avg ", ES9.2, "s) ETA ", I0, ":", I2.2, ":", I2.2)', &
                    & int(ceiling(100._wp*(real(t_step - n_start)/(n_save)))), t_step, n_save, wall_time, wall_time_avg, eta_hh, &
                    & eta_mm, eta_ss
            else
                eta_sec = wall_time_avg*real((t_step_stop - t_step)/t_step_save, wp)
                eta_hh = int(eta_sec)/3600
                eta_mm = mod(int(eta_sec), 3600)/60
                eta_ss = mod(int(eta_sec), 60)
                print '(" [", I3, "%] Saving ", I0, " of ", I0, " (t_step ", I0, ") t/step ", ES9.2, "s (avg ", ES9.2, "s) ETA ", I0, ":", I2.2, ":", I2.2)', &
                    & int(ceiling(100._wp*(real(t_step - t_step_start)/(t_step_stop - t_step_start + 1)))), &
                    & (t_step - t_step_start)/t_step_save + 1, (t_step_stop - t_step_start)/t_step_save + 1, t_step, wall_time, &
                    & wall_time_avg, eta_hh, eta_mm, eta_ss
            end if
        end if
 
        call s_read_data_files(t_step)
 
        ! seed the chemistry temperature over the INTERIOR only (mirrors the simulation,
        ! m_start_up): the ghost q_cons is unread at this point, so a ghost-inclusive sweep
        ! would Newton-iterate on garbage (NaN under NaN-init builds) at rank seams and
        ! physical boundaries; s_populate_variables_buffers below extends q_T into the ghosts
        if (chemistry) call s_compute_q_t_sf(q_t_sf, q_cons_vf, idwint)
 
        if (buff_size > 0) then
            if (n == 0) then
                call s_populate_grid_variables_buffers(x_cb, x_cc, dx, offset_x, offset_y, offset_z)
            else if (p == 0) then
                call s_populate_grid_variables_buffers(x_cb, x_cc, dx, offset_x, offset_y, offset_z, y_cb, y_cc, dy)
            else
                call s_populate_grid_variables_buffers(x_cb, x_cc, dx, offset_x, offset_y, offset_z, y_cb, y_cc, dy, z_cb, z_cc, dz)
            end if
            call s_populate_variables_buffers(bc_type, q_cons_vf, q_t_sf=q_t_sf)
        end if
 
        call s_convert_conservative_to_primitive_variables(q_cons_vf, q_t_sf, q_prim_vf, idwbuff)
 
    end subroutine s_perform_time_step
 
    !> Derive requested flow quantities from primitive variables and write them to the formatted database files.
    impure subroutine s_save_data(t_step, varname, pres, c)
 
        integer, intent(inout)                 :: t_step
        character(LEN=name_len), intent(inout) :: varname
        real(wp), intent(inout)                :: pres, c
 
        real(wp), dimension(-offset_x%beg:m + offset_x%end,-offset_y%beg:n + offset_y%end, &
             & -offset_z%beg:p + offset_z%end) :: liutex_mag
        real(wp), dimension(-offset_x%beg:m + offset_x%end,-offset_y%beg:n + offset_y%end,-offset_z%beg:p + offset_z%end, &
             & 3) :: liutex_axis
        integer                         :: i, j, k, l, kx, ky, kz, kf, j_glb, k_glb, l_glb
        character(50)                   :: filename
        logical                         :: file_exists
        real(wp), dimension(num_fluids) :: alpha_rho
        real(wp)                        :: t
        integer                         :: x_beg, x_end, y_beg, y_end, z_beg, z_end
 
        if (output_partial_domain) then
            call s_define_output_region
            x_beg = -offset_x%beg + x_output_idx%beg
            x_end = offset_x%end + x_output_idx%end
            y_beg = -offset_y%beg + y_output_idx%beg
            y_end = offset_y%end + y_output_idx%end
            z_beg = -offset_z%beg + z_output_idx%beg
            z_end = offset_z%end + z_output_idx%end
        else
            x_beg = -offset_x%beg
            x_end = offset_x%end + m
            y_beg = -offset_y%beg
            y_end = offset_y%end + n
            z_beg = -offset_z%beg
            z_end = offset_z%end + p
        end if
 
        call s_open_formatted_database_file(t_step)
 
        if (sim_data .and. proc_rank == 0) then
            call s_open_intf_data_file()
            call s_open_energy_data_file()
        end if
 
        if (sim_data) then
            call s_write_intf_data_file(q_prim_vf)
            call s_write_energy_data_file(q_prim_vf, q_cons_vf)
        end if
 
        call s_write_grid_to_formatted_database_file(t_step)
 
        if (omega_wrt(2) .or. omega_wrt(3) .or. qm_wrt .or. liutex_wrt .or. schlieren_wrt) then
            call s_compute_finite_difference_coefficients(m, x_cc, fd%fd_coeff_x, buff_size, fd_number, fd_order, offset_x)
        end if
 
        if (omega_wrt(1) .or. omega_wrt(3) .or. qm_wrt .or. liutex_wrt .or. (n > 0 .and. schlieren_wrt)) then
            call s_compute_finite_difference_coefficients(n, y_cc, fd%fd_coeff_y, buff_size, fd_number, fd_order, offset_y)
        end if
 
        if (omega_wrt(1) .or. omega_wrt(2) .or. qm_wrt .or. liutex_wrt .or. (p > 0 .and. schlieren_wrt)) then
            call s_compute_finite_difference_coefficients(p, z_cc, fd%fd_coeff_z, buff_size, fd_number, fd_order, offset_z)
        end if
 
        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
                    write (varname, '(A,I0)') 'alpha_rho', i
                    call s_write_field(varname, t_step, q_cons_vf(i), x_beg, x_end, y_beg, y_end, z_beg, z_end)
                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%q_sf(:,:,:) = rho_sf(x_beg:x_end,y_beg:y_end,z_beg:z_end)
            write (varname, '(A)') 'rho'
            call s_write_field(varname, t_step)
        end if
 
        if (relativity .and. (rho_wrt .or. prim_vars_wrt)) then
            write (varname, '(A)') 'rho'
            call s_write_field(varname, t_step, q_prim_vf(1), x_beg, x_end, y_beg, y_end, z_beg, z_end)
        end if
 
        if (relativity .and. (rho_wrt .or. cons_vars_wrt)) then
            ! For relativistic flow, conservative and primitive densities are different Hard-coded single-component for now
            write (varname, '(A)') 'D'
            call s_write_field(varname, t_step, q_cons_vf(1), x_beg, x_end, y_beg, y_end, z_beg, z_end)
        end if
 
        do i = 1, eqn_idx%E - eqn_idx%mom%beg
            if (mom_wrt(i) .or. cons_vars_wrt) then
                write (varname, '(A,I0)') 'mom', i
                call s_write_field(varname, t_step, q_cons_vf(i + eqn_idx%cont%end), x_beg, x_end, y_beg, y_end, z_beg, z_end)
            end if
        end do
 
        do i = 1, eqn_idx%E - eqn_idx%mom%beg
            if (vel_wrt(i) .or. prim_vars_wrt) then
                write (varname, '(A,I0)') 'vel', i
                call s_write_field(varname, t_step, q_prim_vf(i + eqn_idx%cont%end), x_beg, x_end, y_beg, y_end, z_beg, z_end)
            end if
        end do
 
        if (chemistry) then
            do i = 1, num_species
                if (chem_wrt_y(i) .or. prim_vars_wrt) then
                    write (varname, '(A,A)') 'Y_', trim(species_names(i))
                    call s_write_field(varname, t_step, q_prim_vf(eqn_idx%species%beg + i - 1), x_beg, x_end, y_beg, y_end, &
                                       & z_beg, z_end)
                end if
            end do
 
            if (chem_wrt_t) then
                out%q_sf(:,:,:) = q_t_sf%sf(x_beg:x_end,y_beg:y_end,z_beg:z_end)
                write (varname, '(A)') 'T'
                call s_write_field(varname, t_step)
            end if
        end if
 
        do i = 1, eqn_idx%E - eqn_idx%mom%beg
            if (flux_wrt(i)) then
                call s_derive_flux_limiter(i, q_prim_vf, out%q_sf)
                write (varname, '(A,I0)') 'flux', i
                call s_write_field(varname, t_step)
            end if
        end do
 
        if (e_wrt .or. cons_vars_wrt) then
            write (varname, '(A)') 'E'
            call s_write_field(varname, t_step, q_cons_vf(eqn_idx%E), x_beg, x_end, y_beg, y_end, z_beg, z_end)
        end if
 
        if (model_eqns == model_eqns_6eq) then
            do i = 1, num_fluids
                if (alpha_rho_e_wrt(i) .or. cons_vars_wrt) then
                    write (varname, '(A,I0)') 'alpha_rho_e', i
                    call s_write_field(varname, t_step, q_cons_vf(i + eqn_idx%int_en%beg - 1), x_beg, x_end, y_beg, y_end, z_beg, &
                                       & z_end)
                end if
            end do
        end if
 
        if (fft_wrt) then
            do l = 0, p
                do k = 0, n
                    do j = 0, m
                        data_cmplx(j + 1, k + 1, l + 1) = cmplx(q_cons_vf(eqn_idx%mom%beg)%sf(j, k, l)/q_cons_vf(1)%sf(j, k, l), &
                                   & 0._wp)
                    end do
                end do
            end do
 
            call s_mpi_fft_fwd()
 
            en_real = 0.5_wp*abs(data_cmplx_z)**2._wp/(1._wp*nx*ny*nz)**2._wp
 
            do l = 0, p
                do k = 0, n
                    do j = 0, m
                        data_cmplx(j + 1, k + 1, l + 1) = cmplx(q_cons_vf(eqn_idx%mom%beg + 1)%sf(j, k, l)/q_cons_vf(1)%sf(j, k, &
                                   & l), 0._wp)
                    end do
                end do
            end do
 
            call s_mpi_fft_fwd()
 
            en_real = en_real + 0.5_wp*abs(data_cmplx_z)**2._wp/(1._wp*nx*ny*nz)**2._wp
 
            do l = 0, p
                do k = 0, n
                    do j = 0, m
                        data_cmplx(j + 1, k + 1, l + 1) = cmplx(q_cons_vf(eqn_idx%mom%beg + 2)%sf(j, k, l)/q_cons_vf(1)%sf(j, k, &
                                   & l), 0._wp)
                    end do
                end do
            end do
 
            call s_mpi_fft_fwd()
 
            en_real = en_real + 0.5_wp*abs(data_cmplx_z)**2._wp/(1._wp*nx*ny*nz)**2._wp
 
            do kf = 1, nf
                en(kf) = 0._wp
            end do
 
            do l = 1, nz
                do k = 1, nyloc2
                    do j = 1, nxloc
                        j_glb = j + cart3d_coords(2)*nxloc
                        k_glb = k + cart3d_coords(3)*nyloc2
                        l_glb = l
 
                        if (j_glb >= (m_glb + 1)/2) then
                            kx = (j_glb - 1) - (m_glb + 1)
                        else
                            kx = j_glb - 1
                        end if
 
                        if (k_glb >= (n_glb + 1)/2) then
                            ky = (k_glb - 1) - (n_glb + 1)
                        else
                            ky = k_glb - 1
                        end if
 
                        if (l_glb >= (p_glb + 1)/2) then
                            kz = (l_glb - 1) - (p_glb + 1)
                        else
                            kz = l_glb - 1
                        end if
 
                        kf = nint(sqrt(kx**2._wp + ky**2._wp + kz**2._wp)) + 1
 
                        en(kf) = en(kf) + en_real(j, k, l)
                    end do
                end do
            end do
 
#ifdef MFC_MPI
            call mpi_allreduce(mpi_in_place, en, nf, mpi_p, mpi_sum, mpi_comm_world, ierr)
#endif
 
            if (proc_rank == 0) then
                call s_create_directory('En_FFT_DATA')
                write (filename, '(a,i0,a)') 'En_FFT_DATA/En_tot', t_step, '.dat'
                inquire (file=filename, exist=file_exists)
                if (file_exists) then
                    call s_delete_file(trim(filename))
                end if
            end if
 
            do kf = 1, nf
                if (proc_rank == 0) then
                    write (filename, '(a,i0,a)') 'En_FFT_DATA/En_tot', t_step, '.dat'
                    inquire (file=filename, exist=file_exists)
                    if (file_exists) then
                        open (1, file=filename, position='append', status='old')
                        write (1, *) en(kf), t_step
                        close (1)
                    else
                        open (1, file=filename, status='new')
                        write (1, *) en(kf), t_step
                        close (1)
                    end if
                end if
            end do
        end if
 
        if (mhd .and. prim_vars_wrt) then
            do i = eqn_idx%B%beg, eqn_idx%B%end
                ! 1D: output By, Bz
                if (n == 0) then
                    if (i == eqn_idx%B%beg) then
                        write (varname, '(A)') 'By'
                    else
                        write (varname, '(A)') 'Bz'
                    end if
                    ! 2D/3D: output Bx, By, Bz
                else
                    if (i == eqn_idx%B%beg) then
                        write (varname, '(A)') 'Bx'
                    else if (i == eqn_idx%B%beg + 1) then
                        write (varname, '(A)') 'By'
                    else
                        write (varname, '(A)') 'Bz'
                    end if
                end if
                call s_write_field(varname, t_step, q_prim_vf(i), x_beg, x_end, y_beg, y_end, z_beg, z_end)
            end do
        end if
 
        if (hypoelasticity) then
            do i = 1, eqn_idx%stress%end - eqn_idx%stress%beg + 1
                if (prim_vars_wrt) then
                    write (varname, '(A,I0)') 'tau', i
                    call s_write_field(varname, t_step, q_prim_vf(i - 1 + eqn_idx%stress%beg), x_beg, x_end, y_beg, y_end, z_beg, &
                                       & z_end)
                end if
            end do
        end if
 
        if (cont_damage) then
            write (varname, '(A)') 'damage_state'
            call s_write_field(varname, t_step, q_prim_vf(eqn_idx%damage), x_beg, x_end, y_beg, y_end, z_beg, z_end)
        end if
 
        if (hyper_cleaning) then
            write (varname, '(A)') 'psi'
            call s_write_field(varname, t_step, q_cons_vf(eqn_idx%psi), x_beg, x_end, y_beg, y_end, z_beg, z_end)
        end if
 
        if (pres_wrt .or. prim_vars_wrt) then
            write (varname, '(A)') 'pres'
            call s_write_field(varname, t_step, q_prim_vf(eqn_idx%E), x_beg, x_end, y_beg, y_end, z_beg, z_end)
        end if
 
        if (((model_eqns == model_eqns_5eq) .and. (bubbles_euler .neqv. .true.)) .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
                    write (varname, '(A,I0)') 'alpha', i
                    call s_write_field(varname, t_step, q_cons_vf(i + eqn_idx%E), x_beg, x_end, y_beg, y_end, z_beg, z_end)
                end if
            end do
 
            if (alpha_wrt(num_fluids) .or. (cons_vars_wrt .or. prim_vars_wrt)) then
                if (igr) then
                    do k = z_beg, z_end
                        do j = y_beg, y_end
                            do i = x_beg, x_end
                                out%q_sf(i, j, k) = 1._wp
                                do l = 1, num_fluids - 1
                                    out%q_sf(i, j, k) = out%q_sf(i, j, k) - q_cons_vf(eqn_idx%E + l)%sf(i, j, k)
                                end do
                            end do
                        end do
                    end do
                else
                    out%q_sf(:,:,:) = q_cons_vf(eqn_idx%adv%end)%sf(x_beg:x_end,y_beg:y_end,z_beg:z_end)
                end if
                write (varname, '(A,I0)') 'alpha', num_fluids
                call s_write_field(varname, t_step)
            end if
        end if
 
        if (gamma_wrt .or. (model_eqns == model_eqns_gamma_law .and. (cons_vars_wrt .or. prim_vars_wrt))) then
            out%q_sf(:,:,:) = gamma_sf(x_beg:x_end,y_beg:y_end,z_beg:z_end)
            write (varname, '(A)') 'gamma'
            call s_write_field(varname, t_step)
        end if
 
        if (heat_ratio_wrt) then
            call s_derive_specific_heat_ratio(out%q_sf)
            write (varname, '(A)') 'heat_ratio'
            call s_write_field(varname, t_step)
        end if
 
        if (pi_inf_wrt .or. (model_eqns == model_eqns_gamma_law .and. (cons_vars_wrt .or. prim_vars_wrt))) then
            out%q_sf(:,:,:) = pi_inf_sf(x_beg:x_end,y_beg:y_end,z_beg:z_end)
            write (varname, '(A)') 'pi_inf'
            call s_write_field(varname, t_step)
        end if
 
        if (pres_inf_wrt) then
            call s_derive_liquid_stiffness(out%q_sf)
            write (varname, '(A)') 'pres_inf'
            call s_write_field(varname, t_step)
        end if
 
        if (c_wrt) then
            do k = -offset_z%beg, p + offset_z%end
                do j = -offset_y%beg, n + offset_y%end
                    do i = -offset_x%beg, m + offset_x%end
                        do l = 1, eqn_idx%adv%end - eqn_idx%E
                            adv(l) = q_prim_vf(eqn_idx%E + l)%sf(i, j, k)
                            alpha_rho(l) = q_prim_vf(eqn_idx%cont%beg + l - 1)%sf(i, j, k)
                        end do
 
                        pres = q_prim_vf(eqn_idx%E)%sf(i, j, k)
 
                        call s_compute_speed_of_sound(pres, rho_sf(i, j, k), gamma_sf(i, j, k), pi_inf_sf(i, j, k), adv, c, &
                                                      & alpha_rho)
 
                        out%q_sf(i, j, k) = c
                    end do
                end do
            end do
 
            write (varname, '(A)') 'c'
            call s_write_field(varname, t_step)
        end if
 
        if (t_wrt) then
            do l = 1, num_fluids
                do k = -offset_z%beg, p + offset_z%end
                    do j = -offset_y%beg, n + offset_y%end
                        do i = -offset_x%beg, m + offset_x%end
                            call s_phase_temperature(q_prim_vf(eqn_idx%cont%beg + l - 1)%sf(i, j, &
                                                     & k)/max(q_prim_vf(eqn_idx%E + l)%sf(i, j, k), sgm_eps), &
                                                     & q_prim_vf(eqn_idx%E)%sf(i, j, k), l, t)
                            out%q_sf(i, j, k) = t
                        end do
                    end do
                end do
                write (varname, '(A,I0)') 'T', l
                call s_write_field(varname, t_step)
            end do
        end if
 
        do i = 1, 3
            if (omega_wrt(i)) then
                call s_derive_vorticity_component(i, q_prim_vf, out%q_sf)
                write (varname, '(A,I0)') 'omega', i
                call s_write_field(varname, t_step)
            end if
        end do
 
        if (ib) then
            out%q_sf(:,:,:) = real(ib_markers%sf(-offset_x%beg:m + offset_x%end,-offset_y%beg:n + offset_y%end, &
                     & -offset_z%beg:p + offset_z%end), wp)
            varname = 'ib_markers'
            call s_write_field(varname, t_step)
        end if
 
        if (p > 0 .and. qm_wrt) then
            call s_derive_qm(q_prim_vf, out%q_sf)
            write (varname, '(A)') 'qm'
            call s_write_field(varname, t_step)
        end if
 
        if (liutex_wrt) then
            call s_derive_liutex(q_prim_vf, liutex_mag, liutex_axis)
 
            out%q_sf = liutex_mag
            write (varname, '(A)') 'liutex_mag'
            call s_write_field(varname, t_step)
 
            do i = 1, 3
                out%q_sf = liutex_axis(:,:,:,i)
                write (varname, '(A,I0)') 'liutex_axis', i
                call s_write_field(varname, t_step)
            end do
        end if
 
        if (schlieren_wrt) then
            call s_derive_numerical_schlieren_function(q_cons_vf, out%q_sf)
            write (varname, '(A)') 'schlieren'
            call s_write_field(varname, t_step)
        end if
 
        if (cf_wrt) then
            write (varname, '(A,I0)') 'color_function'
            call s_write_field(varname, t_step, q_cons_vf(eqn_idx%c), x_beg, x_end, y_beg, y_end, z_beg, z_end)
        end if
 
        if (bubbles_euler) then
            do i = eqn_idx%adv%beg, eqn_idx%adv%end
                write (varname, '(A,I0)') 'alpha', i - eqn_idx%E
                call s_write_field(varname, t_step, q_cons_vf(i), x_beg, x_end, y_beg, y_end, z_beg, z_end)
            end do
        end if
 
        if (bubbles_euler) then
            ! nR
            do i = 1, nb
                write (varname, '(A,I3.3)') 'nR', i
                call s_write_field(varname, t_step, q_cons_vf(qbmm_idx%rs(i)), x_beg, x_end, y_beg, y_end, z_beg, z_end)
            end do
 
            ! nRdot
            do i = 1, nb
                write (varname, '(A,I3.3)') 'nV', i
                call s_write_field(varname, t_step, q_cons_vf(qbmm_idx%vs(i)), x_beg, x_end, y_beg, y_end, z_beg, z_end)
            end do
            if ((polytropic .neqv. .true.) .and. (.not. qbmm)) then
                ! nP
                do i = 1, nb
                    write (varname, '(A,I3.3)') 'nP', i
                    call s_write_field(varname, t_step, q_cons_vf(qbmm_idx%ps(i)), x_beg, x_end, y_beg, y_end, z_beg, z_end)
                end do
 
                ! nM
                do i = 1, nb
                    write (varname, '(A,I3.3)') 'nM', i
                    call s_write_field(varname, t_step, q_cons_vf(qbmm_idx%ms(i)), x_beg, x_end, y_beg, y_end, z_beg, z_end)
                end do
            end if
 
            ! number density
            if (adv_n) then
                write (varname, '(A)') 'n'
                call s_write_field(varname, t_step, q_cons_vf(eqn_idx%n), x_beg, x_end, y_beg, y_end, z_beg, z_end)
            end if
        end if
 
        if (bubbles_lagrange) then
            ! Void fraction field
            out%q_sf(:,:,:) = 1._wp - q_cons_vf(beta_idx)%sf(-offset_x%beg:m + offset_x%end,-offset_y%beg:n + offset_y%end, &
                     & -offset_z%beg:p + offset_z%end)
            write (varname, '(A)') 'voidFraction'
            call s_write_field(varname, t_step)
 
            if (lag_txt_wrt) call s_write_lag_bubbles_results_to_text(t_step)  ! text output
            if (lag_db_wrt) call s_write_lag_bubbles_to_formatted_database_file(t_step)  ! silo file output
        end if
 
        if (ib_state_wrt) call s_write_ib_bodies_to_formatted_database_file(t_step)
 
        if (sim_data .and. proc_rank == 0) then
            call s_close_intf_data_file()
            call s_close_energy_data_file()
        end if
 
        call s_close_formatted_database_file()
 
    end subroutine s_save_data
 
    !> Fill out%q_sf from src (if given), write varname to the database, and clear varname.
    !! @param varname  field name (set by caller); blanked on return
    !! @param t_step   current time step
    !! @param src      optional scalar_field to slice into out%q_sf
    !! @param x_beg, x_end, y_beg, y_end, z_beg, z_end  output region bounds (required if src present)
    impure subroutine s_write_field(varname, t_step, src, x_beg, x_end, y_beg, y_end, z_beg, z_end)
 
        character(LEN=name_len), intent(inout)   :: varname
        integer, intent(in)                      :: t_step
        type(scalar_field), intent(in), optional :: src
        integer, intent(in), optional            :: x_beg, x_end, y_beg, y_end, z_beg, z_end
 
        if (present(src)) then
    if (.not. (present(x_beg) .and. present(x_end) .and. present(y_beg) .and. present(y_end) .and. present(z_beg)  .and. present(z_end))) then
# 688 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
        call s_mpi_abort("m_start_up.fpp:688: " // .and..and..and..and..and."Assertion failed: present(x_beg)  present(x_end)  present(y_beg)  present(y_end)  present(z_beg)   present(z_end). " &
# 688 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
                         & // "s_write_field: src requires all six output bounds")
# 688 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
    end if
# 690 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
            out%q_sf(:,:,:) = src%sf(x_beg:x_end,y_beg:y_end,z_beg:z_end)
        end if
        call s_write_variable_to_formatted_database_file(varname, t_step)
        varname(:) = ' '
 
    end subroutine s_write_field
 
    !> Transpose 3-D complex data from x-pencil to y-pencil layout via MPI_Alltoall.
    subroutine s_mpi_transpose_x2y
 
        complex(c_double_complex), allocatable :: sendbuf(:), recvbuf(:)
        integer                                :: dest_rank, src_rank
        integer                                :: i, j, k, l
 
#ifdef MFC_MPI
        allocate (sendbuf(nx*nyloc*nzloc))
        allocate (recvbuf(nx*nyloc*nzloc))
 
        do dest_rank = 0, num_procs_y - 1
            do l = 1, nzloc
                do k = 1, nyloc
                    do j = 1, nxloc
                        sendbuf(j + (k - 1)*nxloc + (l - 1)*nxloc*nyloc + dest_rank*nxloc*nyloc*nzloc) = data_cmplx(j &
                                & + dest_rank*nxloc, k, l)
                    end do
                end do
            end do
        end do
 
        call mpi_alltoall(sendbuf, nxloc*nyloc*nzloc, mpi_c_double_complex, recvbuf, nxloc*nyloc*nzloc, mpi_c_double_complex, &
                          & mpi_comm_cart12, ierr)
 
        do src_rank = 0, num_procs_y - 1
            do l = 1, nzloc
                do k = 1, nyloc
                    do j = 1, nxloc
                        data_cmplx_y(j, k + src_rank*nyloc, &
                                     & l) = recvbuf(j + (k - 1)*nxloc + (l - 1)*nxloc*nyloc + src_rank*nxloc*nyloc*nzloc)
                    end do
                end do
            end do
        end do
 
        deallocate (sendbuf)
        deallocate (recvbuf)
#endif
 
    end subroutine s_mpi_transpose_x2y
 
    !> Transpose 3-D complex data from y-pencil to z-pencil layout via MPI_Alltoall.
    subroutine s_mpi_transpose_y2z
 
        complex(c_double_complex), allocatable :: sendbuf(:), recvbuf(:)
        integer                                :: dest_rank, src_rank
        integer                                :: j, k, l
 
#ifdef MFC_MPI
        allocate (sendbuf(ny*nxloc*nzloc))
        allocate (recvbuf(ny*nxloc*nzloc))
 
        do dest_rank = 0, num_procs_z - 1
            do l = 1, nzloc
                do j = 1, nxloc
                    do k = 1, nyloc2
                        sendbuf(k + (j - 1)*nyloc2 + (l - 1)*(nyloc2*nxloc) + dest_rank*nyloc2*nxloc*nzloc) = data_cmplx_y(j, &
                                & k + dest_rank*nyloc2, l)
                    end do
                end do
            end do
        end do
 
        call mpi_alltoall(sendbuf, nyloc2*nxloc*nzloc, mpi_c_double_complex, recvbuf, nyloc2*nxloc*nzloc, mpi_c_double_complex, &
                          & mpi_comm_cart13, ierr)
 
        do src_rank = 0, num_procs_z - 1
            do l = 1, nzloc
                do j = 1, nxloc
                    do k = 1, nyloc2
                        data_cmplx_z(j, k, &
                                     & l + src_rank*nzloc) = recvbuf(k + (j - 1)*nyloc2 + (l - 1)*(nyloc2*nxloc) &
                                     & + src_rank*nyloc2*nxloc*nzloc)
                    end do
                end do
            end do
        end do
 
        deallocate (sendbuf)
        deallocate (recvbuf)
#endif
 
    end subroutine s_mpi_transpose_y2z
 
    !> Initialize all post-process sub-modules, set up I/O pointers, and prepare FFTW plans and MPI communicators.
    impure subroutine s_initialize_modules
 
        integer :: size_n(1), inembed(1), onembed(1)
 
        call s_initialize_global_parameters_module()
        if (bubbles_euler .or. bubbles_lagrange) then
            call s_initialize_bubbles_model()
        end if
        if (num_procs > 1) then
            call s_initialize_mpi_proxy_module()
            call s_initialize_mpi_common_module(exchange_all_chemistry_temperatures_in=.true., use_rdma_transport_in=.false.)
        end if
        call s_initialize_boundary_common_module()
        call s_initialize_eos_module()
        call s_initialize_variables_conversion_module(store_mixture_fields=.true., lagrange_beta_index=beta_idx)
        call s_initialize_data_input_module()
        call s_initialize_derived_variables_module()
        call s_initialize_data_output_module()
 
        if (parallel_io .neqv. .true.) then
            s_read_data_files => s_read_serial_data_files
        else
            s_read_data_files => s_read_parallel_data_files
        end if
 
#ifdef MFC_MPI
        if (fft_wrt) then
            num_procs_x = (m_glb + 1)/(m + 1)
            num_procs_y = (n_glb + 1)/(n + 1)
            num_procs_z = (p_glb + 1)/(p + 1)
 
            nx = m_glb + 1
            ny = n_glb + 1
            nz = p_glb + 1
 
            nxloc = (m_glb + 1)/num_procs_y
            nyloc = n + 1
            nyloc2 = (n_glb + 1)/num_procs_z
            nzloc = p + 1
 
            nf = max(nx, ny, nz)
 
#ifdef MFC_DEBUG
# 825 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
    block
# 825 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
        use iso_fortran_env, only: output_unit
# 825 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
 
# 825 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
        print *, 'm_start_up.fpp:825: ', '@:ALLOCATE(data_in(Nx*Nyloc*Nzloc))'
# 825 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
 
# 825 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
        call flush (output_unit)
# 825 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
    end block
# 825 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#endif
# 825 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
    allocate (data_in(nx*nyloc*nzloc))
# 825 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
 
# 825 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
 
# 825 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#if defined(MFC_OpenACC)
# 825 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
!$acc enter data create(data_in) 
# 825 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#elif defined(MFC_OpenMP)
# 825 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
!$omp target enter data map(always,alloc:data_in) 
# 825 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#endif
# 826 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#ifdef MFC_DEBUG
# 826 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
    block
# 826 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
        use iso_fortran_env, only: output_unit
# 826 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
 
# 826 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
        print *, 'm_start_up.fpp:826: ', '@:ALLOCATE(data_out(Nx*Nyloc*Nzloc))'
# 826 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
 
# 826 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
        call flush (output_unit)
# 826 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
    end block
# 826 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#endif
# 826 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
    allocate (data_out(nx*nyloc*nzloc))
# 826 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
 
# 826 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
 
# 826 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#if defined(MFC_OpenACC)
# 826 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
!$acc enter data create(data_out) 
# 826 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#elif defined(MFC_OpenMP)
# 826 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
!$omp target enter data map(always,alloc:data_out) 
# 826 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#endif
# 827 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
 
#ifdef MFC_DEBUG
# 828 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
    block
# 828 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
        use iso_fortran_env, only: output_unit
# 828 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
 
# 828 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
        print *, 'm_start_up.fpp:828: ', '@:ALLOCATE(data_cmplx(Nx, Nyloc, Nzloc))'
# 828 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
 
# 828 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
        call flush (output_unit)
# 828 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
    end block
# 828 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#endif
# 828 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
    allocate (data_cmplx(nx, nyloc, nzloc))
# 828 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
 
# 828 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
 
# 828 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#if defined(MFC_OpenACC)
# 828 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
!$acc enter data create(data_cmplx) 
# 828 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#elif defined(MFC_OpenMP)
# 828 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
!$omp target enter data map(always,alloc:data_cmplx) 
# 828 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#endif
# 829 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#ifdef MFC_DEBUG
# 829 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
    block
# 829 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
        use iso_fortran_env, only: output_unit
# 829 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
 
# 829 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
        print *, 'm_start_up.fpp:829: ', '@:ALLOCATE(data_cmplx_y(Nxloc, Ny, Nzloc))'
# 829 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
 
# 829 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
        call flush (output_unit)
# 829 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
    end block
# 829 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#endif
# 829 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
    allocate (data_cmplx_y(nxloc, ny, nzloc))
# 829 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
 
# 829 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
 
# 829 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#if defined(MFC_OpenACC)
# 829 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
!$acc enter data create(data_cmplx_y) 
# 829 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#elif defined(MFC_OpenMP)
# 829 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
!$omp target enter data map(always,alloc:data_cmplx_y) 
# 829 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#endif
# 830 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#ifdef MFC_DEBUG
# 830 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
    block
# 830 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
        use iso_fortran_env, only: output_unit
# 830 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
 
# 830 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
        print *, 'm_start_up.fpp:830: ', '@:ALLOCATE(data_cmplx_z(Nxloc, Nyloc2, Nz))'
# 830 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
 
# 830 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
        call flush (output_unit)
# 830 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
    end block
# 830 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#endif
# 830 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
    allocate (data_cmplx_z(nxloc, nyloc2, nz))
# 830 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
 
# 830 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
 
# 830 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#if defined(MFC_OpenACC)
# 830 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
!$acc enter data create(data_cmplx_z) 
# 830 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#elif defined(MFC_OpenMP)
# 830 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
!$omp target enter data map(always,alloc:data_cmplx_z) 
# 830 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#endif
# 831 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
 
#ifdef MFC_DEBUG
# 832 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
    block
# 832 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
        use iso_fortran_env, only: output_unit
# 832 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
 
# 832 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
        print *, 'm_start_up.fpp:832: ', '@:ALLOCATE(En_real(Nxloc, Nyloc2, Nz))'
# 832 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
 
# 832 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
        call flush (output_unit)
# 832 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
    end block
# 832 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#endif
# 832 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
    allocate (en_real(nxloc, nyloc2, nz))
# 832 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
 
# 832 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
 
# 832 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#if defined(MFC_OpenACC)
# 832 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
!$acc enter data create(En_real) 
# 832 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#elif defined(MFC_OpenMP)
# 832 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
!$omp target enter data map(always,alloc:En_real) 
# 832 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#endif
# 833 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#ifdef MFC_DEBUG
# 833 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
    block
# 833 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
        use iso_fortran_env, only: output_unit
# 833 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
 
# 833 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
        print *, 'm_start_up.fpp:833: ', '@:ALLOCATE(En(Nf))'
# 833 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
 
# 833 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
        call flush (output_unit)
# 833 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
    end block
# 833 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#endif
# 833 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
    allocate (en(nf))
# 833 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
 
# 833 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
 
# 833 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#if defined(MFC_OpenACC)
# 833 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
!$acc enter data create(En) 
# 833 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#elif defined(MFC_OpenMP)
# 833 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
!$omp target enter data map(always,alloc:En) 
# 833 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
#endif
# 834 "/home/runner/work/MFC/MFC/src/post_process/m_start_up.fpp"
 
            size_n(1) = nx
            inembed(1) = nx
            onembed(1) = nx
 
            fwd_plan_x = fftw_plan_many_dft(1, size_n, nyloc*nzloc, data_in, inembed, 1, nx, data_out, onembed, 1, nx, &
                                            & fftw_forward, fftw_measure)
 
            size_n(1) = ny
            inembed(1) = ny
            onembed(1) = ny
 
            fwd_plan_y = fftw_plan_many_dft(1, size_n, nxloc*nzloc, data_out, inembed, 1, ny, data_in, onembed, 1, ny, &
                                            & fftw_forward, fftw_measure)
 
            size_n(1) = nz
            inembed(1) = nz
            onembed(1) = nz
 
            fwd_plan_z = fftw_plan_many_dft(1, size_n, nxloc*nyloc2, data_in, inembed, 1, nz, data_out, onembed, 1, nz, &
                                            & fftw_forward, fftw_measure)
 
            call mpi_cart_create(mpi_comm_world, 3, (/num_procs_x, num_procs_y, num_procs_z/), (/.true., .true., .true./), &
                                 & .false., mpi_comm_cart, ierr)
            call mpi_cart_coords(mpi_comm_cart, proc_rank, 3, cart3d_coords, ierr)
 
            call mpi_cart_sub(mpi_comm_cart, (/.true., .true., .false./), mpi_comm_cart12, ierr)
            call mpi_comm_rank(mpi_comm_cart12, proc_rank12, ierr)
            call mpi_cart_coords(mpi_comm_cart12, proc_rank12, 2, cart2d12_coords, ierr)
 
            call mpi_cart_sub(mpi_comm_cart, (/.true., .false., .true./), mpi_comm_cart13, ierr)
            call mpi_comm_rank(mpi_comm_cart13, proc_rank13, ierr)
            call mpi_cart_coords(mpi_comm_cart13, proc_rank13, 2, cart2d13_coords, ierr)
        end if
#endif
 
    end subroutine s_initialize_modules
 
    !> Perform a distributed forward 3-D FFT using pencil decomposition with FFTW and MPI transposes.
    subroutine s_mpi_fft_fwd
 
        integer :: j, k, l
 
#ifdef MFC_MPI
        do l = 1, nzloc
            do k = 1, nyloc
                do j = 1, nx
                    data_in(j + (k - 1)*nx + (l - 1)*nx*nyloc) = data_cmplx(j, k, l)
                end do
            end do
        end do
 
        call fftw_execute_dft(fwd_plan_x, data_in, data_out)
 
        do l = 1, nzloc
            do k = 1, nyloc
                do j = 1, nx
                    data_cmplx(j, k, l) = data_out(j + (k - 1)*nx + (l - 1)*nx*nyloc)
                end do
            end do
        end do
 
        call s_mpi_transpose_x2y !!Change Pencil from data_cmplx to data_cmpx_y
 
        do l = 1, nzloc
            do k = 1, nxloc
                do j = 1, ny
                    data_out(j + (k - 1)*ny + (l - 1)*ny*nxloc) = data_cmplx_y(k, j, l)
                end do
            end do
        end do
 
        call fftw_execute_dft(fwd_plan_y, data_out, data_in)
 
        do l = 1, nzloc
            do k = 1, nxloc
                do j = 1, ny
                    data_cmplx_y(k, j, l) = data_in(j + (k - 1)*ny + (l - 1)*ny*nxloc)
                end do
            end do
        end do
 
        call s_mpi_transpose_y2z !!Change Pencil from data_cmplx_y to data_cmpx_z
 
        do l = 1, nyloc2
            do k = 1, nxloc
                do j = 1, nz
                    data_in(j + (k - 1)*nz + (l - 1)*nz*nxloc) = data_cmplx_z(k, l, j)
                end do
            end do
        end do
 
        call fftw_execute_dft(fwd_plan_z, data_in, data_out)
 
        do l = 1, nyloc2
            do k = 1, nxloc
                do j = 1, nz
                    data_cmplx_z(k, l, j) = data_out(j + (k - 1)*nz + (l - 1)*nz*nxloc)
                end do
            end do
        end do
#endif
 
    end subroutine s_mpi_fft_fwd
 
    !> Set up the MPI environment, read and broadcast user inputs, and decompose the computational domain.
    impure subroutine s_initialize_mpi_domain
 
        type(int_bounds_info), dimension(3) :: output_offsets
 
        num_dims = 1 + min(1, n) + min(1, p)
 
        call s_mpi_initialize()
 
        if (proc_rank == 0) then
            call s_assign_default_values_to_user_inputs()
            call s_read_input_file()
            call s_check_input_file()
 
            print '(" Post-processing a ", I0, "x", I0, "x", I0, " case on ", I0, " rank(s)")', m, n, p, num_procs
        end if
 
        call s_mpi_bcast_user_inputs()
 
        ! Save original BCs before decomposition overwrites them with MPI neighbor ranks
        ib_bc_x = bc_x
        ib_bc_y = bc_y
        ib_bc_z = bc_z
 
        call s_initialize_parallel_io()
        output_offsets = (/offset_x, offset_y, offset_z/)
        call s_mpi_decompose_computational_domain(write_silo_ghost_offsets=format == format_silo, adjust_local_domains=.false., &
            & output_offsets=output_offsets)
        offset_x = output_offsets(1)
        offset_y = output_offsets(2)
        offset_z = output_offsets(3)
        call s_check_inputs_fft()
 
        bc = bc_xyz_info(bc_x, bc_y, bc_z)
 
    end subroutine s_initialize_mpi_domain
 
    !> Destroy FFTW plans, free MPI communicators, and finalize all post-process sub-modules.
    impure subroutine s_finalize_modules
 
        s_read_data_files => null()
 
        if (fft_wrt) then
            if (c_associated(fwd_plan_x)) call fftw_destroy_plan(fwd_plan_x)
            if (c_associated(fwd_plan_y)) call fftw_destroy_plan(fwd_plan_y)
            if (c_associated(fwd_plan_z)) call fftw_destroy_plan(fwd_plan_z)
            if (allocated(data_in)) deallocate (data_in)
            if (allocated(data_out)) deallocate (data_out)
            if (allocated(data_cmplx)) deallocate (data_cmplx)
            if (allocated(data_cmplx_y)) deallocate (data_cmplx_y)
            if (allocated(data_cmplx_z)) deallocate (data_cmplx_z)
            if (allocated(en_real)) deallocate (en_real)
            if (allocated(en)) deallocate (en)
            call fftw_cleanup()
        end if
 
#ifdef MFC_MPI
        if (fft_wrt) then
            if (mpi_comm_cart12 /= mpi_comm_null) call mpi_comm_free(mpi_comm_cart12, ierr)
            if (mpi_comm_cart13 /= mpi_comm_null) call mpi_comm_free(mpi_comm_cart13, ierr)
            if (mpi_comm_cart /= mpi_comm_null) call mpi_comm_free(mpi_comm_cart, ierr)
        end if
#endif
 
        call s_finalize_data_output_module()
        call s_finalize_derived_variables_module()
        call s_finalize_data_input_module()
        call s_finalize_variables_conversion_module()
        call s_finalize_eos_module()
        if (num_procs > 1) then
            call s_finalize_mpi_proxy_module()
            call s_finalize_mpi_common_module()
        end if
        call s_finalize_global_parameters_module()
 
        call s_mpi_finalize()
 
    end subroutine s_finalize_modules
 
end module m_start_up