diff --git a/CMakeLists.txt b/CMakeLists.txt
index 751a690b5..246f03ab7 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -460,8 +460,10 @@ function(MFC_SETUP_TARGET)
endif()
if (ARGS_HDF5)
- find_package(HDF5 REQUIRED)
+ find_package(HDF5 REQUIRED COMPONENTS Fortran)
target_link_libraries(${a_target} PRIVATE HDF5::HDF5)
+ target_include_directories(${a_target} PRIVATE ${HDF5_Fortran_INCLUDE_DIRS})
+ target_link_libraries(${a_target} PRIVATE ${HDF5_Fortran_LIBRARIES})
endif()
if (ARGS_FFTW)
diff --git a/docs/documentation/case.md b/docs/documentation/case.md
index 98e75b326..4fa420320 100644
--- a/docs/documentation/case.md
+++ b/docs/documentation/case.md
@@ -638,7 +638,10 @@ To restart the simulation from $k$-th time step, see [Restarting Cases](running.
The table lists formatted database output parameters. The parameters define variables that are outputted from simulation and file types and formats of data as well as options for post-processing.
-- `format` specifies the choice of the file format of data file outputted by MFC by an integer of 1 and 2. `format = 1` and `2` correspond to Silo-HDF5 format and binary format, respectively.
+- `format` specifies the choice of the file format of data file outputted by MFC by an integer of 1, 2, or 3:
+ - `format = 1`: Silo-HDF5 format (default). Works with VisIt and ParaView 5.11 and earlier.
+ - `format = 2`: Binary format. Lightweight output for custom analysis pipelines.
+ - `format = 3`: Native HDF5 + XDMF format. Recommended for ParaView 5.12 and later compatibility. Produces portable, self-describing output without external library dependencies beyond HDF5.
- `precision` specifies the choice of the floating-point format of the data file outputted by MFC by an integer of 1 and 2. `precision = 1` and `2` correspond to single-precision and double-precision formats, respectively.
diff --git a/docs/documentation/testing.md b/docs/documentation/testing.md
index e139c5b20..1907d57e1 100644
--- a/docs/documentation/testing.md
+++ b/docs/documentation/testing.md
@@ -97,8 +97,21 @@ To test the post-processing code, append the `-a` or `--test-all` option:
./mfc.sh test -a -j 8
```
-This argument will re-run the test stack with `parallel_io='T'`, which generates silo_hdf5 files.
+This argument will re-run the test stack with `parallel_io='T'`, which generates output files depending on the `format` parameter.
It will also turn most write parameters (`*_wrt`) on.
-Then, it searches through the silo files using `h5dump` to ensure that there are no `NaN`s or `Infinity`s.
-Although adding this option does not guarantee that accurate `.silo` files are generated, it does ensure that the post-process code does not fail or produce malformed data.
+
+The test framework validates post-process output for the following formats:
+
+| Format | Output Directory | Validation |
+|--------|------------------|------------|
+| `format=1` | `silo_hdf5/` | Checks HDF5 files for NaN/Infinity values |
+| `format=3` | `hdf5_xdmf/` | Checks HDF5 files for NaN/Infinity, then validates XDMF metadata structure |
+
+For Silo-HDF5 output (`format=1`), the test searches through the `.silo` files using `h5dump` to ensure that there are no `NaN`s or `Infinity`s.
+
+For native HDF5+XDMF output (`format=3`), the test:
+1. Validates all `.h5` files contain no `NaN` or `Infinity` values
+2. Validates `.xdmf` metadata files have correct XML structure with required elements (`Domain`, `Grid`, `Topology`, `Geometry`, `Attribute`)
+
+Although adding this option does not guarantee that fully accurate output files are generated, it does ensure that the post-process code does not fail or produce malformed data.
diff --git a/docs/documentation/visualization.md b/docs/documentation/visualization.md
index dd9c5a295..86fbe7f86 100644
--- a/docs/documentation/visualization.md
+++ b/docs/documentation/visualization.md
@@ -1,5 +1,70 @@
# Flow visualization
+MFC supports three output formats for post-processed simulation data. You may choose the format that best matches your visualization workflow:
+
+| Format | Parameter | Output Directory | Best For |
+|--------|-----------|------------------|----------|
+| Silo-HDF5 | `format = 1` | `silo_hdf5/` | VisIt, ParaView ≤ 5.11 |
+| Binary | `format = 2` | `D/` | Custom analysis scripts |
+| HDF5 + XDMF | `format = 3` | `hdf5_xdmf/` | ParaView 5.12-6.0+, portable workflows |
+
+## Output Format Details
+
+### Format 1: Silo-HDF5 (Default)
+
+Silo-HDF5 is the traditional MFC output format. It requires the Silo library and produces databases compatible with VisIt and older versions of ParaView.
+
+**Directory structure:**
+```
+case_dir/silo_hdf5/
+├── root/ # Index files for visualization
+│ └── *.silo # One per time step
+└── p0/, p1/, ... # Per-processor data files
+```
+
+**Pros:** Well-established, works with VisIt
+**Cons:** Requires separate Silo library installation; ParaView 5.12-6.0.1 Silo reader is broken
+
+### Format 2: Binary
+
+Binary format writes raw floating-point data for custom analysis pipelines.
+
+**Directory structure:**
+```
+case_dir/D/
+└── *.dat # Raw data files
+```
+
+**Pros:** Lightweight, easy to parse programmatically
+**Cons:** No built-in visualization support; requires custom readers
+
+### Format 3: HDF5 + XDMF (Recommended for ParaView 5.12+)
+
+Native HDF5 output with XDMF metadata provides a portable, self-describing format that works reliably with modern visualization tools without external library dependencies beyond HDF5.
+
+**Directory structure:**
+```
+case_dir/hdf5_xdmf/
+├── root/ # XDMF metadata files
+│ ├── timestep_0.xdmf # Open these in ParaView
+│ ├── timestep_100.xdmf
+│ └── ...
+└── p0/, p1/, ... # Per-processor HDF5 data
+ └── *.h5
+```
+
+**Pros:**
+- Works with ParaView 5.12+ and later without Silo dependencies
+- Fully portable across systems and software versions
+- Self-describing via XDMF metadata
+- Supports all MFC features: grid data, MHD fields, and Lagrangian bubble particles
+
+**Cons:** Larger file sizes than Silo-HDF5 (internal compression not implemented yet)
+
+**Usage:** Set `format = 3` in your post_process input file, then open the `.xdmf` files from `hdf5_xdmf/root/` in ParaView.
+
+---
+
A post-processed database in Silo-HDF5 format can be visualized and analyzed using Paraview and VisIt.
After the post-processing of simulation data (see section [Running](running.md)), a directory named `silo_hdf5` contains a silo-HDF5 database.
Here, `silo_hdf5/` includes a directory named `root/` that contains index files for flow field data at each saved time step.
@@ -17,6 +82,27 @@ The figure below shows the iso-contour of the liquid void fraction (`alpha1`) in
+### Visualizing HDF5+XDMF with ParaView 5.12+
+
+ParaView 5.12-6.0.1 have issues with their ViSiT Silo plugin. Use `format = 3` (HDF5+XDMF) for post-processing:
+
+1. Run post-processing with `format = 3` in your input file
+2. Navigate to the `hdf5_xdmf/root/` directory
+3. Open the `.xdmf` file for your desired time step (e.g., `timestep_0.xdmf`), or open the master time series file for aggregated time steps
+4. ParaView will automatically load the grid and all available flow variables
+
+**Opening Multiple Time Steps:**
+Each XDMF file represents a single time step. To animate through time:
+- Open the master time series XDMF file if generated
+- Or, use the XDMF reader's time support if your ParaView version supports it
+
+**Lagrangian Bubble Visualization:**
+When simulations include Lagrangian bubbles (`lagrange_bubbles = .true.`), the XDMF file contains both:
+- **Grid data:** Eulerian field variables on the computational mesh
+- **Particle data:** Bubble positions and properties as point clouds
+
+In ParaView, the particle data appears as a separate "Particles" grid that can be visualized using the Glyph filter to render spheres at bubble locations.
+
### Visualizing data in cylindrical coordinates
Visualizing data in cylindrical coordinates requires a coordinate transformation of the raw data in the database file.
diff --git a/src/post_process/m_data_output.fpp b/src/post_process/m_data_output.fpp
index 410ab9ebc..2867996b4 100644
--- a/src/post_process/m_data_output.fpp
+++ b/src/post_process/m_data_output.fpp
@@ -4,9 +4,10 @@
!> @brief This module enables the restructuring of the raw simulation data
!! file(s) into formatted database file(s). The formats that may be
-!! chosen from include Silo-HDF5 and Binary. Each of these database
-!! structures contains information about the grid as well as each of
-!! the flow variable(s) that were chosen by the user to be included.
+!! chosen from include Silo-HDF5, Binary, and HDF5+XDMF. Each of
+!! these database structures contains information about the grid as
+!! well as each of the flow variable(s) that were chosen by the user
+!! to be included.
module m_data_output
use m_derived_types ! Definitions of the derived types
@@ -23,6 +24,8 @@ module m_data_output
use m_variables_conversion
+ use m_hdf5_xdmf_output ! Native HDF5+XDMF output support
+
implicit none
private; public :: s_initialize_data_output_module, &
@@ -37,6 +40,7 @@ module m_data_output
s_write_intf_data_file, &
s_write_energy_data_file, &
s_close_formatted_database_file, &
+ s_write_xdmf_metadata_file, &
s_close_intf_data_file, &
s_close_energy_data_file, &
s_finalize_data_output_module
@@ -253,7 +257,7 @@ contains
! Generating Binary Directory Tree
- else
+ elseif (format == 2) then
! Creating the directory associated with the local process
dbdir = trim(case_dir)//'/binary'
@@ -285,6 +289,13 @@ contains
end if
+ ! Generating HDF5+XDMF Directory Tree
+
+ elseif (format == 3) then
+
+ ! Initialize the HDF5+XDMF module (creates directory structure)
+ call s_initialize_hdf5_xdmf_module()
+
end if
if (bubbles_lagrange) then !Lagrangian solver
@@ -539,7 +550,7 @@ contains
! Binary Database Format
- else
+ elseif (format == 2) then
! Generating the relative path to the formatted database slave
! file, that is to be opened for the current time-step, t_step
@@ -597,6 +608,13 @@ contains
end if
+ ! HDF5+XDMF Database Format
+
+ elseif (format == 3) then
+
+ ! Open new HDF5 file for this time-step
+ call s_open_hdf5_file(t_step)
+
end if
end subroutine s_open_formatted_database_file
@@ -836,6 +854,13 @@ contains
end if
+ ! HDF5+XDMF Database Format
+
+ elseif (format == 3) then
+
+ ! Write grid coordinates to HDF5 file
+ call s_write_grid_to_hdf5()
+
end if
end subroutine s_write_grid_to_formatted_database_file
@@ -1003,7 +1028,7 @@ contains
! Binary Database Format
- else
+ elseif (format == 2) then
! Writing the name of the flow variable and its data, associated
! with the local processor, to the formatted database slave file
@@ -1034,6 +1059,13 @@ contains
end if
+ ! HDF5+XDMF Database Format
+
+ elseif (format == 3) then
+
+ ! Write the flow variable to HDF5 file
+ call s_write_variable_to_hdf5(varname, q_sf)
+
end if
end subroutine s_write_variable_to_formatted_database_file
@@ -1338,23 +1370,31 @@ contains
! extents in the formatted database master file. In addition, it
! also records a sub-domain connectivity map so that the entire
! grid may be reassembled by looking at the master file.
- if (proc_rank == 0) then
- do i = 1, num_procs
- write (meshnames(i), '(A,I0,A,I0,A)') '../p', i - 1, &
- '/', t_step, '.silo:lag_bubbles'
- meshtypes(i) = DB_POINTMESH
- end do
- err = DBSET2DSTRLEN(len(meshnames(1)))
- err = DBPUTMMESH(dbroot, 'lag_bubbles', 16, &
- num_procs, meshnames, &
- len_trim(meshnames), &
- meshtypes, DB_F77NULL, ierr)
- end if
+ ! Format-specific particle mesh output
+ if (format == 1) then
+ ! Silo-HDF5 format
+ if (proc_rank == 0) then
+ do i = 1, num_procs
+ write (meshnames(i), '(A,I0,A,I0,A)') '../p', i - 1, &
+ '/', t_step, '.silo:lag_bubbles'
+ meshtypes(i) = DB_POINTMESH
+ end do
+ err = DBSET2DSTRLEN(len(meshnames(1)))
+ err = DBPUTMMESH(dbroot, 'lag_bubbles', 16, &
+ num_procs, meshnames, &
+ len_trim(meshnames), &
+ meshtypes, DB_F77NULL, ierr)
+ end if
+
+ err = DBPUTPM(dbfile, 'lag_bubbles', 11, 3, &
+ px, py, pz, nBub, &
+ DB_DOUBLE, DB_F77NULL, ierr)
- err = DBPUTPM(dbfile, 'lag_bubbles', 11, 3, &
- px, py, pz, nBub, &
- DB_DOUBLE, DB_F77NULL, ierr)
+ elseif (format == 3) then
+ ! HDF5+XDMF format - write particle mesh to HDF5
+ call s_write_particle_mesh_to_hdf5(px, py, pz, nBub)
+ end if
if (lag_id_wrt) call s_write_lag_variable_to_formatted_database_file('part_id', t_step, bub_id, nBub)
if (lag_vel_wrt) then
@@ -1395,24 +1435,28 @@ contains
call MPI_FILE_CLOSE(ifile, ierr)
call MPI_TYPE_FREE(view, ierr)
- if (proc_rank == 0) then
+ ! Format-specific empty particle mesh handling
+ if (format == 1) then
+ ! Silo-HDF5 format
+ if (proc_rank == 0) then
+ do i = 1, num_procs
+ write (meshnames(i), '(A,I0,A,I0,A)') '../p', i - 1, &
+ '/', t_step, '.silo:lag_bubbles'
+ meshtypes(i) = DB_POINTMESH
+ end do
+ err = DBSET2DSTRLEN(len(meshnames(1)))
+ err = DBPUTMMESH(dbroot, 'lag_bubbles', 16, &
+ num_procs, meshnames, &
+ len_trim(meshnames), &
+ meshtypes, DB_F77NULL, ierr)
+ end if
- do i = 1, num_procs
- write (meshnames(i), '(A,I0,A,I0,A)') '../p', i - 1, &
- '/', t_step, '.silo:lag_bubbles'
- meshtypes(i) = DB_POINTMESH
- end do
- err = DBSET2DSTRLEN(len(meshnames(1)))
- err = DBPUTMMESH(dbroot, 'lag_bubbles', 16, &
- num_procs, meshnames, &
- len_trim(meshnames), &
- meshtypes, DB_F77NULL, ierr)
+ err = DBSETEMPTYOK(1)
+ err = DBPUTPM(dbfile, 'lag_bubbles', 11, 3, &
+ dummy_data, dummy_data, dummy_data, 0, &
+ DB_DOUBLE, DB_F77NULL, ierr)
end if
-
- err = DBSETEMPTYOK(1)
- err = DBPUTPM(dbfile, 'lag_bubbles', 11, 3, &
- dummy_data, dummy_data, dummy_data, 0, &
- DB_DOUBLE, DB_F77NULL, ierr)
+ ! Note: For format=3 with no particles, we don't write particle data
if (lag_id_wrt) call s_write_lag_variable_to_formatted_database_file('part_id', t_step)
if (lag_vel_wrt) then
@@ -1454,39 +1498,55 @@ contains
dummy_data = 0._wp
if (present(nBubs) .and. present(data)) then
- if (proc_rank == 0) then
- do i = 1, num_procs
- write (var_names(i), '(A,I0,A,I0,A)') '../p', i - 1, &
- '/', t_step, '.silo:'//trim(varname)
- var_types(i) = DB_POINTVAR
- end do
- err = DBSET2DSTRLEN(len(var_names(1)))
- err = DBPUTMVAR(dbroot, trim(varname), len_trim(varname), &
- num_procs, var_names, &
- len_trim(var_names), &
- var_types, DB_F77NULL, ierr)
+ ! Has particle data
+
+ if (format == 1) then
+ ! Silo-HDF5 format
+ if (proc_rank == 0) then
+ do i = 1, num_procs
+ write (var_names(i), '(A,I0,A,I0,A)') '../p', i - 1, &
+ '/', t_step, '.silo:'//trim(varname)
+ var_types(i) = DB_POINTVAR
+ end do
+ err = DBSET2DSTRLEN(len(var_names(1)))
+ err = DBPUTMVAR(dbroot, trim(varname), len_trim(varname), &
+ num_procs, var_names, &
+ len_trim(var_names), &
+ var_types, DB_F77NULL, ierr)
+ end if
+
+ err = DBPUTPV1(dbfile, trim(varname), len_trim(varname), &
+ 'lag_bubbles', 11, data, nBubs, DB_DOUBLE, DB_F77NULL, ierr)
+
+ elseif (format == 3) then
+ ! HDF5+XDMF format
+ call s_write_particle_variable_to_hdf5(varname, data, nBubs)
end if
- err = DBPUTPV1(dbfile, trim(varname), len_trim(varname), &
- 'lag_bubbles', 11, data, nBubs, DB_DOUBLE, DB_F77NULL, ierr)
else
- if (proc_rank == 0) then
- do i = 1, num_procs
- write (var_names(i), '(A,I0,A,I0,A)') '../p', i - 1, &
- '/', t_step, '.silo:'//trim(varname)
- var_types(i) = DB_POINTVAR
- end do
- err = DBSET2DSTRLEN(len(var_names(1)))
+ ! No particle data (empty case)
+
+ if (format == 1) then
+ ! Silo-HDF5 format
+ if (proc_rank == 0) then
+ do i = 1, num_procs
+ write (var_names(i), '(A,I0,A,I0,A)') '../p', i - 1, &
+ '/', t_step, '.silo:'//trim(varname)
+ var_types(i) = DB_POINTVAR
+ end do
+ err = DBSET2DSTRLEN(len(var_names(1)))
+ err = DBSETEMPTYOK(1)
+ err = DBPUTMVAR(dbroot, trim(varname), len_trim(varname), &
+ num_procs, var_names, &
+ len_trim(var_names), &
+ var_types, DB_F77NULL, ierr)
+ end if
+
err = DBSETEMPTYOK(1)
- err = DBPUTMVAR(dbroot, trim(varname), len_trim(varname), &
- num_procs, var_names, &
- len_trim(var_names), &
- var_types, DB_F77NULL, ierr)
+ err = DBPUTPV1(dbfile, trim(varname), len_trim(varname), &
+ 'lag_bubbles', 11, dummy_data, 0, DB_DOUBLE, DB_F77NULL, ierr)
end if
-
- err = DBSETEMPTYOK(1)
- err = DBPUTPV1(dbfile, trim(varname), len_trim(varname), &
- 'lag_bubbles', 11, dummy_data, 0, DB_DOUBLE, DB_F77NULL, ierr)
+ ! Note: For format=3 with no particles, no data written
end if
end subroutine s_write_lag_variable_to_formatted_database_file
@@ -1702,10 +1762,14 @@ contains
if (proc_rank == 0) ierr = DBCLOSE(dbroot)
! Binary database format
- else
+ elseif (format == 2) then
close (dbfile)
if (n == 0 .and. proc_rank == 0) close (dbroot)
+ ! HDF5+XDMF database format
+ elseif (format == 3) then
+ call s_close_hdf5_file()
+
end if
end subroutine s_close_formatted_database_file
@@ -1722,6 +1786,21 @@ contains
end subroutine s_close_energy_data_file
+ !> Subroutine to write XDMF metadata file for HDF5+XDMF format
+ !! @param t_step Current time step
+ !! @param time Physical time value
+ impure subroutine s_write_xdmf_metadata_file(t_step, time)
+
+ integer, intent(in) :: t_step
+ real(wp), intent(in) :: time
+
+ ! Only write XDMF for format=3
+ if (format == 3) then
+ call s_write_xdmf_file(t_step, time)
+ end if
+
+ end subroutine s_write_xdmf_metadata_file
+
impure subroutine s_finalize_data_output_module()
! Description: Deallocation procedures for the module
@@ -1746,6 +1825,13 @@ contains
deallocate (dims)
end if
+ ! Finalize HDF5+XDMF module if used
+ if (format == 3) then
+ ! Write master time series XDMF file that references all timesteps
+ call s_write_xdmf_time_series(t_step_start, t_step_stop, t_step_save)
+ call s_finalize_hdf5_xdmf_module()
+ end if
+
end subroutine s_finalize_data_output_module
end module m_data_output
diff --git a/src/post_process/m_hdf5_xdmf_output.fpp b/src/post_process/m_hdf5_xdmf_output.fpp
new file mode 100644
index 000000000..205a9ef02
--- /dev/null
+++ b/src/post_process/m_hdf5_xdmf_output.fpp
@@ -0,0 +1,745 @@
+!>
+!! @file m_hdf5_xdmf_output.fpp
+!! @brief Contains module m_hdf5_xdmf_output
+!!
+!! @details This module implements native HDF5 output with XDMF metadata for
+!! visualization in ParaView 5.12+ and other XDMF-compatible tools.
+!! It provides a portable alternative to Silo-HDF5 format that works
+!! reliably across different visualization software versions.
+!!
+!! @par Output Structure
+!! The module creates the following directory structure:
+!! @verbatim
+!! case_dir/hdf5_xdmf/
+!! ├── p0/ # Per-processor HDF5 files
+!! │ ├── 0.h5 # Time step 0
+!! │ └── 100.h5 # Time step 100
+!! ├── p1/ # Processor 1 files (for parallel runs)
+!! │ └── ...
+!! └── root/ # XDMF metadata files (written by rank 0)
+!! ├── timestep_0.xdmf
+!! └── timestep_100.xdmf
+!! @endverbatim
+!!
+!! @par Supported Data Types
+!! - Structured grid (rectilinear mesh) data for 1D/2D/3D simulations
+!! - Particle/point mesh data for Lagrangian bubble tracking
+!! - All flow variables (density, velocity, pressure, etc.)
+!! - MHD magnetic field components (Bx, By, Bz)
+!! - Chemistry species mass fractions
+!!
+!! @par Usage
+!! Set format=3 in the post_process input file to enable HDF5+XDMF output.
+!! Open the .xdmf files from the root/ folder in ParaView to visualize.
+!!
+
+!> @brief This module provides native HDF5 + XDMF output support for MFC.
+!!
+!! XDMF (eXtensible Data Model and Format) is an XML-based metadata format that
+!! describes HDF5 data for visualization tools like ParaView.
+!!
+!! @note This module requires HDF5 with Fortran bindings (hdf5-mpi recommended).
+module m_hdf5_xdmf_output
+
+ use m_derived_types
+ use m_global_parameters
+ use m_mpi_proxy
+ use m_compile_specific
+ use hdf5
+
+ implicit none
+
+ private; public :: s_initialize_hdf5_xdmf_module, &
+ s_open_hdf5_file, &
+ s_write_grid_to_hdf5, &
+ s_write_variable_to_hdf5, &
+ s_write_particle_mesh_to_hdf5, &
+ s_write_particle_variable_to_hdf5, &
+ s_close_hdf5_file, &
+ s_write_xdmf_file, &
+ s_write_xdmf_time_series, &
+ s_finalize_hdf5_xdmf_module
+
+ !> @name HDF5 Identifiers
+ !! @{
+ integer(HID_T) :: h5file_id !< HDF5 file identifier
+ integer(HID_T) :: h5dspace_id !< HDF5 dataspace identifier
+ integer(HID_T) :: h5dset_id !< HDF5 dataset identifier
+ integer(HID_T) :: h5plist_id !< HDF5 property list identifier
+ integer(HID_T) :: h5type_id !< HDF5 datatype (single or double precision)
+ !> @}
+
+ !> @name Directory Paths
+ !! @{
+ character(LEN=path_len + name_len) :: hdf5_dir !< Base HDF5 output directory
+ character(LEN=path_len + 2*name_len) :: proc_rank_dir_hdf5 !< Per-processor directory
+ character(LEN=path_len + 2*name_len) :: rootdir_hdf5 !< Root directory for XDMF files
+ !> @}
+
+ integer :: h5err !< HDF5 error flag
+
+ !> @name Variable Tracking for XDMF
+ !! These arrays store variable names and particle information for XDMF generation.
+ !! @{
+ character(LEN=name_len), allocatable, dimension(:) :: stored_varnames !< Grid variable names
+ integer :: num_stored_vars !< Number of stored grid variables
+
+ character(LEN=name_len), allocatable, dimension(:) :: stored_particle_varnames !< Particle variable names
+ integer :: num_stored_particle_vars !< Number of stored particle variables
+ integer :: stored_num_particles !< Number of particles for current timestep
+ logical :: has_particle_data !< Flag indicating particle data exists
+ !> @}
+
+contains
+
+ !> @brief Initialize the HDF5+XDMF output module
+ !!
+ !! This subroutine initializes the HDF5 library and creates the directory
+ !! structure for output files. It must be called before any other routines
+ !! in this module.
+ !!
+ !! @par Directory Structure Created
+ !! - hdf5_xdmf/p{rank}/ for each MPI rank
+ !! - hdf5_xdmf/root/ for XDMF metadata (rank 0 only)
+ subroutine s_initialize_hdf5_xdmf_module()
+
+ character(LEN=len_trim(case_dir) + 2*name_len) :: file_loc
+ logical :: dir_check
+ integer :: i
+
+ ! Initialize HDF5 library
+ call h5open_f(h5err)
+ if (h5err /= 0) then
+ call s_mpi_abort('Failed to initialize HDF5 library. Exiting.')
+ end if
+
+ ! Create HDF5 output directory structure
+ hdf5_dir = trim(case_dir)//'/hdf5_xdmf'
+
+ write (proc_rank_dir_hdf5, '(A,I0)') '/p', proc_rank
+ proc_rank_dir_hdf5 = trim(hdf5_dir)//trim(proc_rank_dir_hdf5)
+
+ file_loc = trim(proc_rank_dir_hdf5)//'/.'
+ call my_inquire(file_loc, dir_check)
+ if (.not. dir_check) then
+ call s_create_directory(trim(proc_rank_dir_hdf5))
+ end if
+
+ ! Root process creates root directory for XDMF files
+ if (proc_rank == 0) then
+ rootdir_hdf5 = trim(hdf5_dir)//'/root'
+ file_loc = trim(rootdir_hdf5)//'/.'
+ call my_inquire(file_loc, dir_check)
+ if (.not. dir_check) then
+ call s_create_directory(trim(rootdir_hdf5))
+ end if
+ end if
+
+ ! Allocate variable name storage for XDMF generation
+ allocate (stored_varnames(200)) ! Max 200 grid variables
+ num_stored_vars = 0
+
+ ! Allocate particle variable name storage
+ allocate (stored_particle_varnames(50)) ! Max 50 particle variables
+ num_stored_particle_vars = 0
+ stored_num_particles = 0
+ has_particle_data = .false.
+
+ ! Set HDF5 datatype based on precision setting
+ ! precision=1 -> single (32-bit), precision=2 -> double (64-bit)
+ if (precision == 1) then
+ h5type_id = H5T_NATIVE_REAL
+ else
+ h5type_id = H5T_NATIVE_DOUBLE
+ end if
+
+ end subroutine s_initialize_hdf5_xdmf_module
+
+ !> @brief Open an HDF5 file for the current time step
+ !!
+ !! Creates a new HDF5 file for the specified time step in the processor's
+ !! output directory. Resets variable counters for XDMF metadata generation.
+ !!
+ !! @param[in] t_step Current time step number
+ subroutine s_open_hdf5_file(t_step)
+
+ integer, intent(in) :: t_step
+ character(LEN=len_trim(case_dir) + 3*name_len) :: file_loc
+
+ ! Generate file path
+ write (file_loc, '(A,I0,A)') '/', t_step, '.h5'
+ file_loc = trim(proc_rank_dir_hdf5)//trim(file_loc)
+
+ ! Create HDF5 file with default properties
+ call h5fcreate_f(trim(file_loc), H5F_ACC_TRUNC_F, h5file_id, h5err)
+ if (h5err /= 0) then
+ call s_mpi_abort('Unable to create HDF5 file '//trim(file_loc)//'. Exiting.')
+ end if
+
+ ! Reset variable counters for this timestep
+ num_stored_vars = 0
+ num_stored_particle_vars = 0
+ stored_num_particles = 0
+ has_particle_data = .false.
+
+ end subroutine s_open_hdf5_file
+
+ !> @brief Write grid coordinates to HDF5 file
+ !!
+ !! Writes the cell boundary coordinates (x_cb, y_cb, z_cb) to the HDF5 file.
+ !! These define the rectilinear mesh geometry for visualization.
+ !! Grid coordinates are always written in double precision for accuracy.
+ !! Also writes spatial extent metadata as HDF5 attributes for optimization.
+ subroutine s_write_grid_to_hdf5()
+
+ integer(HSIZE_T), dimension(1) :: dims_1d
+ real(wp) :: spatial_extent(6) ! xmin, xmax, ymin, ymax, zmin, zmax
+ integer :: i
+
+ ! Write X coordinates (always double precision for grid accuracy)
+ dims_1d(1) = size(x_cb)
+ call h5screate_simple_f(1, dims_1d, h5dspace_id, h5err)
+ call h5dcreate_f(h5file_id, 'x_cb', H5T_NATIVE_DOUBLE, h5dspace_id, h5dset_id, h5err)
+ call h5dwrite_f(h5dset_id, H5T_NATIVE_DOUBLE, x_cb, dims_1d, h5err)
+ call h5dclose_f(h5dset_id, h5err)
+ call h5sclose_f(h5dspace_id, h5err)
+
+ ! Write Y coordinates if 2D or 3D
+ if (n > 0) then
+ dims_1d(1) = size(y_cb)
+ call h5screate_simple_f(1, dims_1d, h5dspace_id, h5err)
+ call h5dcreate_f(h5file_id, 'y_cb', H5T_NATIVE_DOUBLE, h5dspace_id, h5dset_id, h5err)
+ call h5dwrite_f(h5dset_id, H5T_NATIVE_DOUBLE, y_cb, dims_1d, h5err)
+ call h5dclose_f(h5dset_id, h5err)
+ call h5sclose_f(h5dspace_id, h5err)
+ end if
+
+ ! Write Z coordinates if 3D
+ if (p > 0) then
+ dims_1d(1) = size(z_cb)
+ call h5screate_simple_f(1, dims_1d, h5dspace_id, h5err)
+ call h5dcreate_f(h5file_id, 'z_cb', H5T_NATIVE_DOUBLE, h5dspace_id, h5dset_id, h5err)
+ call h5dwrite_f(h5dset_id, H5T_NATIVE_DOUBLE, z_cb, dims_1d, h5err)
+ call h5dclose_f(h5dset_id, h5err)
+ call h5sclose_f(h5dspace_id, h5err)
+ end if
+
+ ! Write spatial extent as file-level attribute for visualization optimization
+ ! Format: [xmin, xmax, ymin, ymax, zmin, zmax]
+ spatial_extent(1) = x_cb(-1) ! xmin
+ spatial_extent(2) = x_cb(m) ! xmax
+ if (n > 0) then
+ spatial_extent(3) = y_cb(-1) ! ymin
+ spatial_extent(4) = y_cb(n) ! ymax
+ else
+ spatial_extent(3) = 0._wp
+ spatial_extent(4) = 0._wp
+ end if
+ if (p > 0) then
+ spatial_extent(5) = z_cb(-1) ! zmin
+ spatial_extent(6) = z_cb(p) ! zmax
+ else
+ spatial_extent(5) = 0._wp
+ spatial_extent(6) = 0._wp
+ end if
+ call s_write_attribute_real_array(h5file_id, 'spatial_extent', spatial_extent, 6)
+
+ end subroutine s_write_grid_to_hdf5
+
+ !> @brief Write a flow variable to HDF5 file
+ !!
+ !! Writes cell-centered flow variable data to the HDF5 file. The data
+ !! is stored in a dataset named after the variable. Respects the precision
+ !! setting (precision=1 for single, precision=2 for double).
+ !!
+ !! @param[in] varname Name of the variable (used as HDF5 dataset name)
+ !! @param[in] q_sf 3D array containing the variable data
+ subroutine s_write_variable_to_hdf5(varname, q_sf)
+
+ character(LEN=*), intent(in) :: varname
+ real(wp), dimension(-offset_x%beg:, -offset_y%beg:, -offset_z%beg:), intent(in) :: q_sf
+
+ integer(HSIZE_T), dimension(3) :: dims_3d
+ integer(HSIZE_T), dimension(2) :: dims_2d
+ integer(HSIZE_T), dimension(1) :: dims_1d
+
+ real(wp) :: data_range(2) ! min, max values for this variable
+
+ ! Determine dimensions based on problem dimensionality
+ if (p > 0) then
+ ! 3D case
+ dims_3d(1) = m + offset_x%beg + offset_x%end + 1
+ dims_3d(2) = n + offset_y%beg + offset_y%end + 1
+ dims_3d(3) = p + offset_z%beg + offset_z%end + 1
+ call h5screate_simple_f(3, dims_3d, h5dspace_id, h5err)
+ call h5dcreate_f(h5file_id, trim(varname), h5type_id, h5dspace_id, h5dset_id, h5err)
+ call h5dwrite_f(h5dset_id, h5type_id, q_sf, dims_3d, h5err)
+ ! Compute data range for this block
+ data_range(1) = minval(q_sf)
+ data_range(2) = maxval(q_sf)
+ elseif (n > 0) then
+ ! 2D case
+ dims_2d(1) = m + offset_x%beg + offset_x%end + 1
+ dims_2d(2) = n + offset_y%beg + offset_y%end + 1
+ call h5screate_simple_f(2, dims_2d, h5dspace_id, h5err)
+ call h5dcreate_f(h5file_id, trim(varname), h5type_id, h5dspace_id, h5dset_id, h5err)
+ call h5dwrite_f(h5dset_id, h5type_id, q_sf(:, :, 0), dims_2d, h5err)
+ ! Compute data range for this block
+ data_range(1) = minval(q_sf(:, :, 0))
+ data_range(2) = maxval(q_sf(:, :, 0))
+ else
+ ! 1D case
+ dims_1d(1) = m + offset_x%beg + offset_x%end + 1
+ call h5screate_simple_f(1, dims_1d, h5dspace_id, h5err)
+ call h5dcreate_f(h5file_id, trim(varname), h5type_id, h5dspace_id, h5dset_id, h5err)
+ call h5dwrite_f(h5dset_id, h5type_id, q_sf(:, 0, 0), dims_1d, h5err)
+ ! Compute data range for this block
+ data_range(1) = minval(q_sf(:, 0, 0))
+ data_range(2) = maxval(q_sf(:, 0, 0))
+ end if
+
+ ! Write data range as dataset attribute for visualization optimization
+ call s_write_attribute_real_array(h5dset_id, 'data_range', data_range, 2)
+
+ call h5dclose_f(h5dset_id, h5err)
+ call h5sclose_f(h5dspace_id, h5err)
+
+ ! Store variable name for XDMF generation
+ num_stored_vars = num_stored_vars + 1
+ stored_varnames(num_stored_vars) = trim(varname)
+
+ end subroutine s_write_variable_to_hdf5
+
+ !> @brief Write particle/point mesh coordinates to HDF5 file
+ !!
+ !! Creates datasets for Lagrangian particle positions (px, py, pz).
+ !! This data is used by XDMF to define a Polyvertex topology for
+ !! particle visualization. Uses h5type_id for precision-aware output.
+ !!
+ !! @param[in] px Array of particle x-coordinates
+ !! @param[in] py Array of particle y-coordinates
+ !! @param[in] pz Array of particle z-coordinates
+ !! @param[in] nBub Number of particles/bubbles on this processor
+ subroutine s_write_particle_mesh_to_hdf5(px, py, pz, nBub)
+
+ real(wp), dimension(:), intent(in) :: px, py, pz
+ integer, intent(in) :: nBub
+
+ integer(HSIZE_T), dimension(1) :: dims_1d
+ integer(HID_T) :: grp_id
+
+ ! Skip if no particles
+ if (nBub <= 0) then
+ stored_num_particles = 0
+ has_particle_data = .false.
+ return
+ end if
+
+ ! Create a group for particle data
+ call h5gcreate_f(h5file_id, 'particles', grp_id, h5err)
+
+ dims_1d(1) = nBub
+
+ ! Write X coordinates
+ call h5screate_simple_f(1, dims_1d, h5dspace_id, h5err)
+ call h5dcreate_f(grp_id, 'px', h5type_id, h5dspace_id, h5dset_id, h5err)
+ call h5dwrite_f(h5dset_id, h5type_id, px, dims_1d, h5err)
+ call h5dclose_f(h5dset_id, h5err)
+ call h5sclose_f(h5dspace_id, h5err)
+
+ ! Write Y coordinates
+ call h5screate_simple_f(1, dims_1d, h5dspace_id, h5err)
+ call h5dcreate_f(grp_id, 'py', h5type_id, h5dspace_id, h5dset_id, h5err)
+ call h5dwrite_f(h5dset_id, h5type_id, py, dims_1d, h5err)
+ call h5dclose_f(h5dset_id, h5err)
+ call h5sclose_f(h5dspace_id, h5err)
+
+ ! Write Z coordinates
+ call h5screate_simple_f(1, dims_1d, h5dspace_id, h5err)
+ call h5dcreate_f(grp_id, 'pz', h5type_id, h5dspace_id, h5dset_id, h5err)
+ call h5dwrite_f(h5dset_id, h5type_id, pz, dims_1d, h5err)
+ call h5dclose_f(h5dset_id, h5err)
+ call h5sclose_f(h5dspace_id, h5err)
+
+ call h5gclose_f(grp_id, h5err)
+
+ ! Store particle count for XDMF generation
+ stored_num_particles = nBub
+ has_particle_data = .true.
+ num_stored_particle_vars = 0
+
+ end subroutine s_write_particle_mesh_to_hdf5
+
+ !> @brief Write a particle variable to HDF5 file
+ !!
+ !! Writes scalar data associated with particles (e.g., radius, velocity
+ !! components, pressure) to the HDF5 file in the particles group.
+ !! Uses h5type_id for precision-aware output.
+ !!
+ !! @param[in] varname Name of the variable (used as HDF5 dataset name)
+ !! @param[in] data Array of variable values (one per particle)
+ !! @param[in] nBub Number of particles/bubbles
+ subroutine s_write_particle_variable_to_hdf5(varname, data, nBub)
+
+ character(LEN=*), intent(in) :: varname
+ real(wp), dimension(:), intent(in) :: data
+ integer, intent(in) :: nBub
+
+ integer(HSIZE_T), dimension(1) :: dims_1d
+ integer(HID_T) :: grp_id
+ character(LEN=name_len) :: dset_path
+
+ ! Skip if no particles
+ if (nBub <= 0) return
+
+ ! Open particles group
+ call h5gopen_f(h5file_id, 'particles', grp_id, h5err)
+
+ dims_1d(1) = nBub
+
+ ! Write variable data
+ call h5screate_simple_f(1, dims_1d, h5dspace_id, h5err)
+ call h5dcreate_f(grp_id, trim(varname), h5type_id, h5dspace_id, h5dset_id, h5err)
+ call h5dwrite_f(h5dset_id, h5type_id, data, dims_1d, h5err)
+ call h5dclose_f(h5dset_id, h5err)
+ call h5sclose_f(h5dspace_id, h5err)
+
+ call h5gclose_f(grp_id, h5err)
+
+ ! Store variable name for XDMF generation
+ num_stored_particle_vars = num_stored_particle_vars + 1
+ stored_particle_varnames(num_stored_particle_vars) = trim(varname)
+
+ end subroutine s_write_particle_variable_to_hdf5
+
+ !> @brief Close the current HDF5 file
+ subroutine s_close_hdf5_file()
+
+ call h5fclose_f(h5file_id, h5err)
+ if (h5err /= 0) then
+ call s_mpi_abort('Failed to close HDF5 file. Exiting.')
+ end if
+
+ end subroutine s_close_hdf5_file
+
+ !> @brief Write XDMF metadata file for a time step (root process only)
+ !!
+ !! Creates an XDMF file that describes the structure and location of data
+ !! in the HDF5 files. This file can be opened directly in ParaView.
+ !! The precision attribute in XDMF is set based on the global precision
+ !! setting: 4 bytes for single precision, 8 bytes for double.
+ !!
+ !! @param[in] t_step Current time step number
+ !! @param[in] time Physical simulation time
+ subroutine s_write_xdmf_file(t_step, time)
+
+ integer, intent(in) :: t_step
+ real(wp), intent(in) :: time
+
+ character(LEN=len_trim(case_dir) + 3*name_len) :: xdmf_file
+ character(LEN=256) :: dims_str, h5_path
+ character(LEN=1) :: prec_str !< Precision string for XDMF ("4" or "8")
+ integer :: unit_num, i, ip
+ integer :: nx, ny, nz
+
+ ! Only root process writes XDMF
+ if (proc_rank /= 0) return
+
+ ! Set precision string based on global precision setting
+ ! Grid coordinates always use double precision ("8") for accuracy
+ ! Variable data uses the user-specified precision
+ if (precision == 1) then
+ prec_str = '4' ! Single precision (32-bit float)
+ else
+ prec_str = '8' ! Double precision (64-bit float)
+ end if
+
+ ! Calculate grid dimensions (cell-centered values)
+ nx = m + offset_x%beg + offset_x%end + 1
+ if (n > 0) then
+ ny = n + offset_y%beg + offset_y%end + 1
+ else
+ ny = 1
+ end if
+ if (p > 0) then
+ nz = p + offset_z%beg + offset_z%end + 1
+ else
+ nz = 1
+ end if
+
+ ! Generate XDMF file path
+ write (xdmf_file, '(A,I0,A)') '/timestep_', t_step, '.xdmf'
+ xdmf_file = trim(rootdir_hdf5)//trim(xdmf_file)
+
+ unit_num = 100
+ open (unit=unit_num, file=trim(xdmf_file), status='replace', action='write')
+
+ ! Write XDMF header
+ write (unit_num, '(A)') ''
+ write (unit_num, '(A)') ''
+ write (unit_num, '(A)') ''
+ write (unit_num, '(A)') ' '
+
+ ! Write grid collection for all processors
+ write (unit_num, '(A)') ' '
+ write (unit_num, '(A,E15.8,A)') ' '
+
+ ! Write grid for each processor
+ do ip = 0, num_procs - 1
+ write (h5_path, '(A,I0,A,I0,A)') '../p', ip, '/', t_step, '.h5'
+
+ write (unit_num, '(A,I0,A)') ' '
+
+ ! Topology (structured grid)
+ if (p > 0) then
+ write (dims_str, '(I0,A,I0,A,I0)') nz + 1, ' ', ny + 1, ' ', nx + 1
+ write (unit_num, '(A,A,A)') ' '
+ elseif (n > 0) then
+ write (dims_str, '(I0,A,I0)') ny + 1, ' ', nx + 1
+ write (unit_num, '(A,A,A)') ' '
+ else
+ write (dims_str, '(I0)') nx + 1
+ write (unit_num, '(A,A,A)') ' '
+ end if
+
+ ! Geometry (rectilinear mesh coordinates)
+ if (p > 0) then
+ write (unit_num, '(A)') ' '
+ ! X coordinates
+ write (unit_num, '(A,I0,A)') ' '
+ write (unit_num, '(A,A,A)') ' ', trim(h5_path), ':/x_cb'
+ write (unit_num, '(A)') ' '
+ ! Y coordinates
+ write (unit_num, '(A,I0,A)') ' '
+ write (unit_num, '(A,A,A)') ' ', trim(h5_path), ':/y_cb'
+ write (unit_num, '(A)') ' '
+ ! Z coordinates
+ write (unit_num, '(A,I0,A)') ' '
+ write (unit_num, '(A,A,A)') ' ', trim(h5_path), ':/z_cb'
+ write (unit_num, '(A)') ' '
+ write (unit_num, '(A)') ' '
+ elseif (n > 0) then
+ write (unit_num, '(A)') ' '
+ ! X coordinates
+ write (unit_num, '(A,I0,A)') ' '
+ write (unit_num, '(A,A,A)') ' ', trim(h5_path), ':/x_cb'
+ write (unit_num, '(A)') ' '
+ ! Y coordinates
+ write (unit_num, '(A,I0,A)') ' '
+ write (unit_num, '(A,A,A)') ' ', trim(h5_path), ':/y_cb'
+ write (unit_num, '(A)') ' '
+ write (unit_num, '(A)') ' '
+ else
+ write (unit_num, '(A)') ' '
+ ! X coordinates
+ write (unit_num, '(A,I0,A)') ' '
+ write (unit_num, '(A,A,A)') ' ', trim(h5_path), ':/x_cb'
+ write (unit_num, '(A)') ' '
+ ! Y = 0 for 1D (single value)
+ write (unit_num, '(A)') ' '
+ write (unit_num, '(A)') ' 0.0'
+ write (unit_num, '(A)') ' '
+ write (unit_num, '(A)') ' '
+ end if
+
+ ! Attributes (flow variables) - use precision-aware output
+ do i = 1, num_stored_vars
+ write (unit_num, '(A,A,A)') ' '
+ if (p > 0) then
+ write (dims_str, '(I0,A,I0,A,I0)') nz, ' ', ny, ' ', nx
+ elseif (n > 0) then
+ write (dims_str, '(I0,A,I0)') ny, ' ', nx
+ else
+ write (dims_str, '(I0)') nx
+ end if
+ write (unit_num, '(A,A,A,A,A)') ' '
+ write (unit_num, '(A,A,A,A)') ' ', trim(h5_path), ':/', trim(stored_varnames(i))
+ write (unit_num, '(A)') ' '
+ write (unit_num, '(A)') ' '
+ end do
+
+ write (unit_num, '(A)') ' '
+ end do
+
+ ! Write particle data grids if particle data exists
+ if (has_particle_data .and. stored_num_particles > 0) then
+ do ip = 0, num_procs - 1
+ write (h5_path, '(A,I0,A,I0,A)') '../p', ip, '/', t_step, '.h5'
+
+ write (unit_num, '(A,I0,A)') ' '
+
+ ! Topology (polyvertex for points)
+ write (unit_num, '(A,I0,A)') ' '
+
+ ! Geometry (XYZ coordinates) - use precision-aware output
+ write (unit_num, '(A)') ' '
+ write (unit_num, '(A,I0,A,A,A)') ' '
+
+ ! Join px, py, pz into XYZ using DataItem function
+ write (unit_num, '(A)') ' '
+ write (unit_num, '(A,I0,A,A,A)') ' '
+ write (unit_num, '(A,A,A)') ' ', trim(h5_path), ':/particles/px'
+ write (unit_num, '(A)') ' '
+ write (unit_num, '(A,I0,A,A,A)') ' '
+ write (unit_num, '(A,A,A)') ' ', trim(h5_path), ':/particles/py'
+ write (unit_num, '(A)') ' '
+ write (unit_num, '(A,I0,A,A,A)') ' '
+ write (unit_num, '(A,A,A)') ' ', trim(h5_path), ':/particles/pz'
+ write (unit_num, '(A)') ' '
+ write (unit_num, '(A)') ' '
+
+ write (unit_num, '(A)') ' '
+ write (unit_num, '(A)') ' '
+
+ ! Particle attributes - use precision-aware output
+ do i = 1, num_stored_particle_vars
+ write (unit_num, '(A,A,A)') ' '
+ write (unit_num, '(A,I0,A,A,A)') ' '
+ write (unit_num, '(A,A,A,A)') ' ', trim(h5_path), ':/particles/', &
+ trim(stored_particle_varnames(i))
+ write (unit_num, '(A)') ' '
+ write (unit_num, '(A)') ' '
+ end do
+
+ write (unit_num, '(A)') ' '
+ end do
+ end if
+
+ write (unit_num, '(A)') ' '
+ write (unit_num, '(A)') ' '
+ write (unit_num, '(A)') ''
+
+ close (unit_num)
+
+ end subroutine s_write_xdmf_file
+
+ !> @brief Write master XDMF time series file (root process only)
+ !!
+ !! Creates a master XDMF file that references all individual timestep
+ !! XDMF files using XI:Include. This allows ParaView to load all timesteps
+ !! as a single time-varying dataset using the temporal collection.
+ !!
+ !! @param[in] t_step_start First time step to include
+ !! @param[in] t_step_stop Last time step to include
+ !! @param[in] t_step_save Time step interval between saves
+ subroutine s_write_xdmf_time_series(t_step_start, t_step_stop, t_step_save)
+
+ integer, intent(in) :: t_step_start
+ integer, intent(in) :: t_step_stop
+ integer, intent(in) :: t_step_save
+
+ character(LEN=len_trim(case_dir) + 3*name_len) :: xdmf_file
+ integer :: unit_num, t_step
+
+ ! Only root process writes XDMF
+ if (proc_rank /= 0) return
+
+ ! Generate time series XDMF file path
+ xdmf_file = trim(rootdir_hdf5)//'/time_series.xdmf'
+
+ unit_num = 101
+ open (unit=unit_num, file=trim(xdmf_file), status='replace', action='write')
+
+ ! Write XDMF header with XInclude namespace
+ write (unit_num, '(A)') ''
+ write (unit_num, '(A)') ''
+ write (unit_num, '(A)') ''
+ write (unit_num, '(A)') ' '
+ write (unit_num, '(A)') ' '
+
+ ! Include each timestep XDMF file
+ do t_step = t_step_start, t_step_stop, t_step_save
+ write (unit_num, '(A,I0,A)') ' '
+ end do
+
+ write (unit_num, '(A)') ' '
+ write (unit_num, '(A)') ' '
+ write (unit_num, '(A)') '