diff --git a/Makefile b/Makefile
index 82e4dca..b899e11 100644
--- a/Makefile
+++ b/Makefile
@@ -20,6 +20,8 @@ LDFLAGS ?= # user can set default linker flags
 FFLAGS = $(FCFLAGS)
 LFLAGS = $(LDFLAGS)
 
+NTHREADS ?= 1 # gfortran requires number of threads specified at compile time...
+
 #######CMP settings###########
 ifeq ($(CMP),intel)
   FC = mpiifort
@@ -38,7 +40,7 @@ else ifeq ($(CMP),gcc)
     FFLAGS += -ffpe-trap=invalid,zero -fcheck=all -fimplicit-none
   else
     FFLAGS += -O3 -march=native
-    FFLAGS += -fopenmp -ftree-parallelize-loops=12
+    FFLAGS += -fopenmp -ftree-parallelize-loops=$(NTHREADS)
     LFLAGS += -fopenmp
   endif
 else ifeq ($(CMP),nagfor)
@@ -59,10 +61,10 @@ DECOMPDIR = ./decomp2d
 SRCDIR = ./src
 
 ### List of files for the main code
-SRCDECOMP = $(DECOMPDIR)/decomp_2d.f90 $(DECOMPDIR)/glassman.f90 $(DECOMPDIR)/fft_$(FFT).f90 
+SRCDECOMP = $(DECOMPDIR)/decomp_2d.f90 $(DECOMPDIR)/glassman.f90 $(DECOMPDIR)/fft_$(FFT).f90 $(DECOMPDIR)/io.f90 
 OBJDECOMP = $(SRCDECOMP:%.f90=%.o)
 OBJ = $(SRC:%.f90=%.o)
-SRC = $(SRCDIR)/x3d_precision.f90 $(SRCDIR)/module_param.f90 $(SRCDIR)/time_integrators.f90 $(SRCDIR)/tools.f90 $(SRCDIR)/x3d_transpose.f90 $(SRCDIR)/var.f90 $(SRCDIR)/thomas.f90 $(SRCDIR)/x3d_operator_x_data.f90 $(SRCDIR)/x3d_operator_y_data.f90 $(SRCDIR)/x3d_operator_z_data.f90 $(SRCDIR)/x3d_operator_1d.f90 $(SRCDIR)/poisson.f90 $(SRCDIR)/x3d_derive.f90 $(SRCDIR)/x3d_staggered.f90 $(SRCDIR)/x3d_filters.f90 $(SRCDIR)/bc_tgv.f90 $(SRCDIR)/navier.f90 $(SRCDIR)/parameters.f90 $(SRCDIR)/bc_tgv2d.f90 $(SRCDIR)/case.f90 $(SRCDIR)/transeq.f90 $(SRCDIR)/x3d_tools.f90 $(SRCDIR)/xcompact3d.f90
+SRC = $(SRCDIR)/x3d_precision.f90 $(SRCDIR)/module_param.f90 $(SRCDIR)/time_integrators.f90 $(SRCDIR)/tools.f90 $(SRCDIR)/x3d_transpose.f90 $(SRCDIR)/var.f90 $(SRCDIR)/thomas.f90 $(SRCDIR)/x3d_operator_x_data.f90 $(SRCDIR)/x3d_operator_y_data.f90 $(SRCDIR)/x3d_operator_z_data.f90 $(SRCDIR)/x3d_operator_1d.f90 $(SRCDIR)/poisson.f90 $(SRCDIR)/x3d_derive.f90 $(SRCDIR)/x3d_staggered.f90 $(SRCDIR)/x3d_filters.f90 $(SRCDIR)/visu.f90 $(SRCDIR)/bc_tgv.f90 $(SRCDIR)/navier.f90 $(SRCDIR)/parameters.f90 $(SRCDIR)/bc_tgv2d.f90 $(SRCDIR)/case.f90 $(SRCDIR)/transeq.f90 $(SRCDIR)/x3d_tools.f90 $(SRCDIR)/xcompact3d.f90
 
 #######FFT settings##########
 ifeq ($(FFT),fftw3)
diff --git a/decomp2d/CMakeLists.txt b/decomp2d/CMakeLists.txt
index 8ba49d1..44bb5c8 100644
--- a/decomp2d/CMakeLists.txt
+++ b/decomp2d/CMakeLists.txt
@@ -1,4 +1,5 @@
 file(GLOB files_decomp decomp_2d.f90
+                       io.f90
                        glassman.f90)
 include_directories(${CMAKE_SOURCE_DIR}/decomp2d)
 
diff --git a/decomp2d/Makefile.am b/decomp2d/Makefile.am
index 2530f42..7983b7b 100644
--- a/decomp2d/Makefile.am
+++ b/decomp2d/Makefile.am
@@ -9,6 +9,7 @@ endif
 lib_LTLIBRARIES = libdecomp2d.la
 libdecomp2d_la_SOURCES = \
 decomp_2d.f90 \
+io.f90 \
 glassman.f90 \
 fft_@FFTENGINE@.f90
 libdecomp2d_la_LDFLAGS = -no-undefined
@@ -18,6 +19,7 @@ endif
 
 # Allow parallel make
 glassman.lo: decomp_2d.lo
+io.lo: decomp_2d.lo
 if HAVE_MKL
 fft_@FFTENGINE@.lo: glassman.lo mkl_dfti.lo
 else
diff --git a/decomp2d/io.f90 b/decomp2d/io.f90
new file mode 100644
index 0000000..ae14ad9
--- /dev/null
+++ b/decomp2d/io.f90
@@ -0,0 +1,1198 @@
+!=======================================================================
+! This is part of the 2DECOMP&FFT library
+! 
+! 2DECOMP&FFT is a software framework for general-purpose 2D (pencil) 
+! decomposition. It also implements a highly scalable distributed
+! three-dimensional Fast Fourier Transform (FFT).
+!
+! Copyright (C) 2009-2013 Ning Li, the Numerical Algorithms Group (NAG)
+! Copyright (C) 2021               the University of Edinburgh (UoE)
+!
+!=======================================================================
+
+! This module provides parallel IO facilities for applications based on
+! 2D decomposition.
+
+module decomp_2d_io
+
+  use decomp_2d
+  use MPI
+#ifdef T3PIO
+  use t3pio
+#endif
+
+#ifdef ADIOS2
+  use adios2
+#endif
+
+  implicit none
+
+  private        ! Make everything private unless declared public
+
+  public :: decomp_2d_write_one, decomp_2d_read_one, &
+       decomp_2d_write_var, decomp_2d_read_var, &
+       decomp_2d_write_scalar, decomp_2d_read_scalar, &
+       decomp_2d_write_plane, decomp_2d_write_every, &
+       decomp_2d_write_subdomain, &
+       decomp_2d_write_outflow, decomp_2d_read_inflow
+#ifdef ADIOS2
+  public :: adios2_register_variable
+#endif
+
+  ! Generic interface to handle multiple data types
+
+  interface decomp_2d_write_one
+     module procedure write_one_real
+     module procedure write_one_complex
+     module procedure mpiio_write_real_coarse
+     module procedure mpiio_write_real_probe
+#ifdef ADIOS2
+     module procedure adios2_write_real_coarse
+#endif
+  end interface decomp_2d_write_one
+
+  interface decomp_2d_read_one
+     module procedure read_one_real
+     module procedure read_one_complex
+  end interface decomp_2d_read_one
+
+  interface decomp_2d_write_var
+     module procedure write_var_real
+     module procedure write_var_complex
+  end interface decomp_2d_write_var
+
+  interface decomp_2d_read_var
+     module procedure read_var_real
+     module procedure read_var_complex
+  end interface decomp_2d_read_var
+
+  interface decomp_2d_write_scalar
+     module procedure write_scalar_real
+     module procedure write_scalar_complex
+     module procedure write_scalar_integer
+     module procedure write_scalar_logical
+  end interface decomp_2d_write_scalar
+
+  interface decomp_2d_read_scalar
+     module procedure read_scalar_real
+     module procedure read_scalar_complex
+     module procedure read_scalar_integer
+     module procedure read_scalar_logical
+  end interface decomp_2d_read_scalar
+
+  interface decomp_2d_write_plane
+     module procedure write_plane_3d_real
+     module procedure write_plane_3d_complex
+     !     module procedure write_plane_2d
+  end interface decomp_2d_write_plane
+
+  interface decomp_2d_write_every
+     module procedure write_every_real
+     module procedure write_every_complex
+  end interface decomp_2d_write_every
+
+  interface decomp_2d_write_subdomain
+     module procedure write_subdomain
+  end interface decomp_2d_write_subdomain
+
+  interface decomp_2d_write_outflow
+     module procedure write_outflow
+  end interface decomp_2d_write_outflow
+
+  interface decomp_2d_read_inflow
+     module procedure read_inflow
+  end interface decomp_2d_read_inflow
+
+contains
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  ! Using MPI-IO library to write a single 3D array to a file
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  subroutine write_one_real(ipencil,var,filename,opt_decomp)
+
+    implicit none
+
+    integer, intent(IN) :: ipencil
+    real(mytype), dimension(:,:,:), intent(IN) :: var
+    character(len=*), intent(IN) :: filename
+    TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp
+
+    TYPE(DECOMP_INFO) :: decomp
+    integer(kind=MPI_OFFSET_KIND) :: filesize, disp
+    integer, dimension(3) :: sizes, subsizes, starts
+    integer :: ierror, newtype, fh, data_type, info, gs
+
+    data_type = real_type
+
+#include "io_write_one.inc"
+
+    return
+  end subroutine write_one_real
+
+
+  subroutine write_one_complex(ipencil,var,filename,opt_decomp)
+
+    implicit none
+
+    integer, intent(IN) :: ipencil
+    complex(mytype), dimension(:,:,:), intent(IN) :: var
+    character(len=*), intent(IN) :: filename
+    TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp
+
+    TYPE(DECOMP_INFO) :: decomp
+    integer(kind=MPI_OFFSET_KIND) :: filesize, disp
+    integer, dimension(3) :: sizes, subsizes, starts
+    integer :: ierror, newtype, fh, data_type, info, gs
+
+    data_type = complex_type
+
+#include "io_write_one.inc"
+
+    return
+  end subroutine write_one_complex
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  ! Using MPI-IO library to read from a file a single 3D array
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  subroutine read_one_real(ipencil,var,filename,opt_decomp)
+
+    implicit none
+
+    integer, intent(IN) :: ipencil
+    real(mytype), dimension(:,:,:), intent(INOUT) :: var
+    character(len=*), intent(IN) :: filename
+    TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp
+
+    TYPE(DECOMP_INFO) :: decomp
+    integer(kind=MPI_OFFSET_KIND) :: disp
+    integer, dimension(3) :: sizes, subsizes, starts
+    integer :: ierror, newtype, fh, data_type
+    real(mytype_single), allocatable, dimension(:,:,:) :: varsingle
+
+    data_type = real_type_single
+    allocate (varsingle(xstV(1):xenV(1),xstV(2):xenV(2),xstV(3):xenV(3)))
+
+    if (present(opt_decomp)) then
+       decomp = opt_decomp
+    else
+       call get_decomp_info(decomp)
+    end if
+
+    ! determine subarray parameters
+    sizes(1) = decomp%xsz(1)
+    sizes(2) = decomp%ysz(2)
+    sizes(3) = decomp%zsz(3)
+
+    if (ipencil == 1) then
+       subsizes(1) = decomp%xsz(1)
+       subsizes(2) = decomp%xsz(2)
+       subsizes(3) = decomp%xsz(3)
+       starts(1) = decomp%xst(1)-1  ! 0-based index
+       starts(2) = decomp%xst(2)-1
+       starts(3) = decomp%xst(3)-1
+    else if (ipencil == 2) then
+       subsizes(1) = decomp%ysz(1)
+       subsizes(2) = decomp%ysz(2)
+       subsizes(3) = decomp%ysz(3)
+       starts(1) = decomp%yst(1)-1
+       starts(2) = decomp%yst(2)-1
+       starts(3) = decomp%yst(3)-1
+    else if (ipencil == 3) then
+       subsizes(1) = decomp%zsz(1)
+       subsizes(2) = decomp%zsz(2)
+       subsizes(3) = decomp%zsz(3)
+       starts(1) = decomp%zst(1)-1
+       starts(2) = decomp%zst(2)-1
+       starts(3) = decomp%zst(3)-1
+    endif
+
+    call MPI_TYPE_CREATE_SUBARRAY(3, sizes, subsizes, starts,  &
+         MPI_ORDER_FORTRAN, data_type, newtype, ierror)
+    call MPI_TYPE_COMMIT(newtype,ierror)
+    call MPI_FILE_OPEN(MPI_COMM_WORLD, filename, &
+         MPI_MODE_RDONLY, MPI_INFO_NULL, &
+         fh, ierror)
+    disp = 0_MPI_OFFSET_KIND
+    call MPI_FILE_SET_VIEW(fh,disp,data_type, &
+         newtype,'native',MPI_INFO_NULL,ierror)
+    call MPI_FILE_READ_ALL(fh, varsingle, &
+         subsizes(1)*subsizes(2)*subsizes(3), &
+         data_type, MPI_STATUS_IGNORE, ierror)
+    call MPI_FILE_CLOSE(fh,ierror)
+    call MPI_TYPE_FREE(newtype,ierror)
+    var = real(varsingle,mytype)
+    deallocate(varsingle)
+
+    return
+  end subroutine read_one_real
+
+
+  subroutine read_one_complex(ipencil,var,filename,opt_decomp)
+
+    implicit none
+
+    integer, intent(IN) :: ipencil
+    complex(mytype), dimension(:,:,:), intent(INOUT) :: var
+    character(len=*), intent(IN) :: filename
+    TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp
+
+    TYPE(DECOMP_INFO) :: decomp
+    integer(kind=MPI_OFFSET_KIND) :: disp
+    integer, dimension(3) :: sizes, subsizes, starts
+    integer :: ierror, newtype, fh, data_type
+
+    data_type = complex_type
+
+#include "io_read_one.inc"
+
+    return
+  end subroutine read_one_complex
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  ! Write a 3D array as part of a big MPI-IO file, starting from 
+  !  displacement 'disp'; 'disp' will be updated after the writing
+  !  operation to prepare the writing of next chunk of data.
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  subroutine write_var_real(fh,disp,ipencil,var,opt_decomp)
+
+    implicit none
+
+    integer, intent(IN) :: fh
+    integer(KIND=MPI_OFFSET_KIND), intent(INOUT) :: disp
+    integer, intent(IN) :: ipencil
+    real(mytype), dimension(:,:,:), intent(IN) :: var
+    TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp
+
+    TYPE(DECOMP_INFO) :: decomp
+    integer, dimension(3) :: sizes, subsizes, starts
+    integer :: ierror, newtype, data_type
+
+    data_type = real_type
+
+#include "io_write_var.inc"
+
+    return
+  end subroutine write_var_real
+
+
+  subroutine write_var_complex(fh,disp,ipencil,var,opt_decomp)
+
+    implicit none
+
+    integer, intent(IN) :: fh
+    integer(KIND=MPI_OFFSET_KIND), intent(INOUT) :: disp
+    integer, intent(IN) :: ipencil
+    complex(mytype), dimension(:,:,:), intent(IN) :: var
+    TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp
+
+    TYPE(DECOMP_INFO) :: decomp
+    integer, dimension(3) :: sizes, subsizes, starts
+    integer :: ierror, newtype, data_type
+
+    data_type = complex_type
+
+#include "io_write_var.inc"
+
+    return
+  end subroutine write_var_complex
+
+
+  subroutine write_outflow(fh,disp,ntimesteps,var,opt_decomp)
+
+    implicit none
+
+    integer, intent(IN) :: fh
+    integer(KIND=MPI_OFFSET_KIND), intent(INOUT) :: disp
+    integer, intent(IN) :: ntimesteps 
+    real(mytype), dimension(:,:,:), intent(IN) :: var
+    TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp
+
+    TYPE(DECOMP_INFO) :: decomp
+    integer, dimension(3) :: sizes, subsizes, starts
+    integer :: ierror, newtype, data_type
+
+    data_type = real_type
+
+#include "io_write_outflow.f90"
+
+    return
+  end subroutine write_outflow
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  ! Read a 3D array as part of a big MPI-IO file, starting from 
+  !  displacement 'disp'; 'disp' will be updated after the reading
+  !  operation to prepare the reading of next chunk of data.
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  subroutine read_var_real(fh,disp,ipencil,var,opt_decomp)
+
+    implicit none
+
+    integer, intent(IN) :: fh
+    integer(KIND=MPI_OFFSET_KIND), intent(INOUT) :: disp
+    integer, intent(IN) :: ipencil
+    real(mytype), dimension(:,:,:), intent(INOUT) :: var
+    TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp
+
+    TYPE(DECOMP_INFO) :: decomp
+    integer, dimension(3) :: sizes, subsizes, starts
+    integer :: ierror, newtype, data_type
+
+    data_type = real_type
+
+#include "io_read_var.inc"
+
+    return
+  end subroutine read_var_real
+
+
+  subroutine read_var_complex(fh,disp,ipencil,var,opt_decomp)
+
+    implicit none
+
+    integer, intent(IN) :: fh
+    integer(KIND=MPI_OFFSET_KIND), intent(INOUT) :: disp
+    integer, intent(IN) :: ipencil
+    complex(mytype), dimension(:,:,:), intent(INOUT) :: var
+    TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp
+
+    TYPE(DECOMP_INFO) :: decomp
+    integer, dimension(3) :: sizes, subsizes, starts
+    integer :: ierror, newtype, data_type
+
+    data_type = complex_type
+
+#include "io_read_var.inc"
+
+    return
+  end subroutine read_var_complex
+
+
+  subroutine read_inflow(fh,disp,ntimesteps,var,opt_decomp)
+
+    implicit none
+
+    integer, intent(IN) :: fh
+    integer(KIND=MPI_OFFSET_KIND), intent(INOUT) :: disp
+    integer, intent(IN) :: ntimesteps
+    real(mytype), dimension(:,:,:), intent(INOUT) :: var
+    TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp
+
+    TYPE(DECOMP_INFO) :: decomp
+    integer, dimension(3) :: sizes, subsizes, starts
+    integer :: ierror, newtype, data_type
+
+    data_type = real_type
+
+#include "io_read_inflow.f90"
+
+    return
+  end subroutine read_inflow
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  ! Write scalar variables as part of a big MPI-IO file, starting from 
+  !  displacement 'disp'; 'disp' will be updated after the reading
+  !  operation to prepare the reading of next chunk of data.
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  subroutine write_scalar_real(fh,disp,n,var)
+
+    implicit none
+
+    integer, intent(IN) :: fh             ! file handle
+    integer(KIND=MPI_OFFSET_KIND), &
+         intent(INOUT) :: disp            ! displacement
+    integer, intent(IN) :: n              ! number of scalars
+    real(mytype), dimension(n), &
+         intent(IN) :: var                ! array of scalars
+
+    integer :: m, ierror
+
+    call MPI_FILE_SET_VIEW(fh,disp,real_type, &
+         real_type,'native',MPI_INFO_NULL,ierror)
+    if (nrank==0) then
+       m = n ! only one rank needs to write
+    else
+       m = 0
+    end if
+    call MPI_FILE_WRITE_ALL(fh, var, m, real_type, &
+         MPI_STATUS_IGNORE, ierror)
+    disp = disp + n*mytype_bytes
+
+    return
+  end subroutine write_scalar_real
+
+
+  subroutine write_scalar_complex(fh,disp,n,var)
+
+    implicit none
+
+    integer, intent(IN) :: fh
+    integer(KIND=MPI_OFFSET_KIND), intent(INOUT) :: disp
+    integer, intent(IN) :: n
+    complex(mytype), dimension(n), intent(IN) :: var
+
+    integer :: m, ierror
+
+    call MPI_FILE_SET_VIEW(fh,disp,complex_type, &
+         complex_type,'native',MPI_INFO_NULL,ierror)
+    if (nrank==0) then
+       m = n
+    else
+       m = 0
+    end if
+    call MPI_FILE_WRITE_ALL(fh, var, m, complex_type, &
+         MPI_STATUS_IGNORE, ierror)
+    disp = disp + n*mytype_bytes*2
+
+    return
+  end subroutine write_scalar_complex
+
+
+  subroutine write_scalar_integer(fh,disp,n,var)
+
+    implicit none
+
+    integer, intent(IN) :: fh
+    integer(KIND=MPI_OFFSET_KIND), intent(INOUT) :: disp
+    integer, intent(IN) :: n
+    integer, dimension(n), intent(IN) :: var
+
+    integer :: m, ierror
+
+    call MPI_FILE_SET_VIEW(fh,disp,MPI_INTEGER, &
+         MPI_INTEGER,'native',MPI_INFO_NULL,ierror)
+    if (nrank==0) then
+       m = n
+    else
+       m = 0
+    end if
+    call MPI_FILE_WRITE_ALL(fh, var, m, MPI_INTEGER, &
+         MPI_STATUS_IGNORE, ierror)
+    call MPI_TYPE_SIZE(MPI_INTEGER,m,ierror)
+    disp = disp + n*m
+
+    return
+  end subroutine write_scalar_integer
+
+
+  subroutine write_scalar_logical(fh,disp,n,var)
+
+    implicit none
+
+    integer, intent(IN) :: fh
+    integer(KIND=MPI_OFFSET_KIND), intent(INOUT) :: disp
+    integer, intent(IN) :: n
+    logical, dimension(n), intent(IN) :: var
+
+    integer :: m, ierror
+
+    call MPI_FILE_SET_VIEW(fh,disp,MPI_LOGICAL, &
+         MPI_LOGICAL,'native',MPI_INFO_NULL,ierror)
+    if (nrank==0) then
+       m = n
+    else
+       m = 0
+    end if
+    call MPI_FILE_WRITE_ALL(fh, var, m, MPI_LOGICAL, &
+         MPI_STATUS_IGNORE, ierror)
+    call MPI_TYPE_SIZE(MPI_LOGICAL,m,ierror)
+    disp = disp + n*m
+
+    return
+  end subroutine write_scalar_logical
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  ! Read scalar variables as part of a big MPI-IO file, starting from 
+  !  displacement 'disp'; 'disp' will be updated after the reading
+  !  operation to prepare the reading of next chunk of data.
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  subroutine read_scalar_real(fh,disp,n,var)
+
+    implicit none
+
+    integer, intent(IN) :: fh             ! file handle
+    integer(KIND=MPI_OFFSET_KIND), &
+         intent(INOUT) :: disp            ! displacement
+    integer, intent(IN) :: n              ! number of scalars
+    real(mytype), dimension(n), &
+         intent(INOUT) :: var             ! array of scalars
+
+    integer :: ierror
+
+    call MPI_FILE_SET_VIEW(fh,disp,real_type, &
+         real_type,'native',MPI_INFO_NULL,ierror)
+    call MPI_FILE_READ_ALL(fh, var, n, real_type, &
+         MPI_STATUS_IGNORE, ierror)
+    disp = disp + n*mytype_bytes
+
+    return
+  end subroutine read_scalar_real
+
+
+  subroutine read_scalar_complex(fh,disp,n,var)
+
+    implicit none
+
+    integer, intent(IN) :: fh
+    integer(KIND=MPI_OFFSET_KIND), intent(INOUT) :: disp
+    integer, intent(IN) :: n
+    complex(mytype), dimension(n), intent(INOUT) :: var
+
+    integer :: ierror
+
+    call MPI_FILE_SET_VIEW(fh,disp,complex_type, &
+         complex_type,'native',MPI_INFO_NULL,ierror)
+    call MPI_FILE_READ_ALL(fh, var, n, complex_type, &
+         MPI_STATUS_IGNORE, ierror)
+    disp = disp + n*mytype_bytes*2
+
+    return
+  end subroutine read_scalar_complex
+
+
+  subroutine read_scalar_integer(fh,disp,n,var)
+
+    implicit none
+
+    integer, intent(IN) :: fh
+    integer(KIND=MPI_OFFSET_KIND), intent(INOUT) :: disp
+    integer, intent(IN) :: n
+    integer, dimension(n), intent(INOUT) :: var
+
+    integer :: m, ierror
+
+    call MPI_FILE_SET_VIEW(fh,disp,MPI_INTEGER, &
+         MPI_INTEGER,'native',MPI_INFO_NULL,ierror)
+    call MPI_FILE_READ_ALL(fh, var, n, MPI_INTEGER, &
+         MPI_STATUS_IGNORE, ierror)
+    call MPI_TYPE_SIZE(MPI_INTEGER,m,ierror)
+    disp = disp + n*m
+
+    return
+  end subroutine read_scalar_integer
+
+
+  subroutine read_scalar_logical(fh,disp,n,var)
+
+    implicit none
+
+    integer, intent(IN) :: fh
+    integer(KIND=MPI_OFFSET_KIND), intent(INOUT) :: disp
+    integer, intent(IN) :: n
+    logical, dimension(n), intent(INOUT) :: var
+
+    integer :: m, ierror
+
+    call MPI_FILE_SET_VIEW(fh,disp,MPI_LOGICAL, &
+         MPI_LOGICAL,'native',MPI_INFO_NULL,ierror)
+    call MPI_FILE_READ_ALL(fh, var, n, MPI_LOGICAL, &
+         MPI_STATUS_IGNORE, ierror)
+    call MPI_TYPE_SIZE(MPI_LOGICAL,m,ierror)
+    disp = disp + n*m
+
+    return
+  end subroutine read_scalar_logical
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  ! Write a 2D slice of the 3D data to a file
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  subroutine write_plane_3d_real(ipencil,var,iplane,n,filename, &
+       opt_decomp)
+
+    implicit none
+
+    integer, intent(IN) :: ipencil !(x-pencil=1; y-pencil=2; z-pencil=3)
+    real(mytype), dimension(:,:,:), intent(IN) :: var
+    integer, intent(IN) :: iplane !(x-plane=1; y-plane=2; z-plane=3)
+    integer, intent(IN) :: n ! which plane to write (global coordinate)
+    character(len=*), intent(IN) :: filename
+    TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp
+
+    real(mytype), allocatable, dimension(:,:,:) :: wk, wk2
+    real(mytype), allocatable, dimension(:,:,:) :: wk2d
+    TYPE(DECOMP_INFO) :: decomp
+    integer(kind=MPI_OFFSET_KIND) :: filesize, disp
+    integer, dimension(3) :: sizes, subsizes, starts
+    integer :: i,j,k, ierror, newtype, fh, data_type
+
+    data_type = real_type
+
+#include "io_write_plane.inc"
+
+    return
+  end subroutine write_plane_3d_real
+
+
+  subroutine write_plane_3d_complex(ipencil,var,iplane,n, &
+       filename,opt_decomp)
+
+    implicit none
+
+    integer, intent(IN) :: ipencil !(x-pencil=1; y-pencil=2; z-pencil=3)
+    complex(mytype), dimension(:,:,:), intent(IN) :: var
+    integer, intent(IN) :: iplane !(x-plane=1; y-plane=2; z-plane=3)
+    integer, intent(IN) :: n ! which plane to write (global coordinate)
+    character(len=*), intent(IN) :: filename
+    TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp
+
+    complex(mytype), allocatable, dimension(:,:,:) :: wk, wk2
+    complex(mytype), allocatable, dimension(:,:,:) :: wk2d
+    TYPE(DECOMP_INFO) :: decomp
+    integer(kind=MPI_OFFSET_KIND) :: filesize, disp
+    integer, dimension(3) :: sizes, subsizes, starts
+    integer :: i,j,k, ierror, newtype, fh, data_type
+
+    data_type = complex_type
+
+#include "io_write_plane.inc"
+
+    return
+  end subroutine write_plane_3d_complex
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  ! Write a 2D array to a file
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !************** TO DO ***************
+  !* Consider handling distributed 2D data set
+  !  subroutine write_plane_2d(ipencil,var,filename)
+  !    integer, intent(IN) :: ipencil !(x-pencil=1; y-pencil=2; z-pencil=3)
+  !    real(mytype), dimension(:,:), intent(IN) :: var ! 2D array
+  !    character(len=*), intent(IN) :: filename
+  !
+  !    if (ipencil==1) then
+  !       ! var should be defined as var(xsize(2)
+  !
+  !    else if (ipencil==2) then
+  !
+  !    else if (ipencil==3) then
+  !
+  !    end if
+  !
+  !    return
+  !  end subroutine write_plane_2d
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  ! Write 3D array data for every specified mesh point
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  subroutine write_every_real(ipencil,var,iskip,jskip,kskip, &
+       filename, from1)
+
+    implicit none
+
+    integer, intent(IN) :: ipencil !(x-pencil=1; y-pencil=2; z-pencil=3)
+    real(mytype), dimension(:,:,:), intent(IN) :: var
+    integer, intent(IN) :: iskip,jskip,kskip 
+    character(len=*), intent(IN) :: filename
+    logical, intent(IN) :: from1  ! .true.  - save 1,n+1,2n+1...
+    ! .false. - save n,2n,3n...
+
+    real(mytype), allocatable, dimension(:,:,:) :: wk, wk2
+    integer(kind=MPI_OFFSET_KIND) :: filesize, disp
+    integer, dimension(3) :: sizes, subsizes, starts
+    integer :: i,j,k, ierror, newtype, fh, key,color,newcomm, data_type
+    integer, dimension(3) :: xsz,ysz,zsz,xst,yst,zst,xen,yen,zen,skip
+
+    data_type = real_type
+
+#include "io_write_every.inc"
+
+    return
+  end subroutine write_every_real
+
+
+  subroutine write_every_complex(ipencil,var,iskip,jskip,kskip, &
+       filename, from1)
+
+    implicit none
+
+    integer, intent(IN) :: ipencil !(x-pencil=1; y-pencil=2; z-pencil=3)
+    complex(mytype), dimension(:,:,:), intent(IN) :: var
+    integer, intent(IN) :: iskip,jskip,kskip 
+    character(len=*), intent(IN) :: filename
+    logical, intent(IN) :: from1  ! .true.  - save 1,n+1,2n+1...
+    ! .false. - save n,2n,3n...
+
+    complex(mytype), allocatable, dimension(:,:,:) :: wk, wk2
+    integer(kind=MPI_OFFSET_KIND) :: filesize, disp
+    integer, dimension(3) :: sizes, subsizes, starts
+    integer :: i,j,k, ierror, newtype, fh, key,color,newcomm, data_type
+    integer, dimension(3) :: xsz,ysz,zsz,xst,yst,zst,xen,yen,zen,skip
+
+    data_type = complex_type
+
+#include "io_write_every.inc"
+
+    return
+  end subroutine write_every_complex
+
+  subroutine coarse_extents(ipencil, icoarse, sizes, subsizes, starts)
+
+    implicit none
+    
+    integer, intent(IN) :: ipencil !(x-pencil=1; y-pencil=2; z-pencil=3)
+    integer, intent(IN) :: icoarse !(nstat=1; nvisu=2)
+
+    integer, dimension(3) :: sizes, subsizes, starts
+    integer :: ierror
+
+    if ((icoarse.lt.1).or.(icoarse.gt.2)) then
+       print *, "Error invalid value of icoarse: ", icoarse
+       call MPI_ABORT(MPI_COMM_WORLD, -1, ierror)
+    endif
+    if ((ipencil.lt.1).or.(ipencil.gt.3)) then
+       print *, "Error invalid value of ipencil: ", ipencil
+       call MPI_ABORT(MPI_COMM_WORLD, -1, ierror)
+    endif
+    
+    if (icoarse==1) then
+       sizes(1) = xszS(1)
+       sizes(2) = yszS(2)
+       sizes(3) = zszS(3)
+
+       if (ipencil == 1) then
+          subsizes(1) = xszS(1)
+          subsizes(2) = xszS(2)
+          subsizes(3) = xszS(3)
+          starts(1) = xstS(1)-1  ! 0-based index
+          starts(2) = xstS(2)-1
+          starts(3) = xstS(3)-1
+       else if (ipencil == 2) then
+          subsizes(1) = yszS(1)
+          subsizes(2) = yszS(2)
+          subsizes(3) = yszS(3)
+          starts(1) = ystS(1)-1
+          starts(2) = ystS(2)-1
+          starts(3) = ystS(3)-1
+       else if (ipencil == 3) then
+          subsizes(1) = zszS(1)
+          subsizes(2) = zszS(2)
+          subsizes(3) = zszS(3)
+          starts(1) = zstS(1)-1
+          starts(2) = zstS(2)-1
+          starts(3) = zstS(3)-1
+       endif
+    elseif (icoarse==2) then
+       sizes(1) = xszV(1)
+       sizes(2) = yszV(2)
+       sizes(3) = zszV(3)
+
+       if (ipencil == 1) then
+          subsizes(1) = xszV(1)
+          subsizes(2) = xszV(2)
+          subsizes(3) = xszV(3)
+          starts(1) = xstV(1)-1  ! 0-based index
+          starts(2) = xstV(2)-1
+          starts(3) = xstV(3)-1
+       else if (ipencil == 2) then
+          subsizes(1) = yszV(1)
+          subsizes(2) = yszV(2)
+          subsizes(3) = yszV(3)
+          starts(1) = ystV(1)-1
+          starts(2) = ystV(2)-1
+          starts(3) = ystV(3)-1
+       else if (ipencil == 3) then
+          subsizes(1) = zszV(1)
+          subsizes(2) = zszV(2)
+          subsizes(3) = zszV(3)
+          starts(1) = zstV(1)-1
+          starts(2) = zstV(2)-1
+          starts(3) = zstV(3)-1
+       endif
+    endif
+    
+  end subroutine coarse_extents
+
+  subroutine mpiio_write_real_coarse(ipencil,var,filename,icoarse)
+
+    ! USE param
+    ! USE variables
+
+    implicit none
+
+    integer, intent(IN) :: ipencil !(x-pencil=1; y-pencil=2; z-pencil=3)
+    integer, intent(IN) :: icoarse !(nstat=1; nvisu=2)
+    real(mytype), dimension(:,:,:), intent(IN) :: var
+    real(mytype_single), allocatable, dimension(:,:,:) :: varsingle
+    character(len=*) :: filename
+
+    integer (kind=MPI_OFFSET_KIND) :: filesize, disp
+    integer, dimension(3) :: sizes, subsizes, starts
+    integer :: i,j,k, ierror, newtype, fh
+
+    call coarse_extents(ipencil, icoarse, sizes, subsizes, starts)
+    allocate (varsingle(xstV(1):xenV(1),xstV(2):xenV(2),xstV(3):xenV(3)))
+    varsingle=real(var, mytype_single)
+
+    call MPI_TYPE_CREATE_SUBARRAY(3, sizes, subsizes, starts,  &
+         MPI_ORDER_FORTRAN, real_type_single, newtype, ierror)
+    call MPI_TYPE_COMMIT(newtype,ierror)
+    call MPI_FILE_OPEN(MPI_COMM_WORLD, filename, &
+         MPI_MODE_CREATE+MPI_MODE_WRONLY, MPI_INFO_NULL, &
+         fh, ierror)
+    filesize = 0_MPI_OFFSET_KIND
+    call MPI_FILE_SET_SIZE(fh,filesize,ierror)  ! guarantee overwriting
+    disp = 0_MPI_OFFSET_KIND
+    call MPI_FILE_SET_VIEW(fh,disp,real_type_single, &
+         newtype,'native',MPI_INFO_NULL,ierror)
+    call MPI_FILE_WRITE_ALL(fh, varsingle, &
+         subsizes(1)*subsizes(2)*subsizes(3), &
+         real_type_single, MPI_STATUS_IGNORE, ierror)
+    call MPI_FILE_CLOSE(fh,ierror)
+    call MPI_TYPE_FREE(newtype,ierror)
+
+    deallocate(varsingle)
+
+    return
+  end subroutine mpiio_write_real_coarse
+
+#ifdef ADIOS2
+  subroutine adios2_register_variable(io, varname, ipencil, icoarse)
+
+    implicit none
+
+    integer, intent(IN) :: ipencil !(x-pencil=1; y-pencil=2; z-pencil=3)
+    integer, intent(IN) :: icoarse !(nstat=1; nvisu=2)
+    type(adios2_io), intent(in) :: io
+    
+    character*(*), intent(in) :: varname
+
+    integer, dimension(3) :: sizes, subsizes, starts
+    type(adios2_variable) :: var_handle
+    integer, parameter :: ndims = 3
+    logical, parameter :: adios2_constant_dims = .true.
+    integer :: data_type
+    integer :: ierror
+
+    call coarse_extents(ipencil, icoarse, sizes, subsizes, starts)
+
+    ! Check if variable already exists, if not create it
+    call adios2_inquire_variable(var_handle, io, varname, ierror)
+    if (.not.var_handle % valid) then
+       !! New variable
+       if(nrank==0) print *, "Registering variable for IO: ", varname
+   
+       ! Need to set the ADIOS2 data type
+       if (mytype_single.eq.kind(0.0d0)) then
+          !! Double
+          data_type = adios2_type_dp
+       else if (mytype_single.eq.kind(0.0)) then
+          !! Single
+          data_type = adios2_type_real
+       else
+          print *, "Trying to write unknown data type!"
+          call MPI_ABORT(MPI_COMM_WORLD, -1, ierror)
+       endif
+
+       call adios2_define_variable(var_handle, io, varname, data_type, &
+            ndims, int(sizes, kind=8), int(starts, kind=8), int(subsizes, kind=8), &
+            adios2_constant_dims, ierror)
+    endif
+    
+  end subroutine adios2_register_variable
+  subroutine adios2_write_real_coarse(ipencil,var,varname,icoarse,adios,engine,io)
+
+    ! USE param
+    ! USE variables
+
+    implicit none
+
+    integer, intent(IN) :: ipencil !(x-pencil=1; y-pencil=2; z-pencil=3)
+    integer, intent(IN) :: icoarse !(nstat=1; nvisu=2)
+    real(mytype), dimension(:,:,:), intent(IN) :: var
+    type(adios2_adios), intent(in) :: adios
+    type(adios2_engine), intent(in) :: engine
+    type(adios2_io), intent(in) :: io
+    character*(*), intent(in) :: varname
+    
+    integer (kind=MPI_OFFSET_KIND) :: filesize, disp
+    integer :: i,j,k, ierror, newtype, fh
+    type(adios2_variable) :: var_handle
+
+    call adios2_inquire_variable(var_handle, io, varname, ierror)
+    if (.not.var_handle % valid) then
+       call adios2_register_variable(io, varname, ipencil, icoarse)
+    endif
+
+    call adios2_put(engine, var_handle, var, adios2_mode_deferred, ierror)
+
+    return
+  end subroutine adios2_write_real_coarse
+#endif
+  
+  subroutine mpiio_write_real_probe(ipencil,var,filename,nlength)
+
+    ! USE param
+    ! USE variables
+
+    implicit none
+
+    integer, intent(IN) :: ipencil !(x-pencil=1; y-pencil=2; z-pencil=3)
+    integer, intent(in) :: nlength
+    real(mytype), dimension(:,:,:,:), intent(IN) :: var
+
+    character(len=*) :: filename
+
+    integer (kind=MPI_OFFSET_KIND) :: filesize, disp
+    integer, dimension(4) :: sizes, subsizes, starts
+    integer :: i,j,k, ierror, newtype, fh
+
+    sizes(1) = xszP(1)
+    sizes(2) = yszP(2)
+    sizes(3) = zszP(3)
+    sizes(4) = nlength
+    if (ipencil == 1) then
+       subsizes(1) = xszP(1)
+       subsizes(2) = xszP(2)
+       subsizes(3) = xszP(3)
+       subsizes(4) = nlength
+       starts(1) = xstP(1)-1  ! 0-based index
+       starts(2) = xstP(2)-1
+       starts(3) = xstP(3)-1
+       starts(4) = 0
+    else if (ipencil == 2) then
+       subsizes(1) = yszP(1)
+       subsizes(2) = yszP(2)
+       subsizes(3) = yszP(3)
+       starts(1) = ystP(1)-1
+       starts(2) = ystP(2)-1
+       starts(3) = ystP(3)-1
+    else if (ipencil == 3) then
+       subsizes(1) = zszP(1)
+       subsizes(2) = zszP(2)
+       subsizes(3) = zszP(3)
+       starts(1) = zstP(1)-1
+       starts(2) = zstP(2)-1
+       starts(3) = zstP(3)-1
+    endif
+    !   print *,nrank,starts(1),starts(2),starts(3),starts(4)
+    call MPI_TYPE_CREATE_SUBARRAY(4, sizes, subsizes, starts,  &
+         MPI_ORDER_FORTRAN, real_type, newtype, ierror)
+    call MPI_TYPE_COMMIT(newtype,ierror)
+    call MPI_FILE_OPEN(MPI_COMM_WORLD, filename, &
+         MPI_MODE_CREATE+MPI_MODE_WRONLY, MPI_INFO_NULL, &
+         fh, ierror)
+    filesize = 0_MPI_OFFSET_KIND
+    call MPI_FILE_SET_SIZE(fh,filesize,ierror)  ! guarantee overwriting
+    disp = 0_MPI_OFFSET_KIND
+    call MPI_FILE_SET_VIEW(fh,disp,real_type, &
+         newtype,'native',MPI_INFO_NULL,ierror)
+    call MPI_FILE_WRITE_ALL(fh, var, &
+         subsizes(1)*subsizes(2)*subsizes(3)*subsizes(4), &
+         real_type, MPI_STATUS_IGNORE, ierror)
+    call MPI_FILE_CLOSE(fh,ierror)
+    call MPI_TYPE_FREE(newtype,ierror)
+
+
+    return
+  end subroutine mpiio_write_real_probe
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  ! Write a 3D data set covering a smaller sub-domain only
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  subroutine write_subdomain(ipencil,var,is,ie,js,je,ks,ke,filename)
+
+    implicit none
+
+    integer, intent(IN) :: ipencil !(x-pencil=1; y-pencil=2; z-pencil=3)
+    real(mytype), dimension(:,:,:), intent(IN) :: var
+    integer, intent(IN) :: is, ie, js, je, ks, ke
+    character(len=*), intent(IN) :: filename
+
+    real(mytype), allocatable, dimension(:,:,:) :: wk, wk2
+    integer(kind=MPI_OFFSET_KIND) :: filesize, disp
+    integer, dimension(3) :: sizes, subsizes, starts
+    integer :: color, key, errorcode, newcomm, ierror
+    integer :: newtype, fh, data_type, i, j, k
+    integer :: i1, i2, j1, j2, k1, k2
+
+    data_type = real_type
+
+    ! validate the input paramters
+    if (is<1 .OR. ie>nx_global .OR. js<1 .OR. je>ny_global .OR. &
+         ks<1 .OR. ke>nz_global) then
+       errorcode = 10
+       call decomp_2d_abort(errorcode, &
+            'Invalid subdomain specified in I/O')
+    end if
+
+    ! create a communicator for all those MPI ranks containing the subdomain
+    color = 1
+    key = 1
+    if (ipencil==1) then
+       if (xstart(1)>ie .OR. xend(1)<is .OR. xstart(2)>je .OR. xend(2)<js &
+            .OR. xstart(3)>ke .OR. xend(3)<ks) then
+          color = 2
+       end if
+    else if (ipencil==2) then
+       if (ystart(1)>ie .OR. yend(1)<is .OR. ystart(2)>je .OR. yend(2)<js &
+            .OR. ystart(3)>ke .OR. yend(3)<ks) then
+          color = 2
+       end if
+    else if (ipencil==3) then
+       if (zstart(1)>ie .OR. zend(1)<is .OR. zstart(2)>je .OR. zend(2)<js &
+            .OR. zstart(3)>ke .OR. zend(3)<ks) then
+          color = 2
+       end if
+    end if
+    call MPI_COMM_SPLIT(MPI_COMM_WORLD,color,key,newcomm,ierror)
+
+    if (color==1) then ! only ranks in this group do IO collectively
+
+       ! generate MPI-IO subarray information
+
+       ! global size of the sub-domain to write
+       sizes(1) = ie - is + 1
+       sizes(2) = je - js + 1
+       sizes(3) = ke - ks + 1
+
+       ! 'subsizes' & 'starts' as required by MPI_TYPE_CREATE_SUBARRAY
+       ! note the special code whe subdomain only occupy part of the pencil
+       if (ipencil==1) then
+
+          subsizes(1) = xsize(1)
+          starts(1) = xstart(1) - is
+          if (xend(1)>ie .AND. xstart(1)<is) then
+             subsizes(1) = ie - is + 1
+             starts(1) = 0
+          else if (xstart(1)<is) then
+             subsizes(1) = xend(1) - is + 1
+             starts(1) = 0
+          else if (xend(1)>ie) then
+             subsizes(1) = ie - xstart(1) + 1
+          end if
+          subsizes(2) = xsize(2)
+          starts(2) = xstart(2) - js
+          if (xend(2)>je .AND. xstart(2)<js) then
+             subsizes(2) = je - js + 1
+             starts(2) = 0
+          else if (xstart(2)<js) then
+             subsizes(2) = xend(2) - js + 1
+             starts(2) = 0
+          else if (xend(2)>je) then
+             subsizes(2) = je - xstart(2) + 1
+          end if
+          subsizes(3) = xsize(3)
+          starts(3) = xstart(3) - ks
+          if (xend(3)>ke .AND. xstart(3)<ks) then
+             subsizes(3) = ke - ks + 1
+             starts(3) = 0
+          else if (xstart(3)<ks) then
+             subsizes(3) = xend(3) - ks + 1
+             starts(3) = 0
+          else if (xend(3)>ke) then
+             subsizes(3) = ke - xstart(3) + 1
+          end if
+
+       else if (ipencil==2) then
+
+          ! TODO 
+
+       else if (ipencil==3) then
+
+          ! TODO
+
+       end if
+
+
+       ! copy data from orginal to a temp array
+       ! pay attention to blocks only partially cover the sub-domain
+       if (ipencil==1) then
+
+          if (xend(1)>ie .AND. xstart(1)<is) then
+             i1 = is
+             i2 = ie
+          else if (xend(1)>ie) then
+             i1 = xstart(1)
+             i2 = ie
+          else if (xstart(1)<is) then
+             i1 = is
+             i2 = xend(1)
+          else
+             i1 = xstart(1)
+             i2 = xend(1)
+          end if
+
+          if (xend(2)>je .AND. xstart(2)<js) then
+             j1 = js
+             j2 = je
+          else if (xend(2)>je) then
+             j1 = xstart(2)
+             j2 = je
+          else if (xstart(2)<js) then
+             j1 = js
+             j2 = xend(2)
+          else
+             j1 = xstart(2)
+             j2 = xend(2)
+          end if
+
+          if (xend(3)>ke .AND. xstart(3)<ks) then
+             k1 = ks
+             k2 = ke
+          else if (xend(3)>ke) then
+             k1 = xstart(3)
+             k2 = ke
+          else if (xstart(3)<ks) then
+             k1 = ks
+             k2 = xend(3)
+          else
+             k1 = xstart(3)
+             k2 = xend(3)
+          end if
+
+          allocate(wk(i1:i2, j1:j2, k1:k2))
+          allocate(wk2(xstart(1):xend(1),xstart(2):xend(2),xstart(3):xend(3)))
+          wk2 = var
+          do k=k1,k2
+             do j=j1,j2
+                do i=i1,i2
+                   wk(i,j,k) = wk2(i,j,k)
+                end do
+             end do
+          end do
+
+       else if (ipencil==2) then
+
+          ! TODO
+
+       else if (ipencil==3) then
+
+          ! TODO
+
+       end if
+
+       deallocate(wk2)
+
+       ! MPI-IO
+       call MPI_TYPE_CREATE_SUBARRAY(3, sizes, subsizes, starts,  &
+            MPI_ORDER_FORTRAN, data_type, newtype, ierror)
+       call MPI_TYPE_COMMIT(newtype,ierror)
+       call MPI_FILE_OPEN(newcomm, filename, &
+            MPI_MODE_CREATE+MPI_MODE_WRONLY, MPI_INFO_NULL, &
+            fh, ierror)
+       filesize = 0_MPI_OFFSET_KIND
+       call MPI_FILE_SET_SIZE(fh,filesize,ierror)  ! guarantee overwriting
+       disp = 0_MPI_OFFSET_KIND
+       call MPI_FILE_SET_VIEW(fh,disp,data_type, &
+            newtype,'native',MPI_INFO_NULL,ierror)
+       call MPI_FILE_WRITE_ALL(fh, wk, &
+            subsizes(1)*subsizes(2)*subsizes(3), &
+            data_type, MPI_STATUS_IGNORE, ierror)
+       call MPI_FILE_CLOSE(fh,ierror)
+       call MPI_TYPE_FREE(newtype,ierror)
+
+       deallocate(wk)
+
+    end if
+
+    return
+  end subroutine write_subdomain
+
+
+end module decomp_2d_io
diff --git a/decomp2d/io_read_inflow.f90 b/decomp2d/io_read_inflow.f90
new file mode 100644
index 0000000..f7f310b
--- /dev/null
+++ b/decomp2d/io_read_inflow.f90
@@ -0,0 +1,52 @@
+!=======================================================================
+! This is part of the 2DECOMP&FFT library
+! 
+! 2DECOMP&FFT is a software framework for general-purpose 2D (pencil) 
+! decomposition. It also implements a highly scalable distributed
+! three-dimensional Fast Fourier Transform (FFT).
+!
+! Copyright (C) 2009-2011 Ning Li, the Numerical Algorithms Group (NAG)
+!
+!=======================================================================
+
+! This file contain common code to be included by subroutines 
+! 'read_var_...' in io.f90
+
+  ! Using MPI-IO to read a distributed 3D variable from a file. File 
+  ! operations (open/close) need to be done in calling application. This 
+  ! allows multiple variables to be read from a single file. Together 
+  ! with the corresponding write operation, this is the perfect solution
+  ! for applications to perform restart/checkpointing.
+
+    if (present(opt_decomp)) then
+       decomp = opt_decomp
+    else
+       call get_decomp_info(decomp)
+    end if
+
+    ! Create file type and set file view
+    sizes(1) = ntimesteps 
+    sizes(2) = decomp%ysz(2)
+    sizes(3) = decomp%zsz(3)
+    subsizes(1) = ntimesteps
+    subsizes(2) = decomp%xsz(2)
+    subsizes(3) = decomp%xsz(3)
+    starts(1) = 0  ! 0-based index
+    starts(2) = decomp%xst(2)-1
+    starts(3) = decomp%xst(3)-1
+    
+    call MPI_TYPE_CREATE_SUBARRAY(3, sizes, subsizes, starts,  &
+         MPI_ORDER_FORTRAN, data_type, newtype, ierror)
+    call MPI_TYPE_COMMIT(newtype,ierror)
+    call MPI_FILE_SET_VIEW(fh,disp,data_type, &
+         newtype,'native',MPI_INFO_NULL,ierror)
+    call MPI_FILE_READ_ALL(fh, var, &
+         subsizes(1)*subsizes(2)*subsizes(3), &
+         data_type, MPI_STATUS_IGNORE, ierror)
+    call MPI_TYPE_FREE(newtype,ierror)
+
+    ! update displacement for the next read operation
+    disp = disp + sizes(1)*sizes(2)*sizes(3)*mytype_bytes
+    if (data_type == complex_type) then
+       disp = disp + sizes(1)*sizes(2)*sizes(3)*mytype_bytes
+    end if
diff --git a/decomp2d/io_read_one.inc b/decomp2d/io_read_one.inc
new file mode 100644
index 0000000..aaa7711
--- /dev/null
+++ b/decomp2d/io_read_one.inc
@@ -0,0 +1,65 @@
+!                          -*- mode: f90 -*-                            
+!=======================================================================
+! This is part of the 2DECOMP&FFT library
+! 
+! 2DECOMP&FFT is a software framework for general-purpose 2D (pencil) 
+! decomposition. It also implements a highly scalable distributed
+! three-dimensional Fast Fourier Transform (FFT).
+!
+! Copyright (C) 2009-2011 Ning Li, the Numerical Algorithms Group (NAG)
+!
+!=======================================================================
+
+! This file contain common code to be included by subroutines 
+! 'mpiio_read_one_...' in io.f90
+
+! Using MPI-IO to write a distributed 3D array into a file
+
+if (present(opt_decomp)) then
+   decomp = opt_decomp
+else
+   call get_decomp_info(decomp)
+end if
+
+! determine subarray parameters
+sizes(1) = decomp%xsz(1)
+sizes(2) = decomp%ysz(2)
+sizes(3) = decomp%zsz(3)
+
+if (ipencil == 1) then
+   subsizes(1) = decomp%xsz(1)
+   subsizes(2) = decomp%xsz(2)
+   subsizes(3) = decomp%xsz(3)
+   starts(1) = decomp%xst(1)-1  ! 0-based index
+   starts(2) = decomp%xst(2)-1
+   starts(3) = decomp%xst(3)-1
+else if (ipencil == 2) then
+   subsizes(1) = decomp%ysz(1)
+   subsizes(2) = decomp%ysz(2)
+   subsizes(3) = decomp%ysz(3)
+   starts(1) = decomp%yst(1)-1
+   starts(2) = decomp%yst(2)-1
+   starts(3) = decomp%yst(3)-1
+else if (ipencil == 3) then
+   subsizes(1) = decomp%zsz(1)
+   subsizes(2) = decomp%zsz(2)
+   subsizes(3) = decomp%zsz(3)
+   starts(1) = decomp%zst(1)-1
+   starts(2) = decomp%zst(2)-1
+   starts(3) = decomp%zst(3)-1
+endif
+
+call MPI_TYPE_CREATE_SUBARRAY(3, sizes, subsizes, starts,  &
+     MPI_ORDER_FORTRAN, data_type, newtype, ierror)
+call MPI_TYPE_COMMIT(newtype,ierror)
+call MPI_FILE_OPEN(MPI_COMM_WORLD, filename, &
+     MPI_MODE_RDONLY, MPI_INFO_NULL, &
+     fh, ierror)
+disp = 0_MPI_OFFSET_KIND
+call MPI_FILE_SET_VIEW(fh,disp,data_type, &
+     newtype,'native',MPI_INFO_NULL,ierror)
+call MPI_FILE_READ_ALL(fh, var, &
+     subsizes(1)*subsizes(2)*subsizes(3), &
+     data_type, MPI_STATUS_IGNORE, ierror)
+call MPI_FILE_CLOSE(fh,ierror)
+call MPI_TYPE_FREE(newtype,ierror)
diff --git a/decomp2d/io_read_var.inc b/decomp2d/io_read_var.inc
new file mode 100644
index 0000000..398f46a
--- /dev/null
+++ b/decomp2d/io_read_var.inc
@@ -0,0 +1,69 @@
+!                          -*- mode: f90 -*-                            
+!=======================================================================
+! This is part of the 2DECOMP&FFT library
+! 
+! 2DECOMP&FFT is a software framework for general-purpose 2D (pencil) 
+! decomposition. It also implements a highly scalable distributed
+! three-dimensional Fast Fourier Transform (FFT).
+!
+! Copyright (C) 2009-2011 Ning Li, the Numerical Algorithms Group (NAG)
+!
+!=======================================================================
+
+! This file contain common code to be included by subroutines 
+! 'read_var_...' in io.f90
+
+! Using MPI-IO to read a distributed 3D variable from a file. File 
+! operations (open/close) need to be done in calling application. This 
+! allows multiple variables to be read from a single file. Together 
+! with the corresponding write operation, this is the perfect solution
+! for applications to perform restart/checkpointing.
+
+if (present(opt_decomp)) then
+   decomp = opt_decomp
+else
+   call get_decomp_info(decomp)
+end if
+
+! Create file type and set file view
+sizes(1) = decomp%xsz(1)
+sizes(2) = decomp%ysz(2)
+sizes(3) = decomp%zsz(3)
+if (ipencil == 1) then
+   subsizes(1) = decomp%xsz(1)
+   subsizes(2) = decomp%xsz(2)
+   subsizes(3) = decomp%xsz(3)
+   starts(1) = decomp%xst(1)-1  ! 0-based index
+   starts(2) = decomp%xst(2)-1
+   starts(3) = decomp%xst(3)-1
+else if (ipencil == 2) then
+   subsizes(1) = decomp%ysz(1)
+   subsizes(2) = decomp%ysz(2)
+   subsizes(3) = decomp%ysz(3)
+   starts(1) = decomp%yst(1)-1
+   starts(2) = decomp%yst(2)-1
+   starts(3) = decomp%yst(3)-1
+else if (ipencil == 3) then
+   subsizes(1) = decomp%zsz(1)
+   subsizes(2) = decomp%zsz(2)
+   subsizes(3) = decomp%zsz(3)
+   starts(1) = decomp%zst(1)-1
+   starts(2) = decomp%zst(2)-1
+   starts(3) = decomp%zst(3)-1
+endif
+
+call MPI_TYPE_CREATE_SUBARRAY(3, sizes, subsizes, starts,  &
+     MPI_ORDER_FORTRAN, data_type, newtype, ierror)
+call MPI_TYPE_COMMIT(newtype,ierror)
+call MPI_FILE_SET_VIEW(fh,disp,data_type, &
+     newtype,'native',MPI_INFO_NULL,ierror)
+call MPI_FILE_READ_ALL(fh, var, &
+     subsizes(1)*subsizes(2)*subsizes(3), &
+     data_type, MPI_STATUS_IGNORE, ierror)
+call MPI_TYPE_FREE(newtype,ierror)
+
+! update displacement for the next read operation
+disp = disp + sizes(1)*sizes(2)*sizes(3)*mytype_bytes
+if (data_type == complex_type) then
+   disp = disp + sizes(1)*sizes(2)*sizes(3)*mytype_bytes
+end if
diff --git a/decomp2d/io_write_every.inc b/decomp2d/io_write_every.inc
new file mode 100644
index 0000000..ec55c1c
--- /dev/null
+++ b/decomp2d/io_write_every.inc
@@ -0,0 +1,224 @@
+!                          -*- mode: f90 -*-                            
+!=======================================================================
+! This is part of the 2DECOMP&FFT library
+! 
+! 2DECOMP&FFT is a software framework for general-purpose 2D (pencil) 
+! decomposition. It also implements a highly scalable distributed
+! three-dimensional Fast Fourier Transform (FFT).
+!
+! Copyright (C) 2009-2011 Ning Li, the Numerical Algorithms Group (NAG)
+!
+!=======================================================================
+
+! This file contain common code to be included by subroutines 
+! 'write_every_...' in io.f90
+
+! To write every few points of a 3D array to a file
+
+! work out the distribution parameters, which may be different from 
+! the default distribution used by the decomposition library
+!  For exmample if nx=17 and p_row=4
+!    distribution is: 4 4 4 5
+
+! If writing from the 1st element
+!  If saving every 3 points, then 5 points to be saved (17/3)
+!    default distribution would be 1 1 1 2
+!    However, 1st block (1-4) contains the 3rd point
+!             2nd block (5-8) contains the 6th point
+!             3rd block (9-12) contains the 9th and 12th point
+!             4th block (13-17) contains then 15th point
+!    giving a 1 1 2 1 distribution
+!    So cannot use the base decomposition library for such IO
+
+! If writing from the n-th element (n=?skip)
+!  If saving every 3 points, then 6 points to be saved
+!    However, 1st block (1-4) contains the 1st & 4th point
+!             2nd block (5-8) contains the 7th point
+!             3rd block (9-12) contains the 10th point
+!             4th block (13-17) contains then 12th & 15th point
+!    giving a 1 2 2 1 distribution
+
+skip(1)=iskip
+skip(2)=jskip
+skip(3)=kskip
+
+do i=1,3
+   if (from1) then
+      xst(i) = (xstart(i)+skip(i)-1)/skip(i)
+      if (mod(xstart(i)+skip(i)-1,skip(i))/=0) xst(i)=xst(i)+1
+      xen(i) = (xend(i)+skip(i)-1)/skip(i)
+   else
+      xst(i) = xstart(i)/skip(i)
+      if (mod(xstart(i),skip(i))/=0) xst(i)=xst(i)+1
+      xen(i) = xend(i)/skip(i)
+   end if
+   xsz(i) = xen(i)-xst(i)+1
+end do
+
+do i=1,3
+   if (from1) then
+      yst(i) = (ystart(i)+skip(i)-1)/skip(i)
+      if (mod(ystart(i)+skip(i)-1,skip(i))/=0) yst(i)=yst(i)+1
+      yen(i) = (yend(i)+skip(i)-1)/skip(i)
+   else
+      yst(i) = ystart(i)/skip(i)
+      if (mod(ystart(i),skip(i))/=0) yst(i)=yst(i)+1
+      yen(i) = yend(i)/skip(i)
+   end if
+   ysz(i) = yen(i)-yst(i)+1
+end do
+
+do i=1,3
+   if (from1) then
+      zst(i) = (zstart(i)+skip(i)-1)/skip(i)
+      if (mod(zstart(i)+skip(i)-1,skip(i))/=0) zst(i)=zst(i)+1
+      zen(i) = (zend(i)+skip(i)-1)/skip(i)
+   else
+      zst(i) = zstart(i)/skip(i)
+      if (mod(zstart(i),skip(i))/=0) zst(i)=zst(i)+1
+      zen(i) = zend(i)/skip(i)
+   end if
+   zsz(i) = zen(i)-zst(i)+1
+end do
+
+! if 'skip' value is large it is possible that some ranks do not 
+! contain any points to be written. Subarray constructor requires 
+! nonzero size so it is not possible to use MPI_COMM_WORLD for IO.
+! Create a sub communicator for this...
+color = 1
+key = 0  ! rank order doesn't matter
+if (ipencil==1) then
+   if (xsz(1)==0 .or. xsz(2)==0 .or. xsz(3)==0) then
+      color = 2
+   end if
+else if (ipencil==2) then
+   if (ysz(1)==0 .or. ysz(2)==0 .or. ysz(3)==0) then
+      color = 2
+   end if
+else if (ipencil==3) then
+   if (zsz(1)==0 .or. zsz(2)==0 .or. zsz(3)==0) then
+      color = 2
+   end if
+end if
+call MPI_COMM_SPLIT(MPI_COMM_WORLD,color,key,newcomm,ierror)
+
+if (color==1) then ! only ranks in this group do IO collectively
+
+   ! generate subarray information
+   sizes(1) = xsz(1)
+   sizes(2) = ysz(2)
+   sizes(3) = zsz(3)
+   if (ipencil==1) then
+      subsizes(1) = xsz(1)
+      subsizes(2) = xsz(2)
+      subsizes(3) = xsz(3)
+      starts(1) = xst(1)-1
+      starts(2) = xst(2)-1
+      starts(3) = xst(3)-1
+   else if (ipencil==2) then
+      subsizes(1) = ysz(1)
+      subsizes(2) = ysz(2)
+      subsizes(3) = ysz(3)
+      starts(1) = yst(1)-1
+      starts(2) = yst(2)-1
+      starts(3) = yst(3)-1
+   else if (ipencil==3) then
+      subsizes(1) = zsz(1)
+      subsizes(2) = zsz(2)
+      subsizes(3) = zsz(3)
+      starts(1) = zst(1)-1
+      starts(2) = zst(2)-1
+      starts(3) = zst(3)-1
+   end if
+
+   ! copy data from original array
+   ! needs a copy of original array in global coordinate 
+   if (ipencil==1) then
+      allocate(wk(xst(1):xen(1),xst(2):xen(2),xst(3):xen(3)))
+      allocate(wk2(xstart(1):xend(1),xstart(2):xend(2),xstart(3):xend(3)))
+      wk2=var
+      if (from1) then
+         do k=xst(3),xen(3)
+            do j=xst(2),xen(2)
+               do i=xst(1),xen(1)
+                  wk(i,j,k) = wk2((i-1)*iskip+1,(j-1)*jskip+1,(k-1)*kskip+1)
+               end do
+            end do
+         end do
+      else
+         do k=xst(3),xen(3)
+            do j=xst(2),xen(2)
+               do i=xst(1),xen(1)
+                  wk(i,j,k) = wk2(i*iskip,j*jskip,k*kskip)
+               end do
+            end do
+         end do
+      end if
+   else if (ipencil==2) then
+      allocate(wk(yst(1):yen(1),yst(2):yen(2),yst(3):yen(3)))
+      allocate(wk2(ystart(1):yend(1),ystart(2):yend(2),ystart(3):yend(3)))
+      wk2=var
+      if (from1) then
+         do k=yst(3),yen(3)
+            do j=yst(2),yen(2)
+               do i=yst(1),yen(1)
+                  wk(i,j,k) = wk2((i-1)*iskip+1,(j-1)*jskip+1,(k-1)*kskip+1)
+               end do
+            end do
+         end do
+      else
+         do k=yst(3),yen(3)
+            do j=yst(2),yen(2)
+               do i=yst(1),yen(1)
+                  wk(i,j,k) = wk2(i*iskip,j*jskip,k*kskip)
+               end do
+            end do
+         end do
+      end if
+   else if (ipencil==3) then
+      allocate(wk(zst(1):zen(1),zst(2):zen(2),zst(3):zen(3)))
+      allocate(wk2(zstart(1):zend(1),zstart(2):zend(2),zstart(3):zend(3)))
+      wk2=var
+      if (from1) then
+         do k=zst(3),zen(3)
+            do j=zst(2),zen(2)
+               do i=zst(1),zen(1)
+                  wk(i,j,k) = wk2((i-1)*iskip+1,(j-1)*jskip+1,(k-1)*kskip+1)
+               end do
+            end do
+         end do
+      else
+         do k=zst(3),zen(3)
+            do j=zst(2),zen(2)
+               do i=zst(1),zen(1)
+                  wk(i,j,k) = wk2(i*iskip,j*jskip,k*kskip)
+               end do
+            end do
+         end do
+      end if
+   end if
+   deallocate(wk2)
+
+   ! MPI-IO
+   call MPI_TYPE_CREATE_SUBARRAY(3, sizes, subsizes, starts,  &
+        MPI_ORDER_FORTRAN, data_type, newtype, ierror)
+   call MPI_TYPE_COMMIT(newtype,ierror)
+   call MPI_FILE_OPEN(newcomm, filename, &
+        MPI_MODE_CREATE+MPI_MODE_WRONLY, MPI_INFO_NULL, &
+        fh, ierror)
+   filesize = 0_MPI_OFFSET_KIND
+   call MPI_FILE_SET_SIZE(fh,filesize,ierror)  ! guarantee overwriting
+   disp = 0_MPI_OFFSET_KIND
+   call MPI_FILE_SET_VIEW(fh,disp,data_type, &
+        newtype,'native',MPI_INFO_NULL,ierror)
+   call MPI_FILE_WRITE_ALL(fh, wk, &
+        subsizes(1)*subsizes(2)*subsizes(3), &
+        data_type, MPI_STATUS_IGNORE, ierror)
+   call MPI_FILE_CLOSE(fh,ierror)
+   call MPI_TYPE_FREE(newtype,ierror)
+
+   deallocate(wk)
+
+end if ! color==1
+
+call MPI_BARRIER(MPI_COMM_WORLD, ierror)
diff --git a/decomp2d/io_write_one.inc b/decomp2d/io_write_one.inc
new file mode 100644
index 0000000..f7d8692
--- /dev/null
+++ b/decomp2d/io_write_one.inc
@@ -0,0 +1,84 @@
+!                          -*- mode: f90 -*-                            
+!=======================================================================
+! This is part of the 2DECOMP&FFT library
+! 
+! 2DECOMP&FFT is a software framework for general-purpose 2D (pencil) 
+! decomposition. It also implements a highly scalable distributed
+! three-dimensional Fast Fourier Transform (FFT).
+!
+! Copyright (C) 2009-2012 Ning Li, the Numerical Algorithms Group (NAG) 
+! Copyright (C) 2021               the University of Edinburgh (UoE)
+!
+!=======================================================================
+
+! This file contain common code to be included by subroutines 
+! 'mpiio_write_one_...' in io.f90
+
+! Using MPI-IO to write a distributed 3D array into a file
+
+if (present(opt_decomp)) then
+   decomp = opt_decomp
+else
+   call get_decomp_info(decomp)
+end if
+
+! determine subarray parameters
+sizes(1) = decomp%xsz(1)
+sizes(2) = decomp%ysz(2)
+sizes(3) = decomp%zsz(3)
+
+if (ipencil == 1) then
+   subsizes(1) = decomp%xsz(1)
+   subsizes(2) = decomp%xsz(2)
+   subsizes(3) = decomp%xsz(3)
+   starts(1) = decomp%xst(1)-1  ! 0-based index
+   starts(2) = decomp%xst(2)-1
+   starts(3) = decomp%xst(3)-1
+else if (ipencil == 2) then
+   subsizes(1) = decomp%ysz(1)
+   subsizes(2) = decomp%ysz(2)
+   subsizes(3) = decomp%ysz(3)
+   starts(1) = decomp%yst(1)-1
+   starts(2) = decomp%yst(2)-1
+   starts(3) = decomp%yst(3)-1
+else if (ipencil == 3) then
+   subsizes(1) = decomp%zsz(1)
+   subsizes(2) = decomp%zsz(2)
+   subsizes(3) = decomp%zsz(3)
+   starts(1) = decomp%zst(1)-1
+   starts(2) = decomp%zst(2)-1
+   starts(3) = decomp%zst(3)-1
+endif
+
+#ifdef T3PIO
+call MPI_INFO_CREATE(info, ierror)
+gs = ceiling(real(sizes(1),mytype)*real(sizes(2),mytype)* &
+     real(sizes(3),mytype)/1024./1024.)
+call t3pio_set_info(MPI_COMM_WORLD, info, "./", ierror, &
+     GLOBAL_SIZE=gs, factor=1)
+#endif
+
+call MPI_TYPE_CREATE_SUBARRAY(3, sizes, subsizes, starts,  &
+     MPI_ORDER_FORTRAN, data_type, newtype, ierror)
+call MPI_TYPE_COMMIT(newtype,ierror)
+#ifdef T3PIO
+call MPI_FILE_OPEN(MPI_COMM_WORLD, filename, &
+     MPI_MODE_CREATE+MPI_MODE_WRONLY, info, fh, ierror)
+#else
+call MPI_FILE_OPEN(MPI_COMM_WORLD, filename, &
+     MPI_MODE_CREATE+MPI_MODE_WRONLY, MPI_INFO_NULL, &
+     fh, ierror)
+#endif
+filesize = 0_MPI_OFFSET_KIND
+call MPI_FILE_SET_SIZE(fh,filesize,ierror)  ! guarantee overwriting
+disp = 0_MPI_OFFSET_KIND
+call MPI_FILE_SET_VIEW(fh,disp,data_type, &
+     newtype,'native',MPI_INFO_NULL,ierror)
+call MPI_FILE_WRITE_ALL(fh, var, &
+     subsizes(1)*subsizes(2)*subsizes(3), &
+     data_type, MPI_STATUS_IGNORE, ierror)
+call MPI_FILE_CLOSE(fh,ierror)
+call MPI_TYPE_FREE(newtype,ierror)
+#ifdef T3PIO
+call MPI_INFO_FREE(info,ierror)
+#endif
diff --git a/decomp2d/io_write_outflow.f90 b/decomp2d/io_write_outflow.f90
new file mode 100644
index 0000000..1713848
--- /dev/null
+++ b/decomp2d/io_write_outflow.f90
@@ -0,0 +1,52 @@
+!=======================================================================
+! This is part of the 2DECOMP&FFT library
+! 
+! 2DECOMP&FFT is a software framework for general-purpose 2D (pencil) 
+! decomposition. It also implements a highly scalable distributed
+! three-dimensional Fast Fourier Transform (FFT).
+!
+! Copyright (C) 2009-2011 Ning Li, the Numerical Algorithms Group (NAG)
+!
+!=======================================================================
+
+! This file contain common code to be included by subroutines 
+! 'write_var_...' in io.f90
+
+  ! Using MPI-IO to write a distributed 3D variable to a file. File 
+  ! operations (open/close) need to be done in calling application. This
+  ! allows multiple variables to be written to a single file. Together 
+  ! with the corresponding read operation, this is the perfect solution
+  ! for applications to perform restart/checkpointing.
+
+    if (present(opt_decomp)) then
+       decomp = opt_decomp
+    else
+       call get_decomp_info(decomp)
+    end if
+
+    ! Create file type and set file view
+    sizes(1) = ntimesteps 
+    sizes(2) = decomp%ysz(2)
+    sizes(3) = decomp%zsz(3)
+    subsizes(1) = ntimesteps 
+    subsizes(2) = decomp%xsz(2)
+    subsizes(3) = decomp%xsz(3)
+    starts(1) = 0  ! 0-based index
+    starts(2) = decomp%xst(2)-1
+    starts(3) = decomp%xst(3)-1
+   
+    call MPI_TYPE_CREATE_SUBARRAY(3, sizes, subsizes, starts,  &
+         MPI_ORDER_FORTRAN, data_type, newtype, ierror)
+    call MPI_TYPE_COMMIT(newtype,ierror)
+    call MPI_FILE_SET_VIEW(fh,disp,data_type, &
+         newtype,'native',MPI_INFO_NULL,ierror)
+    call MPI_FILE_WRITE_ALL(fh, var, &
+         subsizes(1)*subsizes(2)*subsizes(3), &
+         data_type, MPI_STATUS_IGNORE, ierror)
+    call MPI_TYPE_FREE(newtype,ierror)
+
+    ! update displacement for the next write operation
+    disp = disp + sizes(1)*sizes(2)*sizes(3)*mytype_bytes
+    if (data_type == complex_type) then
+       disp = disp + sizes(1)*sizes(2)*sizes(3)*mytype_bytes
+    end if
diff --git a/decomp2d/io_write_plane.inc b/decomp2d/io_write_plane.inc
new file mode 100644
index 0000000..ed43fbe
--- /dev/null
+++ b/decomp2d/io_write_plane.inc
@@ -0,0 +1,148 @@
+!                          -*- mode: f90 -*-                            
+!=======================================================================
+! This is part of the 2DECOMP&FFT library
+! 
+! 2DECOMP&FFT is a software framework for general-purpose 2D (pencil) 
+! decomposition. It also implements a highly scalable distributed
+! three-dimensional Fast Fourier Transform (FFT).
+!
+! Copyright (C) 2009-2011 Ning Li, the Numerical Algorithms Group (NAG)
+!
+!=======================================================================
+
+! This file contain common code to be included by subroutines 
+! 'mpiio_write_plane_3d_...' in io.f90
+
+! It is much easier to implement if all mpi ranks participate I/O.
+! Transpose the 3D data if necessary.
+
+if (present(opt_decomp)) then
+   decomp = opt_decomp
+else
+   call get_decomp_info(decomp)
+end if
+
+if (iplane==1) then
+allocate(wk(decomp%xsz(1),decomp%xsz(2),decomp%xsz(3)))
+if (ipencil==1) then
+wk = var
+else if (ipencil==2) then
+call transpose_y_to_x(var,wk,decomp)
+else if (ipencil==3) then
+allocate(wk2(decomp%ysz(1),decomp%ysz(2),decomp%ysz(3)))
+call transpose_z_to_y(var,wk2,decomp)
+call transpose_y_to_x(wk2,wk,decomp)
+deallocate(wk2)
+end if
+allocate(wk2d(1,decomp%xsz(2),decomp%xsz(3)))
+if (n.ge.1) then
+  do k=1,decomp%xsz(3)
+  do j=1,decomp%xsz(2)
+  wk2d(1,j,k)=wk(n,j,k)
+  end do
+  end do
+else
+  do k=1,decomp%xsz(3)
+  do j=1,decomp%xsz(2)
+  wk2d(1,j,k)=sum(wk(:,j,k))/real(decomp%xsz(1),kind=mytype)
+  end do
+  end do
+endif
+sizes(1) = 1
+sizes(2) = decomp%ysz(2)
+sizes(3) = decomp%zsz(3)
+subsizes(1) = 1
+subsizes(2) = decomp%xsz(2)
+subsizes(3) = decomp%xsz(3)
+starts(1) = 0
+starts(2) = decomp%xst(2)-1
+starts(3) = decomp%xst(3)-1
+
+else if (iplane==2) then
+allocate(wk(decomp%ysz(1),decomp%ysz(2),decomp%ysz(3)))
+if (ipencil==1) then
+call transpose_x_to_y(var,wk,decomp)
+else if (ipencil==2) then
+wk = var
+else if (ipencil==3) then
+call transpose_z_to_y(var,wk,decomp)
+end if
+allocate(wk2d(decomp%ysz(1),1,decomp%ysz(3)))
+if (n.ge.1) then
+  do k=1,decomp%ysz(3)
+  do i=1,decomp%ysz(1)
+  wk2d(i,1,k)=wk(i,n,k)
+  end do
+  end do
+else
+  do k=1,decomp%ysz(3)
+  do i=1,decomp%ysz(1)
+    wk2d(i,1,k)=sum(wk(i,:,k))/real(decomp%ysz(2),kind=mytype)
+  end do
+  end do
+endif
+sizes(1) = decomp%xsz(1)
+sizes(2) = 1
+sizes(3) = decomp%zsz(3)
+subsizes(1) = decomp%ysz(1)
+subsizes(2) = 1
+subsizes(3) = decomp%ysz(3)
+starts(1) = decomp%yst(1)-1
+starts(2) = 0
+starts(3) = decomp%yst(3)-1
+
+else if (iplane==3) then
+allocate(wk(decomp%zsz(1),decomp%zsz(2),decomp%zsz(3)))
+if (ipencil==1) then
+allocate(wk2(decomp%ysz(1),decomp%ysz(2),decomp%ysz(3)))
+call transpose_x_to_y(var,wk2,decomp)
+call transpose_y_to_z(wk2,wk,decomp)
+deallocate(wk2)
+else if (ipencil==2) then
+call transpose_y_to_z(var,wk,decomp)
+else if (ipencil==3) then
+wk = var
+end if
+allocate(wk2d(decomp%zsz(1),decomp%zsz(2),1))
+if (n.ge.1) then
+  do j=1,decomp%zsz(2)
+  do i=1,decomp%zsz(1) 
+  wk2d(i,j,1)=wk(i,j,n)
+  end do
+  end do
+else
+  do j=1,decomp%zsz(2)
+  do i=1,decomp%zsz(1)
+  wk2d(i,j,1)=sum(wk(i,j,:))/real(decomp%zsz(3),kind=mytype)
+  end do
+  end do
+endif
+sizes(1) = decomp%xsz(1)
+sizes(2) = decomp%ysz(2)
+sizes(3) = 1
+subsizes(1) = decomp%zsz(1)
+subsizes(2) = decomp%zsz(2)
+subsizes(3) = 1
+starts(1) = decomp%zst(1)-1
+starts(2) = decomp%zst(2)-1
+starts(3) = 0
+end if
+
+call MPI_TYPE_CREATE_SUBARRAY(3, sizes, subsizes, starts,  &
+     MPI_ORDER_FORTRAN, data_type, newtype, ierror)
+call MPI_TYPE_COMMIT(newtype,ierror)
+call MPI_FILE_OPEN(MPI_COMM_WORLD, filename, &
+     MPI_MODE_CREATE+MPI_MODE_WRONLY, MPI_INFO_NULL, &
+     fh, ierror)
+filesize = 0_MPI_OFFSET_KIND
+call MPI_FILE_SET_SIZE(fh,filesize,ierror)  ! guarantee overwriting
+disp = 0_MPI_OFFSET_KIND
+call MPI_FILE_SET_VIEW(fh,disp,data_type, &
+     newtype,'native',MPI_INFO_NULL,ierror)
+call MPI_FILE_WRITE_ALL(fh, wk2d, &
+     subsizes(1)*subsizes(2)*subsizes(3), &
+     data_type, MPI_STATUS_IGNORE, ierror)
+call MPI_FILE_CLOSE(fh,ierror)
+call MPI_TYPE_FREE(newtype,ierror)
+
+deallocate(wk,wk2d)
diff --git a/decomp2d/io_write_var.inc b/decomp2d/io_write_var.inc
new file mode 100644
index 0000000..c511bf4
--- /dev/null
+++ b/decomp2d/io_write_var.inc
@@ -0,0 +1,69 @@
+!                          -*- mode: f90 -*-                            
+!=======================================================================
+! This is part of the 2DECOMP&FFT library
+! 
+! 2DECOMP&FFT is a software framework for general-purpose 2D (pencil) 
+! decomposition. It also implements a highly scalable distributed
+! three-dimensional Fast Fourier Transform (FFT).
+!
+! Copyright (C) 2009-2011 Ning Li, the Numerical Algorithms Group (NAG)
+!
+!=======================================================================
+
+! This file contain common code to be included by subroutines 
+! 'write_var_...' in io.f90
+
+! Using MPI-IO to write a distributed 3D variable to a file. File 
+! operations (open/close) need to be done in calling application. This
+! allows multiple variables to be written to a single file. Together 
+! with the corresponding read operation, this is the perfect solution
+! for applications to perform restart/checkpointing.
+
+if (present(opt_decomp)) then
+   decomp = opt_decomp
+else
+   call get_decomp_info(decomp)
+end if
+
+! Create file type and set file view
+sizes(1) = decomp%xsz(1)
+sizes(2) = decomp%ysz(2)
+sizes(3) = decomp%zsz(3)
+if (ipencil == 1) then
+   subsizes(1) = decomp%xsz(1)
+   subsizes(2) = decomp%xsz(2)
+   subsizes(3) = decomp%xsz(3)
+   starts(1) = decomp%xst(1)-1  ! 0-based index
+   starts(2) = decomp%xst(2)-1
+   starts(3) = decomp%xst(3)-1
+else if (ipencil == 2) then
+   subsizes(1) = decomp%ysz(1)
+   subsizes(2) = decomp%ysz(2)
+   subsizes(3) = decomp%ysz(3)
+   starts(1) = decomp%yst(1)-1
+   starts(2) = decomp%yst(2)-1
+   starts(3) = decomp%yst(3)-1
+else if (ipencil == 3) then
+   subsizes(1) = decomp%zsz(1)
+   subsizes(2) = decomp%zsz(2)
+   subsizes(3) = decomp%zsz(3)
+   starts(1) = decomp%zst(1)-1
+   starts(2) = decomp%zst(2)-1
+   starts(3) = decomp%zst(3)-1
+endif
+
+call MPI_TYPE_CREATE_SUBARRAY(3, sizes, subsizes, starts,  &
+     MPI_ORDER_FORTRAN, data_type, newtype, ierror)
+call MPI_TYPE_COMMIT(newtype,ierror)
+call MPI_FILE_SET_VIEW(fh,disp,data_type, &
+     newtype,'native',MPI_INFO_NULL,ierror)
+call MPI_FILE_WRITE_ALL(fh, var, &
+     subsizes(1)*subsizes(2)*subsizes(3), &
+     data_type, MPI_STATUS_IGNORE, ierror)
+call MPI_TYPE_FREE(newtype,ierror)
+
+! update displacement for the next write operation
+disp = disp + sizes(1)*sizes(2)*sizes(3)*mytype_bytes
+if (data_type == complex_type) then
+   disp = disp + sizes(1)*sizes(2)*sizes(3)*mytype_bytes
+end if
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 38a8ecf..71b5df7 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -17,6 +17,7 @@ set(X3DFILES
 	    tools.f90
             transeq.f90
             var.f90
+	    visu.f90
 	    x3d_derive.f90
 	    x3d_filters.f90
 	    x3d_operator_1d.f90
diff --git a/src/Makefile.am b/src/Makefile.am
index b6440c9..5228064 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -8,6 +8,7 @@ time_integrators.f90 \
 tools.f90 \
 x3d_transpose.f90 \
 var.f90 \
+visu.f90 \
 thomas.f90 \
 x3d_operator_x_data.f90 \
 x3d_operator_y_data.f90 \
@@ -31,8 +32,9 @@ time_integrators.lo: module_param.lo
 tools.lo: module_param.lo
 x3d_transpose.lo: module_param.lo
 var.lo: module_param.lo
+visu.lo: module_param.lo
 thomas.lo: module_param.lo
-x3d_operator_1d.f90: module_param.lo x3d_operator_x_data.lo x3d_operator_y_data.lo x3d_operator_z_data.lo
+x3d_operator_1d.lo: module_param.lo x3d_operator_x_data.lo x3d_operator_y_data.lo x3d_operator_z_data.lo
 poisson.lo: x3d_operator_x_data.lo x3d_operator_y_data.lo x3d_operator_z_data.lo
 x3d_derive.lo: x3d_operator_1d.lo
 x3d_staggered.lo: x3d_operator_1d.lo
@@ -40,10 +42,10 @@ x3d_filters.lo: x3d_operator_1d.lo
 navier.lo: poisson.lo x3d_staggered.lo
 parameters.lo: module_param.lo
 bc_tgv2d.lo: x3d_operator_x_data.lo x3d_operator_y_data.lo
-bc_tgv.lo: x3d_derive.lo
+bc_tgv.lo: x3d_derive.lo visu.lo
 case.lo: bc_tgv2d.lo bc_tgv.lo
 transeq.lo: x3d_derive.lo case.lo
-x3d_tools.lo: x3d_derive.lo case.lo
+x3d_tools.lo: x3d_derive.lo case.lo visu.lo
 
 clean-local:
 	-rm -f *.mod
diff --git a/src/bc_tgv.f90 b/src/bc_tgv.f90
index ec5b9d4..55f132a 100644
--- a/src/bc_tgv.f90
+++ b/src/bc_tgv.f90
@@ -25,6 +25,8 @@ module bc_tgv
             tgv_listing, &
             tgv_init, &
             tgv_postprocess, &
+            tgv_visu_init, &
+            tgv_visu, &
             tgv_finalize
 
 contains
@@ -129,7 +131,7 @@ subroutine tgv_init(ux1, uy1, uz1)
   !
   subroutine tgv_postprocess(ux1, uy1, uz1, ndt)
 
-    use decomp_2d, only : mytype, real_type, decomp_2d_abort, xsize, ysize, zsize
+    use decomp_2d, only : decomp_2d_abort, ysize, zsize
     use param, only : half, two, xnu, dt
     use variables, only : nx, ny, nz, ppy
     use var, only : ux2,uy2,uz2,ux3,uy3,uz3
@@ -254,6 +256,96 @@ subroutine tgv_postprocess(ux1, uy1, uz1, ndt)
 
   end subroutine tgv_postprocess
 
+  !
+  ! Register variables for visualization
+  !
+  subroutine tgv_visu_init()
+
+#ifdef ADIOS2
+    use decomp_2d_io, only : adios2_register_variable
+    use visu, only : io_write_real_coarse
+#endif
+
+    implicit none
+
+#ifdef ADIOS2
+    call adios2_register_variable(io_write_real_coarse, "critq", 1, 2)
+#endif
+
+  end subroutine tgv_visu_init
+
+  !
+  ! Output case-specific variables
+  !
+  subroutine tgv_visu(ux1, uy1, uz1, num)
+
+    use decomp_2d, only : ysize, zsize
+    use variables, only : ppy
+    use param, only : half
+    use var, only : ux2,uy2,uz2,ux3,uy3,uz3
+    use var, only : ta1,tb1,tc1,td1,te1,tf1,tg1,th1,ti1
+    use var, only : ta2,tb2,tc2
+    use var, only : ta3,tb3,tc3
+    use x3d_derive
+    use x3d_transpose
+    use x3d_operator_1d
+    use visu, only : write_field
+
+    implicit none
+
+    ! Arguments
+    real(mytype), dimension(xsize(1),xsize(2),xsize(3)), intent(in) :: ux1, uy1, uz1
+    character(len=32), intent(in) :: num
+
+    ! Local variable
+    integer :: i, j, k
+    real(mytype), dimension(xsize(1),xsize(2),xsize(3)) :: array
+
+    ! This is needed to compute derivatives
+    call x3d_transpose_x_to_y(ux1, ux2)
+    call x3d_transpose_x_to_y(uy1, uy2)
+    call x3d_transpose_x_to_y(uz1, uz2)
+    call x3d_transpose_y_to_z(ux2, ux3)
+    call x3d_transpose_y_to_z(uy2, uy3)
+    call x3d_transpose_y_to_z(uz2, uz3)
+
+    ! Compute X derivative
+    call derx (ta1,ux1,x3d_op_derx, xsize(1),xsize(2),xsize(3))
+    call derx (tb1,uy1,x3d_op_derxp,xsize(1),xsize(2),xsize(3))
+    call derx (tc1,uz1,x3d_op_derxp,xsize(1),xsize(2),xsize(3))
+
+    ! Compute Y derivative and transpose back to X
+    call dery (ta2,ux2,x3d_op_deryp,ppy,ysize(1),ysize(2),ysize(3))
+    call dery (tb2,uy2,x3d_op_dery ,ppy,ysize(1),ysize(2),ysize(3))
+    call dery (tc2,uz2,x3d_op_deryp,ppy,ysize(1),ysize(2),ysize(3))
+    call x3d_transpose_y_to_x(ta2, td1)
+    call x3d_transpose_y_to_x(tb2, te1)
+    call x3d_transpose_y_to_x(tc2, tf1)
+
+    ! Compute Z derivative and transpose back to X
+    call derz (ta3,ux3,x3d_op_derzp,zsize(1),zsize(2),zsize(3))
+    call derz (tb3,uy3,x3d_op_derzp,zsize(1),zsize(2),zsize(3))
+    call derz (tc3,uz3,x3d_op_derz ,zsize(1),zsize(2),zsize(3))
+    call x3d_transpose_z_to_y(ta3, ta2)
+    call x3d_transpose_z_to_y(tb3, tb2)
+    call x3d_transpose_z_to_y(tc3, tc2)
+    call x3d_transpose_y_to_x(ta2, tg1)
+    call x3d_transpose_y_to_x(tb2, th1)
+    call x3d_transpose_y_to_x(tc2, ti1)
+    !du/dx=ta1   du/dy=td1   du/dz=tg1
+    !dv/dx=tb1   dv/dy=te1   dv/dz=th1
+    !dw/dx=tc1   dw/dy=tf1   dw/dz=ti1
+
+    ! Compute Q criterion and write it
+    do concurrent (k=1:xsize(3), j=1:xsize(2), i=1:xsize(1))
+      array(i,j,k) = - half*(ta1(i,j,k)**2+te1(i,j,k)**2+ti1(i,j,k)**2) &
+                     - td1(i,j,k)*tb1(i,j,k)&
+                     - tg1(i,j,k)*tc1(i,j,k)&
+                     - th1(i,j,k)*tf1(i,j,k)
+    enddo
+    call write_field(array, ".", "critq", trim(num))
+
+  end subroutine tgv_visu
   !
   ! Finalize case-specific IO
   !
diff --git a/src/case.f90 b/src/case.f90
index d2d1049..c1f3c11 100644
--- a/src/case.f90
+++ b/src/case.f90
@@ -18,6 +18,7 @@ module case
             case_init, &
             case_bc, &
             case_forcing, &
+            case_visu_init, &
             case_visu, &
             case_postprocess, &
             case_finalize
@@ -94,36 +95,55 @@ subroutine case_forcing(dux1, duy1, duz1)
 
   end subroutine case_forcing
 
+  !
+  ! Register variables for visualization
+  !
+  subroutine case_visu_init()
+
+    implicit none
+
+    if (itype == itype_tgv)  call tgv_visu_init()
+
+  end subroutine case_visu_init
+
   !
   ! Visualization
   ! This is called when itime % ioutput = 0
   !
-  subroutine case_visu()
+  subroutine case_visu(ux1, uy1, uz1, num)
 
     implicit none
 
+    ! Arguments
+    real(mytype), dimension(xsize(1),xsize(2),xsize(3)), intent(in) :: ux1, uy1, uz1
+    character(len=32), intent(in) :: num
+
+    if (itype == itype_tgv) call tgv_visu(ux1, uy1, uz1, num)
+
   end subroutine case_visu
 
   !
   ! Case-specific post-processing
   ! This is called at the end of each time step
   !
-  subroutine case_postprocess(ux1, uy1, uz1, ndt)
+  subroutine case_postprocess(ux1, uy1, uz1, it, num)
 
-    use param, only : ivisu, ioutput
+    use param, only : ioutput
+    use variables, only : nvisu
 
     implicit none
 
     ! Arguments
     real(mytype), dimension(xsize(1),xsize(2),xsize(3)), intent(in) :: ux1, uy1, uz1
-    integer, intent(in) :: ndt
+    integer, intent(in) :: it
+    character(len=32), intent(in) :: num
 
-    if (itype == itype_tgv) call tgv_postprocess(ux1, uy1, uz1, ndt)
+    if (itype == itype_tgv) call tgv_postprocess(ux1, uy1, uz1, it)
 
-    if (itype == itype_tgv2d) call tgv2d_postprocess(ux1, uy1, uz1, ndt)
+    if (itype == itype_tgv2d) call tgv2d_postprocess(ux1, uy1, uz1, it)
 
-    if ((ivisu /= 0).and.(ioutput /= 0)) then
-      if (mod(ndt, ioutput) == 0) call case_visu()
+    if ((nvisu /= 0).and.(ioutput /= 0)) then
+      if (mod(it, ioutput) == 0) call case_visu(ux1, uy1, uz1, num)
     endif
 
   end subroutine case_postprocess
diff --git a/src/visu.f90 b/src/visu.f90
new file mode 100644
index 0000000..70d366d
--- /dev/null
+++ b/src/visu.f90
@@ -0,0 +1,612 @@
+!Copyright (c) 2012-2022, Xcompact3d
+!This file is part of Xcompact3d (xcompact3d.com)
+!SPDX-License-Identifier: BSD 3-Clause
+
+module visu
+
+#ifdef ADIOS2
+  use adios2
+#endif
+  
+  implicit none
+
+  ! True to activate the XDMF output
+  logical, save :: use_xdmf = .true.
+  ! True to use the new enumeration
+  logical, save :: filenamedigits = .false.
+  ! output2D is defined in the input.i3d file
+  !        0 for 3D output (default)
+  !        1 for 2D output with X average
+  !        2 for 2D output with Y average
+  !        3 for 2D output with Z average
+  integer, save :: output2D
+  integer :: ioxdmf
+  character(len=9) :: ifilenameformat = '(I3.3)'
+  real, save :: tstart, tend
+#ifdef ADIOS2
+  type(adios2_adios) :: adios
+  type(adios2_io) :: io_write_real_coarse
+  type(adios2_engine) :: engine_write_real_coarse
+#endif
+
+  private
+  public :: output2D, visu_init, visu_finalise, write_snapshot, end_snapshot, write_field
+
+contains
+
+  !
+  ! Initialize the visu module
+  !
+  subroutine visu_init()
+
+    use MPI
+    use param, only : iscalar, ilast, ifirst, ioutput, istret
+    use variables, only : numscalar, prec, nvisu
+    use decomp_2d, only : nrank, mytype, xszV, yszV, zszV
+#ifdef ADIOS2
+    use decomp_2d_io, only : adios2_register_variable
+#endif
+    
+    implicit none
+
+    ! Local variables
+    logical :: dir_exists
+    integer :: noutput, nsnapout
+    real(mytype) :: memout
+
+#ifdef ADIOS2
+    integer :: ierror
+    logical :: adios2_debug_mode
+    character(len=80) :: config_file="adios2_config.xml"
+    character(len=80) :: outfile
+    integer :: is
+#endif
+
+#ifdef DEBUG
+    if (nrank == 0) then
+       print *, "# visu_init start"
+    endif
+#endif
+    
+    ! Create folder if needed
+    ! XXX: Is this needed for ADIOS2?
+    if (nrank==0) then
+      inquire(file="data", exist=dir_exists)
+      if (.not.dir_exists) then
+        call system("mkdir data 2> /dev/null")
+      end if
+    end if
+
+    ! ! HDD usage of visu module
+    ! if (nrank==0) then
+    !   noutput = (ilast - ifirst +1)/ioutput
+
+    !   nsnapout = 4
+    !   if (iscalar.ne.0) nsnapout = nsnapout + numscalar
+
+    !   memout = prec * nsnapout * noutput
+    !   if (output2D.eq.0) then
+    !     memout = memout * xszV(1) * yszV(2) * zszV(3)
+    !   else if (output2D.eq.1) then
+    !     memout = memout *           yszV(2) * zszV(3)
+    !   else if (output2D.eq.2) then
+    !     memout = memout * xszV(1)           * zszV(3)
+    !   else if (output2D.eq.3) then
+    !     memout = memout * xszV(1) * yszV(2)
+    !   endif
+    !   write(*,*)'==========================================================='
+    !   write(*,*)'Visu module requires ',real(memout*1e-9,4),'GB'
+    !   write(*,*)'==========================================================='
+    ! end if
+
+    ! Safety check
+    if (output2D < 0 .or. output2D > 3 &
+        .or. (output2d == 2.and.istret /= 0)) then
+      if (nrank.eq.0) write(*,*) "Visu module: incorrect value for output2D."
+      call MPI_ABORT(MPI_COMM_WORLD, 0, noutput)
+      stop
+    endif
+
+#ifdef ADIOS2
+    !! TODO: make this a runtime-option
+    adios2_debug_mode = .true.
+
+    call adios2_init(adios, trim(config_file), MPI_COMM_WORLD, adios2_debug_mode, ierror)
+    if (ierror.ne.0) then
+       print *, "Error initialising ADIOS2 - is adios2_config.xml present and valid?"
+       call MPI_ABORT(MPI_COMM_WORLD, -1, ierror)
+    endif
+    call adios2_declare_io(io_write_real_coarse, adios, "solution-io", ierror)
+    if (io_write_real_coarse % engine_type.eq."BP4") then
+       write(outfile, *) "data.bp4"
+    else if (io_write_real_coarse % engine_type.eq."HDF5") then
+       write(outfile, *) "data.hdf5"
+    else
+       print *, "Unknown engine!"
+       call MPI_ABORT(MPI_COMM_WORLD, -1, ierror)
+    endif
+
+    !! Register variables
+    call adios2_register_variable(io_write_real_coarse, "ux", 1, 2)
+    call adios2_register_variable(io_write_real_coarse, "uy", 1, 2)
+    call adios2_register_variable(io_write_real_coarse, "uz", 1, 2)
+    call adios2_register_variable(io_write_real_coarse, "pp", 1, 2)
+    ! if (ilmn) then
+    !    call adios2_register_variable(io_write_real_coarse, "rho", 1, 2)
+    ! endif
+    if (iscalar.ne.0) then
+       do is = 1, numscalar
+          call adios2_register_variable(io_write_real_coarse, "phi"//char(48+is), 1, 2)
+       enddo
+    endif
+    
+    call adios2_open(engine_write_real_coarse, io_write_real_coarse, trim(outfile), adios2_mode_write, ierror)
+#endif
+
+#ifdef DEBUG
+    if (nrank == 0) then
+       print *, "# visu_init end"
+    endif
+#endif
+
+  end subroutine visu_init
+
+  !
+  ! Finalise the visu module
+  ! - Currently only needed to clean up ADIOS2
+  !
+  subroutine visu_finalise()
+    
+#ifdef ADIOS2
+    use adios2
+#endif
+  implicit none
+
+#ifdef ADIOS2
+  integer :: ierr
+#endif
+  
+#ifdef ADIOS2
+    call adios2_close(engine_write_real_coarse, ierr)
+    call adios2_finalize(adios, ierr)
+#endif
+    
+  end subroutine visu_finalise
+
+  !
+  ! Write a snapshot
+  !
+  ! subroutine write_snapshot(rho1, ux1, uy1, uz1, pp3, phi1, ep1, itime, num)
+  subroutine write_snapshot(ux1, uy1, uz1, itime, num)
+
+    use decomp_2d, only : mytype, xsize
+    use decomp_2d, only : nrank
+
+    use param, only : iscalar, ioutput
+
+    ! use variables, only : sx, cifip6, cisip6, ciwip6, cifx6, cisx6, ciwx6
+    ! use variables, only : sy, cifip6y, cisip6y, ciwip6y, cify6, cisy6, ciwy6
+    ! use variables, only : sz, cifip6z, cisip6z, ciwip6z, cifz6, cisz6, ciwz6
+    use variables, only : numscalar
+
+    ! use var, only : pp1, ta1, di1, nxmsize
+    ! use var, only : pp2, ppi2, dip2, ph2, nymsize
+    ! use var, only : ppi3, dip3, ph3, nzmsize
+    ! use var, only : npress
+
+    ! use tools, only : rescale_pressure
+
+    implicit none
+
+    !! inputs
+    real(mytype), dimension(xsize(1), xsize(2), xsize(3)), intent(in) :: ux1, uy1, uz1
+    ! real(mytype), dimension(xsize(1), xsize(2), xsize(3)), intent(in) :: ep1
+    ! real(mytype), dimension(xsize(1), xsize(2), xsize(3), nrhotime), intent(in) :: rho1
+    ! real(mytype), dimension(ph3%zst(1):ph3%zen(1),ph3%zst(2):ph3%zen(2),nzmsize,npress), intent(in) :: pp3
+    ! real(mytype), dimension(xsize(1), xsize(2), xsize(3), numscalar), intent(in) :: phi1
+    integer, intent(in) :: itime
+    character(len=32), intent(out) :: num
+
+    ! Local variables
+    integer :: is
+    integer :: ierr
+    character(len=30) :: scname
+
+#ifdef DEBUG
+    if (nrank == 0) then
+       print *, "# write_snapshot start"
+    endif
+#endif
+    
+    ! Update log file
+    if (nrank.eq.0) then
+      call cpu_time(tstart)
+      print *,'Writing snapshots =>',itime/ioutput
+    end if
+#ifdef ADIOS2
+    call adios2_begin_step(engine_write_real_coarse, adios2_step_mode_append, ierr)
+#endif
+
+    ! Snapshot number
+#ifndef ADIOS2
+    if (filenamedigits) then
+       ! New enumeration system, it works integrated with xcompact3d_toolbox
+       write(num, ifilenameformat) itime
+    else
+       ! Classic enumeration system
+       write(num, ifilenameformat) itime/ioutput
+    endif
+#else
+    ! ADIOS2 is zero-indexed
+    write(num, '(I0)') itime/ioutput - 1
+#endif
+    
+    ! Write XDMF header
+    if (use_xdmf) call write_xdmf_header(".", "snapshot", trim(num))
+
+    ! Write velocity
+    call write_field(ux1, ".", "ux", trim(num))
+    call write_field(uy1, ".", "uy", trim(num))
+    call write_field(uz1, ".", "uz", trim(num))
+
+    ! ! Interpolate pressure
+    ! !WORK Z-PENCILS
+    ! call interzpv(ppi3,pp3(:,:,:,1),dip3,sz,cifip6z,cisip6z,ciwip6z,cifz6,cisz6,ciwz6,&
+    ! (ph3%zen(1)-ph3%zst(1)+1),(ph3%zen(2)-ph3%zst(2)+1),nzmsize,zsize(3),1)
+    ! !WORK Y-PENCILS
+    ! call transpose_z_to_y(ppi3,pp2,ph3) !nxm nym nz
+    ! call interypv(ppi2,pp2,dip2,sy,cifip6y,cisip6y,ciwip6y,cify6,cisy6,ciwy6,&
+    !         (ph3%yen(1)-ph3%yst(1)+1),nymsize,ysize(2),ysize(3),1)
+    ! !WORK X-PENCILS
+    ! call transpose_y_to_x(ppi2,pp1,ph2) !nxm ny nz
+    ! call interxpv(ta1,pp1,di1,sx,cifip6,cisip6,ciwip6,cifx6,cisx6,ciwx6,&
+    !         nxmsize,xsize(1),xsize(2),xsize(3),1)
+
+
+    ! ! Rescale pressure
+    ! call rescale_pressure(ta1)
+
+    ! ! Write pressure
+    ! call write_field(ta1, ".", "pp", trim(num), .true.)
+
+    ! ! LMN - write density
+    ! if (ilmn) call write_field(rho1(:,:,:,1), ".", "rho", trim(num))
+
+    ! ! Write scalars
+    ! if (iscalar.ne.0) then
+    !   do is = 1, numscalar
+    !     write(scname,"('phi',I2.2)") is
+    !     call write_field(phi1(:,:,:,is), ".", trim(scname), trim(num), .true.)
+    !   enddo
+    ! endif
+
+#ifdef DEBUG
+    if (nrank == 0) then
+       print *, "# write_snapshot end"
+    endif
+#endif
+
+  end subroutine write_snapshot
+
+  subroutine end_snapshot(itime, num)
+
+    use decomp_2d, only : nrank
+    use param, only : istret, xlx, yly, zlz
+    use variables, only : nx, ny, nz, beta
+
+    implicit none
+
+    integer, intent(in) :: itime
+    character(len=32), intent(in) :: num
+
+    character(len=32) :: fmt2, fmt3, fmt4
+    integer :: is
+    integer :: ierr
+
+#ifdef DEBUG
+    if (nrank == 0) then
+       print *, "# end_snapshot start"
+    endif
+#endif
+    
+    ! Write XDMF footer
+    if (use_xdmf) call write_xdmf_footer()
+
+    ! Add metadata if XDMF is not used
+    if (.not. use_xdmf) then
+      if (nrank==0) then
+
+        write(fmt2,'("(A,I16)")')
+        write(fmt3,'("(A,F16.4)")')
+        write(fmt4,'("(A,F16.12)")')
+
+        open(newunit=is,file="./data/snap"//trim(num)//".ini",action='write',status='replace')
+        write(is,'(A)')'[domain]'
+        write(is,fmt2) 'nx=      ',nx
+        write(is,fmt2) 'ny=      ',ny
+        write(is,fmt2) 'nz=      ',nz
+        write(is,fmt2) 'istret=  ',istret
+        write(is,fmt4) 'beta=    ',beta
+        write(is,fmt3) 'Lx=      ',xlx
+        write(is,fmt3) 'Ly=      ',yly
+        write(is,fmt3) 'Lz=      ',zlz
+        write(is,'(A)')'[time]'
+        write(is,fmt2) 'itime=   ',itime
+        ! write(is,fmt3) 'dt=      ',dt
+        ! write(is,fmt3) 't =      ',t
+        close(is)
+
+      endif
+    endif
+
+#ifdef ADIOS2
+    call adios2_end_step(engine_write_real_coarse, ierr)
+#endif
+
+    ! Update log file
+    if (nrank.eq.0) then
+      call cpu_time(tend)
+      write(*,'(" Time for writing snapshots (s): ",F12.8)') tend-tstart
+    endif
+
+#ifdef DEBUG
+    if (nrank == 0) then
+       print *, "# end_snapshot end"
+    endif
+#endif
+
+  end subroutine end_snapshot
+
+  !
+  ! Write the header of the XDMF file
+  ! Adapted from https://github.com/fschuch/Xcompact3d/blob/master/src/visu.f90
+  !
+  subroutine write_xdmf_header(pathname, filename, num)
+
+    use variables, only : nvisu, yp
+    use param, only : dx,dy,dz,istret
+    use decomp_2d, only : mytype, nrank, xszV, yszV, zszV, ystV
+
+    implicit none
+
+    ! Arguments
+    character(len=*), intent(in) :: pathname, filename, num
+
+    ! Local variables
+    integer :: i,k
+    real(mytype) :: xp(xszV(1)), zp(zszV(3))
+
+#ifdef DEBUG
+    if (nrank == 0) then
+       print *, "## write_xdmf_header start"
+    end if
+#endif
+    
+    if (nrank.eq.0) then
+      OPEN(newunit=ioxdmf,file="./data/"//pathname//"/"//filename//'-'//num//'.xdmf')
+
+      write(ioxdmf,'(A)')'<?xml version="1.0" ?>'
+      write(ioxdmf,*)'<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>'
+      write(ioxdmf,*)'<Xdmf xmlns:xi="http://www.w3.org/2001/XInclude" Version="2.0">'
+      write(ioxdmf,*)'<Domain>'
+      if (istret.ne.0) then
+        do i=1,xszV(1)
+          xp(i) = real(i-1,mytype)*dx*nvisu
+        enddo
+        do k=1,zszV(3)
+          zp(k) = real(k-1,mytype)*dz*nvisu
+        enddo
+        write(ioxdmf,*)'    <Topology name="topo" TopologyType="3DRectMesh"'
+        if (output2D.eq.0) then
+          write(ioxdmf,*)'        Dimensions="',zszV(3),yszV(2),xszV(1),'">'
+        else if (output2D.eq.1) then
+          write(ioxdmf,*)'        Dimensions="',zszV(3),yszV(2),1,'">'
+        else if (output2D.eq.2) then
+          write(ioxdmf,*)'        Dimensions="',zszV(3),1,xszV(1),'">'
+        else if (output2D.eq.3) then
+          write(ioxdmf,*)'        Dimensions="',1,yszV(2),xszV(1),'">'
+        endif
+        write(ioxdmf,*)'    </Topology>'
+        write(ioxdmf,*)'    <Geometry name="geo" Type="VXVYVZ">'
+        if (output2D.ne.1) then
+          write(ioxdmf,*)'        <DataItem Dimensions="',xszV(1),'" NumberType="Float" Precision="4" Format="XML">'
+          write(ioxdmf,*)'        ',xp(:)
+        else
+          write(ioxdmf,*)'        <DataItem Dimensions="1" NumberType="Float" Precision="4" Format="XML">'
+          write(ioxdmf,*)'        ',xp(1)
+        endif
+        write(ioxdmf,*)'        </DataItem>'
+        if (output2D.ne.2) then
+          write(ioxdmf,*)'        <DataItem Dimensions="',yszV(2),'" NumberType="Float" Precision="4" Format="XML">'
+          write(ioxdmf,*)'        ',yp(ystV(1)::nvisu)
+        else
+          write(ioxdmf,*)'        <DataItem Dimensions="1" NumberType="Float" Precision="4" Format="XML">'
+          write(ioxdmf,*)'        ',yp(1)
+        endif
+        write(ioxdmf,*)'        </DataItem>'
+        if (output2D.ne.3) then
+          write(ioxdmf,*)'        <DataItem Dimensions="',zszV(3),'" NumberType="Float" Precision="4" Format="XML">'
+          write(ioxdmf,*)'        ',zp(:)
+        else
+          write(ioxdmf,*)'        <DataItem Dimensions="1" NumberType="Float" Precision="4" Format="XML">'
+          write(ioxdmf,*)'        ',zp(1)
+        endif
+        write(ioxdmf,*)'        </DataItem>'
+        write(ioxdmf,*)'    </Geometry>'
+      else
+        write(ioxdmf,*)'    <Topology name="topo" TopologyType="3DCoRectMesh"'
+        if (output2D.eq.0) then
+          write(ioxdmf,*)'        Dimensions="',zszV(3),yszV(2),xszV(1),'">'
+        else if (output2D.eq.1) then
+          write(ioxdmf,*)'        Dimensions="',zszV(3),yszV(2),1,'">'
+        else if (output2D.eq.2) then
+          write(ioxdmf,*)'        Dimensions="',zszV(3),1,xszV(1),'">'
+        else if (output2D.eq.3) then
+          write(ioxdmf,*)'        Dimensions="',1,yszV(2),xszV(1),'">'
+        endif
+        write(ioxdmf,*)'    </Topology>'
+        write(ioxdmf,*)'    <Geometry name="geo" Type="ORIGIN_DXDYDZ">'
+        write(ioxdmf,*)'        <!-- Origin -->'
+        write(ioxdmf,*)'        <DataItem Format="XML" Dimensions="3">'
+        write(ioxdmf,*)'        0.0 0.0 0.0'
+        write(ioxdmf,*)'        </DataItem>'
+        write(ioxdmf,*)'        <!-- DxDyDz -->'
+        write(ioxdmf,*)'        <DataItem Format="XML" Dimensions="3">'
+        if (output2D.eq.0) then
+          write(ioxdmf,*)'        ',nvisu*dz,nvisu*dy,nvisu*dx
+        else if (output2D.eq.1) then
+          write(ioxdmf,*)'        ',dz,dy,1.
+        else if (output2D.eq.2) then
+          write(ioxdmf,*)'        ',dz,1.,dx
+        else if (output2D.eq.3) then
+          write(ioxdmf,*)'        ',1.,dy,dx
+        endif
+        write(ioxdmf,*)'        </DataItem>'
+        write(ioxdmf,*)'    </Geometry>'
+      endif
+      write(ioxdmf,*)'    <Grid Name="'//num//'" GridType="Uniform">'
+      write(ioxdmf,*)'        <Topology Reference="/Xdmf/Domain/Topology[1]"/>'
+      write(ioxdmf,*)'        <Geometry Reference="/Xdmf/Domain/Geometry[1]"/>'
+    endif
+
+#ifdef DEBUG
+    if (nrank == 0) then
+       print *, "## write_xdmf_header end"
+    end if
+#endif
+  end subroutine write_xdmf_header
+
+  !
+  ! Write the footer of the XDMF file
+  ! Adapted from https://github.com/fschuch/Xcompact3d/blob/master/src/visu.f90
+  !
+  subroutine write_xdmf_footer()
+
+    use decomp_2d, only : nrank
+
+    implicit none
+
+#ifdef DEBUG
+    if (nrank == 0) then
+       print *, "## write_xdmf_footer start"
+    endif
+#endif
+    
+    if (nrank.eq.0) then
+      write(ioxdmf,'(/)')
+      write(ioxdmf,*) '    </Grid>'
+      write(ioxdmf,*) '</Domain>'
+      write(ioxdmf,'(A)') '</Xdmf>'
+      close(ioxdmf)
+    endif
+
+#ifdef DEBUG
+    if (nrank == 0) then
+       print *, "## write_xdmf_footer end"
+    endif
+#endif
+
+  end subroutine write_xdmf_footer
+
+  !
+  ! Write the given field for visualization
+  ! Adapted from https://github.com/fschuch/Xcompact3d/blob/master/src/visu.f90
+  !
+  subroutine write_field(f1, pathname, filename, num, skip_ibm)
+
+    use mpi
+    
+    use param, only : iibm
+    use decomp_2d, only : mytype, xsize, xszV, yszV, zszV
+    use decomp_2d, only : nrank, fine_to_coarseV
+    use decomp_2d_io, only : decomp_2d_write_one, decomp_2d_write_plane
+
+    implicit none
+
+    real(mytype), intent(in), dimension(xsize(1),xsize(2),xsize(3)) :: f1
+    character(len=*), intent(in) :: pathname, filename, num 
+    logical, optional, intent(in) :: skip_ibm
+
+    real(mytype), dimension(xsize(1),xsize(2),xsize(3)) :: uvisu
+    real(mytype), dimension(xsize(1),xsize(2),xsize(3)) :: local_array
+
+    integer :: ierr
+
+#ifdef DEBUG
+    if (nrank == 0) then
+       print *, "## write_field -- ", filename, " -- start"
+    endif
+#endif
+
+    if (use_xdmf) then
+       if (nrank.eq.0) then
+          write(ioxdmf,*)'        <Attribute Name="'//filename//'" Center="Node">'
+#ifndef ADIOS2
+          write(ioxdmf,*)'           <DataItem Format="Binary"'
+#else
+          write(ioxdmf,*)'           <DataItem Format="HDF"'
+#endif
+#ifdef DOUBLE_PREC
+#ifdef SAVE_SINGLE
+          if (output2D.eq.0) then
+             write(ioxdmf,*)'            DataType="Float" Precision="4" Endian="little" Seek="0"'
+          else
+             write(ioxdmf,*)'            DataType="Float" Precision="8" Endian="little" Seek="0"'
+          endif
+#else
+          write(ioxdmf,*)'            DataType="Float" Precision="8" Endian="little" Seek="0"'
+#endif
+#else
+          write(ioxdmf,*)'            DataType="Float" Precision="4" Endian="little" Seek="0"'
+#endif
+          if (output2D.eq.0) then
+             write(ioxdmf,*)'            Dimensions="',zszV(3),yszV(2),xszV(1),'">'
+          else if (output2D.eq.1) then
+             write(ioxdmf,*)'            Dimensions="',zszV(3),yszV(2),1,'">'
+          else if (output2D.eq.2) then
+             write(ioxdmf,*)'            Dimensions="',zszV(3),1,xszV(1),'">'
+          else if (output2D.eq.3) then
+             write(ioxdmf,*)'            Dimensions="',1,yszV(2),xszV(1),'">'
+          endif
+#ifndef ADIOS2
+          write(ioxdmf,*)'              ./'//pathname//"/"//filename//'-'//num//'.bin'
+#else
+          write(ioxdmf,*)'              ../data.hdf5:/Step'//num//'/'//filename
+#endif
+          write(ioxdmf,*)'           </DataItem>'
+          write(ioxdmf,*)'        </Attribute>'
+       endif
+    endif
+    
+    if (iibm==2 .and. .not.present(skip_ibm)) then
+       ! local_array(:,:,:) = (one - ep1(:,:,:)) * f1(:,:,:)
+       local_array(:,:,:) = f1(:,:,:)
+    else
+       local_array(:,:,:) = f1(:,:,:)
+    endif
+
+    if (output2D.eq.0) then
+#ifndef ADIOS2
+       uvisu = 0.0_mytype
+       call fine_to_coarseV(1,local_array,uvisu)
+       call decomp_2d_write_one(1,uvisu,"./data/"//pathname//'/'//filename//'-'//num//'.bin',2)
+#else
+       if (iibm==2 .and. (.not.present(skip_ibm))) then
+          print *, "Not Implemented: currently ADIOS2 IO doesn't support IBM-blanking"
+          call MPI_ABORT(MPI_COMM_WORLD, -1, ierr)
+       endif
+       call decomp_2d_write_one(1,f1,filename,2,adios,engine_write_real_coarse,io_write_real_coarse)
+#endif
+    else
+       call decomp_2d_write_plane(1,local_array,output2D,-1,"./data/"//pathname//'/'//filename//'-'//num//'.bin')
+    endif
+
+#ifdef DEBUG
+    if (nrank == 0) then
+       print *, "## write_field -- ", filename, " -- end"
+    endif
+#endif
+
+  end subroutine write_field
+
+end module visu
diff --git a/src/x3d_tools.f90 b/src/x3d_tools.f90
index edcaadb..cf1cc42 100644
--- a/src/x3d_tools.f90
+++ b/src/x3d_tools.f90
@@ -91,6 +91,8 @@ subroutine init_xcompact3d(ndt_max)
 
   use var
 
+  use visu, only : visu_init
+  
   use variables, only : nx, ny, nz, nxm, nym, nzm
   use variables, only : p_row, p_col
   use variables, only : test_mode
@@ -180,6 +182,9 @@ subroutine init_xcompact3d(ndt_max)
   call decomp_2d_poisson_init()
   call decomp_info_init(nxm,nym,nzm,phG)
 
+  call visu_init()
+  call case_visu_init()
+  
 endsubroutine init_xcompact3d
 !########################################################################
 !########################################################################
@@ -196,6 +201,8 @@ subroutine finalise_xcompact3d(flag)
   use x3d_derive, only : x3d_derive_finalize
   use var, only : var_finalize
 
+  use visu, only : visu_finalise
+  
   implicit none
 
   logical, intent(in) :: flag
@@ -214,9 +221,11 @@ subroutine finalise_xcompact3d(flag)
   call x3d_operator_z_data_finalize()
   call var_finalize()
 
+  call visu_finalise()
+
   call decomp_2d_finalize()
   if (flag) then
     CALL MPI_FINALIZE(ierr)
   endif
-
+  
 endsubroutine finalise_xcompact3d
diff --git a/src/xcompact3d.f90 b/src/xcompact3d.f90
index 031fa64..2858700 100644
--- a/src/xcompact3d.f90
+++ b/src/xcompact3d.f90
@@ -8,12 +8,13 @@ program xcompact3d
 
   use var
   use decomp_2d, only : nrank, xsize, real_type, decomp_2d_warning
-  use param,   only : dt, zero, itr
+  use param,   only : dt, zero, itr, itime, ioutput
   use transeq, only : calculate_transeq_rhs
   use navier,  only : solve_poisson, cor_vel
   use case
   use time_integrators, only : int_time
-
+  use visu, only : write_snapshot, end_snapshot
+  
   implicit none
 
   double precision :: tstart, tend, telapsed, tmin, tmax
@@ -21,6 +22,7 @@ program xcompact3d
   integer :: i, j, k
   integer :: ndt, ndt_max
   integer :: code
+  character(len=32) :: num
 
   call boot_xcompact3d()
 
@@ -30,8 +32,14 @@ program xcompact3d
 
   telapsed = 0
   tmin = telapsed
-  ndt = 1
-
+  ndt = 0
+
+  ioutput = 5
+  print *, "=== WARNING ==="
+  print *, " Writing with ioutput=",ioutput,", edit xcompact3d.f90 and recompile to change"
+  print *, "=== WARNING ==="
+  
+  itime = 0
   do while(ndt < ndt_max)
 
      itr = 1 ! no inner iterations
@@ -63,8 +71,12 @@ program xcompact3d
         print *, "Elapse time min ", tmin, " max ", tmax
      end if
 
+     !! End of timestep
+     itime = itime + 1
      ndt = ndt + 1
-     call case_postprocess(ux1, uy1, uz1, ndt)
+     if (mod(itime,ioutput) == 0) call write_snapshot(ux1, uy1, uz1, itime, num)
+     call case_postprocess(ux1, uy1, uz1, itime, num)
+     if (mod(itime,ioutput) == 0) call end_snapshot(itime, num)
 
   end do