! $Id$
!
!  This module is a dummy module for non-HDF5 IO.
!  28-Oct-2016/PABoudin: coded
!
module HDF5_IO
!
  use Cdata
  use General, only: keep_compiler_quiet, itoa, numeric_precision
  use Messages, only: fatal_error
  use Mpicomm, only: lroot
!
  implicit none
!
  interface hdf5_input_slice
    module procedure input_slice_real_arr
    module procedure input_slice_scat_arr
  endinterface
!
  interface input_hdf5
    module procedure input_hdf5_int_0D
    module procedure input_hdf5_int_1D
    module procedure input_hdf5_0D
    module procedure input_hdf5_1D
    module procedure input_hdf5_part_2D
    module procedure input_hdf5_profile_1D
    module procedure input_hdf5_3D
    module procedure input_hdf5_4D
  endinterface
!
  interface output_hdf5
    module procedure output_hdf5_string
    module procedure output_hdf5_int_0D
    module procedure output_hdf5_int_1D
    module procedure output_hdf5_0D
    module procedure output_hdf5_1D
    module procedure output_hdf5_part_2D
    module procedure output_hdf5_profile_1D
    module procedure output_local_hdf5_2D
    module procedure output_hdf5_slice_2D
    module procedure output_local_hdf5_3D
    module procedure output_hdf5_3D
    module procedure output_hdf5_4D
  endinterface
!
  interface output_hdf5_double
    module procedure output_hdf5_double_0D
    module procedure output_hdf5_double_1D
  endinterface
!
  include 'hdf5_io.h'
!
  private
!
  integer, parameter :: lun_input = 89, lun_output = 92
!
  contains
!***********************************************************************
    subroutine initialize_hdf5
!
      ! nothing to do
!
    endsubroutine initialize_hdf5
!***********************************************************************
    subroutine finalize_hdf5
!
      call fatal_error ('finalize_hdf5', 'You can not use HDF5 without setting an HDF5_IO module.')
!
    endsubroutine finalize_hdf5
!***********************************************************************
    subroutine file_open_hdf5(file, truncate, global, read_only, write, comm)
!
      character (len=*), intent(in) :: file
      logical, optional, intent(in) :: truncate
      logical, optional, intent(in) :: global
      logical, optional, intent(in) :: read_only
      logical, optional, intent(in) :: write
      integer, optional, intent(in) :: comm
!
      call keep_compiler_quiet(file)
      call keep_compiler_quiet(.true.,truncate, global, read_only)
      call keep_compiler_quiet(.true.,write)
      call keep_compiler_quiet(comm)

      call fatal_error ('file_open_hdf5', 'You can not use HDF5 without setting an HDF5_IO module.')
!
    endsubroutine file_open_hdf5
!***********************************************************************
    subroutine file_close_hdf5
!
      call fatal_error ('file_close_hdf5', 'You can not use HDF5 without setting an HDF5_IO module.')
!
    endsubroutine file_close_hdf5
!***********************************************************************
    subroutine create_group_hdf5(name)
!
      character (len=*), intent(in) :: name
!
      call fatal_error ('create_group_hdf5', 'You can not use HDF5 without setting an HDF5_IO module.')
      call keep_compiler_quiet(name)
!
    endsubroutine create_group_hdf5
!***********************************************************************
    logical function exists_in_hdf5(name)
!
      character (len=*), intent(in) :: name
!
      call fatal_error ('exists_in_hdf5', 'You can not use HDF5 without setting an HDF5_IO module.')
      call keep_compiler_quiet(name)
      exists_in_hdf5 = .false.
!
    endfunction exists_in_hdf5
!***********************************************************************
    subroutine input_hdf5_int_0D(name, data)
!
      character (len=*), intent(in) :: name
      integer, intent(out) :: data
!
      call fatal_error ('input_hdf5_int_0D', 'You can not use HDF5 without setting an HDF5_IO module.')
      call keep_compiler_quiet(name)
      data = -1
!
    endsubroutine input_hdf5_int_0D
!***********************************************************************
    subroutine input_hdf5_int_1D(name, data, nv, same_size)
!
      character (len=*), intent(in) :: name
      integer, intent(in) :: nv
      integer, dimension (nv), intent(out) :: data
      logical, optional, intent(in) :: same_size
!
      call fatal_error ('input_hdf5_int_1D', 'You can not use HDF5 without setting an HDF5_IO module.')
      call keep_compiler_quiet(name)
      data(:) = -1
      call keep_compiler_quiet(.true., same_size)
!
    endsubroutine input_hdf5_int_1D
!***********************************************************************
    subroutine input_hdf5_0D(name, data)
!
      character (len=*), intent(in) :: name
      real, intent(out) :: data
!
      call fatal_error ('input_hdf5_0D', 'You can not use HDF5 without setting an HDF5_IO module.')
      call keep_compiler_quiet(name)
      data = -1.0
!
    endsubroutine input_hdf5_0D
!***********************************************************************
    subroutine input_hdf5_1D(name, data, nv, same_size)
!
      character (len=*), intent(in) :: name
      integer, intent(in) :: nv
      real, dimension (nv), intent(out) :: data
      logical, optional, intent(in) :: same_size
!
      call fatal_error ('input_hdf5_1D', 'You can not use HDF5 without setting an HDF5_IO module.')
      call keep_compiler_quiet(name)
      data(:) = -1.0
      call keep_compiler_quiet(.true., same_size)
!
    endsubroutine input_hdf5_1D
!***********************************************************************
    subroutine input_hdf5_part_2D(name, data, mv, nc, nv)
!
      character (len=*), intent(in) :: name
      integer, intent(in) :: mv, nc
      real, dimension (mv,mparray), intent(out) :: data
      integer, intent(out) :: nv
!
      call fatal_error ('input_hdf5_part_2D', 'You can not use HDF5 without setting an HDF5_IO module.')
      call keep_compiler_quiet(name)
      call keep_compiler_quiet(mv)
      call keep_compiler_quiet(nc)
      data(:,:) = -1.0
      nv = -1
!
    endsubroutine input_hdf5_part_2D
!***********************************************************************
    subroutine input_hdf5_profile_1D(name, data, ldim, gdim, np1, np2)
!
      character (len=*), intent(in) :: name
      real, dimension (:) :: data
      integer, intent(in) :: ldim, gdim, np1, np2
!
      call fatal_error ('input_hdf5_profile_1D', 'You can not use HDF5 without setting an HDF5_IO module.')
      call keep_compiler_quiet(name)
      call keep_compiler_quiet(ldim)
      call keep_compiler_quiet(gdim)
      call keep_compiler_quiet(np1)
      call keep_compiler_quiet(np2)
!
    endsubroutine input_hdf5_profile_1D
!***********************************************************************
    subroutine input_hdf5_3D(name, data)
!
      character (len=*), intent(in) :: name
      real, dimension (mx,my,mz), intent(out) :: data
!
      call fatal_error ('input_hdf5_3D', 'You can not use HDF5 without setting an HDF5_IO module.')
      call keep_compiler_quiet(name)
      data(:,:,:) = -1.0
!
    endsubroutine input_hdf5_3D
!***********************************************************************
    subroutine input_hdf5_4D(name, data, nv)
!
      character (len=*), intent(in) :: name
      integer, intent(in) :: nv
      real, dimension (mx,my,mz,nv), intent(out) :: data
!
      call fatal_error ('input_hdf5_4D', 'You can not use HDF5 without setting an HDF5_IO module.')
      call keep_compiler_quiet(name)
      data(:,:,:,:) = -1.0
!
    endsubroutine input_hdf5_4D
!***********************************************************************
    subroutine output_hdf5_string(name, data)
!
      character (len=*), intent(in) :: name
      character (len=*), intent(in) :: data
!
      call fatal_error ('output_hdf5_string', 'You can not use HDF5 without setting an HDF5_IO module.')
      call keep_compiler_quiet(name,data)
!
    endsubroutine output_hdf5_string
!***********************************************************************
    subroutine output_hdf5_torus_rect(name, data)
!
      use Geometrical_types, only: torus_rect

      character (len=*), intent(in) :: name
      type(torus_rect), intent(in) :: data

      call fatal_error ('output_hdf5_torus_rect', 'You can not use HDF5 without setting an HDF5_IO module.')
      call keep_compiler_quiet(name)
      call keep_compiler_quiet(data%height)
!
    endsubroutine output_hdf5_torus_rect
!***********************************************************************
    subroutine output_hdf5_int_0D(name, data)
!
      character (len=*), intent(in) :: name
      integer, intent(in) :: data
!
      call fatal_error ('output_hdf5_int_0D', 'You can not use HDF5 without setting an HDF5_IO module.')
      call keep_compiler_quiet(name)
      call keep_compiler_quiet(data)
!
    endsubroutine output_hdf5_int_0D
!***********************************************************************
    subroutine output_hdf5_int_1D(name, data, nv, same_size)
!
      character (len=*), intent(in) :: name
      integer, intent(in) :: nv
      integer, dimension(nv), intent(in) :: data
      logical, optional, intent(in) :: same_size
!
      call fatal_error ('output_hdf5_int_1D', 'You can not use HDF5 without setting an HDF5_IO module.')
      call keep_compiler_quiet(name)
      call keep_compiler_quiet(data)
      call keep_compiler_quiet(.true., same_size)
!
    endsubroutine output_hdf5_int_1D
!***********************************************************************
    subroutine output_hdf5_0D(name, data)
!
      character (len=*), intent(in) :: name
      real, intent(in) :: data
!
      call fatal_error ('output_hdf5_0D', 'You can not use HDF5 without setting an HDF5_IO module.')
      call keep_compiler_quiet(name)
      call keep_compiler_quiet(data)
!
    endsubroutine output_hdf5_0D
!***********************************************************************
    subroutine output_hdf5_1D(name, data, nv, same_size)
!
      character (len=*), intent(in) :: name
      integer, intent(in) :: nv
      real, dimension (nv), intent(in) :: data
      logical, optional, intent(in) :: same_size
!
      call fatal_error ('output_hdf5_1D', 'You can not use HDF5 without setting an HDF5_IO module.')
      call keep_compiler_quiet(name)
      call keep_compiler_quiet(data)
      call keep_compiler_quiet(.true., same_size)
!
    endsubroutine output_hdf5_1D
!***********************************************************************
    subroutine output_hdf5_part_2D(name, data, mv, nc, nv)
!
      character (len=*), intent(in) :: name
      integer, intent(in) :: mv, nc
      real, dimension (mv,mparray), intent(in) :: data
      integer, intent(in) :: nv
!
      call fatal_error ('output_hdf5_part_2D', 'You can not use HDF5 without setting an HDF5_IO module.')
      call keep_compiler_quiet(name)
      call keep_compiler_quiet(data)
      call keep_compiler_quiet(mv)
      call keep_compiler_quiet(nc)
      call keep_compiler_quiet(nv)
!
    endsubroutine output_hdf5_part_2D
!***********************************************************************
    subroutine output_hdf5_profile_1D(name, data, ldim, gdim, ip, np1, np2, ng, lhas_data)
!
      character (len=*), intent(in) :: name
      real, dimension (:), intent(in) :: data
      integer, intent(in) :: ldim, gdim, ip, np1, np2, ng
      logical, intent(in) :: lhas_data
!
      call fatal_error ('output_hdf5_profile_1D', 'You can not use HDF5 without setting an HDF5_IO module.')
      call keep_compiler_quiet(name)
      call keep_compiler_quiet(data)
      call keep_compiler_quiet(ldim)
      call keep_compiler_quiet(gdim)
      call keep_compiler_quiet(ip)
      call keep_compiler_quiet(np1)
      call keep_compiler_quiet(np2)
      call keep_compiler_quiet(ng)
      call keep_compiler_quiet(lhas_data)
!
    endsubroutine output_hdf5_profile_1D
!***********************************************************************
    subroutine output_local_hdf5_2D(name, data, dim1, dim2)
!
      character (len=*), intent(in) :: name
      integer, intent(in) :: dim1, dim2
      real, dimension (dim1,dim2), intent(in) :: data
!
      call fatal_error ('output_local_hdf5_2D', 'You can not use HDF5 without setting an HDF5_IO module.')
      call keep_compiler_quiet(name)
      call keep_compiler_quiet(data)
      call keep_compiler_quiet(dim1)
      call keep_compiler_quiet(dim2)
!
    endsubroutine output_local_hdf5_2D
!***********************************************************************
    subroutine output_hdf5_slice_2D(name, data, ldim1, ldim2, gdim1, gdim2, ip1, ip2, lhas_data)
!
      character (len=*), intent(in) :: name
      real, dimension (:,:), pointer :: data
      integer, intent(in) :: ldim1, ldim2, gdim1, gdim2, ip1, ip2
      logical, intent(in) :: lhas_data
!
      call fatal_error ('output_hdf5_slice_2D', 'You can not use HDF5 without setting an HDF5_IO module.')
      call keep_compiler_quiet(name)
      call keep_compiler_quiet(data)
      call keep_compiler_quiet(ldim1)
      call keep_compiler_quiet(ldim2)
      call keep_compiler_quiet(gdim1)
      call keep_compiler_quiet(gdim2)
      call keep_compiler_quiet(ip1)
      call keep_compiler_quiet(ip2)
      call keep_compiler_quiet(lhas_data)
!
    endsubroutine output_hdf5_slice_2D
!***********************************************************************
    subroutine output_local_hdf5_3D(name, data, dim1, dim2, dim3)
!
      character (len=*), intent(in) :: name
      integer, intent(in) :: dim1, dim2, dim3
      real, dimension (dim1,dim2,dim3), intent(in) :: data
!
      call fatal_error ('output_local_hdf5_3D', 'You can not use HDF5 without setting an HDF5_IO module.')
      call keep_compiler_quiet(name)
      call keep_compiler_quiet(data)
      call keep_compiler_quiet(dim1)
      call keep_compiler_quiet(dim2)
      call keep_compiler_quiet(dim3)
!
    endsubroutine output_local_hdf5_3D
!***********************************************************************
    subroutine output_hdf5_3D(name, data)
!
      character (len=*), intent(in) :: name
      real, dimension (mx,my,mz), intent(in) :: data
!
      call fatal_error ('output_hdf5_3D', 'You can not use HDF5 without setting an HDF5_IO module.')
      call keep_compiler_quiet(name)
      call keep_compiler_quiet(data)
!
    endsubroutine output_hdf5_3D
!***********************************************************************
    subroutine output_hdf5_4D(name, data, nv)
!
      character (len=*), intent(in) :: name
      integer, intent(in) :: nv
      real, dimension (mx,my,mz,nv), intent(in) :: data
!
      call fatal_error ('output_hdf5_4D', 'You can not use HDF5 without setting an HDF5_IO module.')
      call keep_compiler_quiet(name)
      call keep_compiler_quiet(data)
!
    endsubroutine output_hdf5_4D
!***********************************************************************
    subroutine input_dim(wrkdir, mx_in, my_in, mz_in, mvar_in, maux_in, mglobal_in, &
                         prec_in, nghost_in, nprocx_in, nprocy_in, nprocz_in, local)
!
!  Read dimensions from dim.dat (local or global).
!
      use General, only: loptest

      character (len=*), intent(in) :: wrkdir
      integer, intent(out) :: mx_in, my_in, mz_in, mvar_in, maux_in, mglobal_in
      integer, intent(out) :: nprocx_in, nprocy_in, nprocz_in, nghost_in
      logical, optional :: local
!
      character (len=fnlen) :: filename
      character :: prec_in
      integer :: iprocz_slowest, ipx_in, ipy_in, ipz_in
!
      if (lroot) then
        ! local or global dimension file
        if (loptest(local)) then
          filename = trim(wrkdir)//'/data/proc0/dim.dat'
        else
          filename = trim(wrkdir)//'/data/dim.dat'
        endif
        open(lun_input,file=filename)
        read(lun_input,'(3i7,3i7)') mx_in, my_in, mz_in, mvar_in, maux_in, mglobal_in
        read(lun_input,'(a)') prec_in
        read(lun_input,'(3i5)') nghost_in
        if (loptest(local)) then
          read(lun_input,'(4i5)') ipx_in, ipy_in, ipz_in
          nprocx_in=ipx_in; nprocy_in=ipy_in; nprocz_in=ipz_in
        else
          read(lun_input,'(4i5)') nprocx_in, nprocy_in, nprocz_in, iprocz_slowest
        endif
        close(lun_input)
      endif

    endsubroutine input_dim
!***********************************************************************
    subroutine output_hdf5_double_0D(name, data)
!
      character (len=*), intent(in) :: name
      double precision, intent(in) :: data
!
      call fatal_error ('output_hdf5_double_0D', 'You can not use HDF5 without setting an HDF5_IO module.')
      call keep_compiler_quiet(name)
!
    endsubroutine output_hdf5_double_0D
!***********************************************************************
    subroutine output_hdf5_double_1D(name, data, nv)
!
      character (len=*), intent(in) :: name
      integer, intent(in) :: nv
      double precision, dimension (nv), intent(in) :: data
!
      call fatal_error ('output_hdf5_double_1D', 'You can not use HDF5 without setting an HDF5_IO module.')
      call keep_compiler_quiet(name)
      call keep_compiler_quiet(nv)
!
    endsubroutine output_hdf5_double_1D
!***********************************************************************
    subroutine output_dim(file, mx_out, my_out, mz_out, mxgrid_out, mygrid_out, mzgrid_out, mvar_out, maux_out, mglobal)
!
!  Write dimension to file.
!
!  02-Nov-2018/PABourdin: moved IO parts from wdim to low-level IO modules
!
      character (len=*), intent(in) :: file
      integer, intent(in) :: mx_out, my_out, mz_out, mxgrid_out, mygrid_out, mzgrid_out, mvar_out, maux_out, mglobal
!
      character (len=fnlen) :: filename
      integer :: iprocz_slowest
!
      if (lroot) then
        ! global dimension file
        filename = trim(datadir)//'/'//trim(file)
        open(lun_output,file=filename)
        write(lun_output,'(3i7,3i7)') mxgrid_out, mygrid_out, mzgrid_out, mvar_out, maux_out, mglobal
        write(lun_output,'(a)') numeric_precision()
        write(lun_output,'(3i5)') nghost, nghost, nghost ! why do we write three always identical numbers?
        iprocz_slowest = 0
        if (lprocz_slowest) iprocz_slowest=1
        write(lun_output,'(4i5)') nprocx, nprocy, nprocz, iprocz_slowest
        close(lun_output)
      endif
!
      if (.not. lmonolithic_io) then
        ! write processor dimension files to "proc#/" for distributed IO cases
        filename = trim(directory)//'/'//trim(file)
        open(lun_output,file=filename)
        write(lun_output,'(3i7,3i7)') mx_out, my_out, mz_out, mvar_out, maux_out, mglobal
        write(lun_output,'(a)') numeric_precision()
        write(lun_output,'(3i5)') nghost, nghost, nghost ! why do we write three always identical numbers?
        write(lun_output,'(3i5)') ipx, ipy, ipz
        close(lun_output)
      endif
!
    endsubroutine output_dim
!***********************************************************************
    subroutine wdim_default_grid(file)
!
!  Write dimension to file.
!
!  02-Nov-2018/PABourdin: redesigned
!
      character (len=*), intent(in) :: file
!
      call output_dim (file, mx, my, mz, mxgrid, mygrid, mzgrid, mvar, maux, mglobal)
!
    endsubroutine wdim_default_grid
!***********************************************************************
    subroutine wdim_default(file, mx_out, my_out, mz_out, mxgrid_out, mygrid_out, mzgrid_out)
!
!  Write dimension to file.
!
!  02-Nov-2018/PABourdin: redesigned
!
      character (len=*), intent(in) :: file
      integer, intent(in) :: mx_out, my_out, mz_out, mxgrid_out, mygrid_out, mzgrid_out
!
      call output_dim (file, mx_out, my_out, mz_out, mxgrid_out, mygrid_out, mzgrid_out, mvar, maux, mglobal)
!
    endsubroutine wdim_default
!***********************************************************************
    subroutine wdim(file, mx_out, my_out, mz_out, mxgrid_out, mygrid_out, mzgrid_out, mvar_out, maux_out)
!
!  Write dimension to file.
!
!   8-sep-01/axel: adapted to take my_out,mz_out
!   4-oct-16/MR: added optional parameters mvar_out,maux_out
!  02-Nov-2018/PABourdin: redesigned, moved to IO modules
!
      character (len=*), intent(in) :: file
      integer, intent(in) :: mx_out, my_out, mz_out, mxgrid_out, mygrid_out, mzgrid_out, mvar_out, maux_out
!
      call output_dim (file, mx_out, my_out, mz_out, mxgrid_out, mygrid_out, mzgrid_out, mvar_out, maux_out, mglobal)
!
    endsubroutine wdim
!***********************************************************************
    subroutine output_timeseries(data, data_im)
!   
!  Append diagnostic data to a binary file.
!   
!  01-Apr-2019/PABourdin: coded
! 
      real, dimension(2*nname), intent(in) :: data, data_im
! 
      ! dummy routine   
!                       
    endsubroutine output_timeseries
!***********************************************************************
    subroutine input_slice_position(directory,ix_bc_,iy_bc_,iy2_bc_,iz_bc_,iz2_bc_,iz3_bc_,iz4_bc_)
!
!  'data/procN/slice_position.dat' is distributed, but may not be synchronized
!  on I/O error (-> dist=0) as this would make it disfunctional; correct a posteriori if necessary.
!
!  24-May-2019/MR: cloned from output_slice_position
!
      character(LEN=*) :: directory
      integer, intent(out) :: ix_bc_,iy_bc_,iy2_bc_,iz_bc_,iz2_bc_,iz3_bc_,iz4_bc_
      logical :: lexist_slice_xy, lexist_slice_xy2, lexist_slice_xy3, lexist_slice_xy4, &
                 lexist_slice_xz, lexist_slice_xz2, lexist_slice_yz

      open (lun_input, file=trim(directory)//'/data/slice_position.dat', STATUS='unknown')
      read (lun_input, '(l5,i5)') lexist_slice_xy, iz_bc_
      read (lun_input, '(l5,i5)') lexist_slice_xy2, iz2_bc_
      read (lun_input, '(l5,i5)') lexist_slice_xy3, iz3_bc_
      read (lun_input, '(l5,i5)') lexist_slice_xy4, iz4_bc_
      read (lun_input, '(l5,i5)') lexist_slice_xz, iy_bc_
      read (lun_input, '(l5,i5)') lexist_slice_xz2, iy2_bc_
      read (lun_input, '(l5,i5)') lexist_slice_yz, ix_bc_
      close(lun_input)
!
    endsubroutine input_slice_position
!***********************************************************************
    subroutine input_slice_real_arr(file, time, pos, data)
!
!  read a slice file
!
!  24-may-19/MR: coded
!
      use File_io, only: file_exists

      character (len=*),   intent(in) :: file
      real,                intent(out):: time
      real,                intent(out):: pos
      real, dimension(:,:,:),intent(out):: data
!
      integer :: nt, ios
!
      if (.not.file_exists(file)) &
        call fatal_error('input_slice', 'no slices file '//trim(file))

      open(lun_input, file=file, form='unformatted')
      nt=0; ios=0
      do while(ios==0)
        read(lun_input,iostat=ios) data(:,:,nt+1), time, pos
        if (ios/=0) exit
        nt=nt+1
      enddo
      close(lun_input)
!
    endsubroutine input_slice_real_arr
!***********************************************************************
    subroutine input_slice_scat_arr(file, pos, data, ivar, nt)
!
!  read a slice file
!
!  24-may-19/MR: coded
!
      use General, only: scattered_array, store_scattered_array
      use File_io, only: file_exists

      character (len=*),   intent(in) :: file
      real,                intent(out):: pos
      type(scattered_array), pointer  :: data   !intent(inout)
      integer,             intent(in) :: ivar
      integer,             intent(out):: nt
!
      integer :: ios
      real :: time
      real,dimension(data%dim1,data%dim2) :: slc
!
      if (.not.file_exists(file)) &
        call fatal_error('input_slice', 'no slices file '//trim(file))

      open(lun_input, file=file, form='unformatted')
      nt=0; ios=0
      do while(ios==0)
        read(lun_input,iostat=ios) slc, time, pos
        if (ios/=0) exit
        nt=nt+1
!if (ivar==1.and.iproc==120) print*, 'nt=', nt
        call store_scattered_array(ivar,nt,slc,data,time)
      enddo
      close(lun_input)
!
    endsubroutine input_slice_scat_arr
!***********************************************************************
    subroutine hdf5_output_slice_position
!
!  'data/procN/slice_position.dat' is distributed, but may not be synchronized
!  on I/O error (-> dist=0) as this would make it disfunctional; correct a posteriori if necessary.
!
!  27-Oct-2018/PABourdin: cleaned up
!
      use Slices_methods, only: write_rslice_position

      open (lun_output, file=trim(directory)//'/slice_position.dat', STATUS='unknown')
      write (lun_output, '(l5,i5," XY")' ) lwrite_slice_xy, iz_loc
      write (lun_output, '(l5,i5," XY2")') lwrite_slice_xy2, iz2_loc
      write (lun_output, '(l5,i5," XY3")') lwrite_slice_xy3, iz3_loc
      write (lun_output, '(l5,i5," XY4")') lwrite_slice_xy4, iz4_loc
      write (lun_output, '(l5,i5," XZ")' ) lwrite_slice_xz, iy_loc
      write (lun_output, '(l5,i5," XZ2")') lwrite_slice_xz2, iy2_loc
      write (lun_output, '(l5,i5," YZ")' ) lwrite_slice_yz, ix_loc
      call write_rslice_position(lun_output)

      close (lun_output)
!
    endsubroutine hdf5_output_slice_position
!***********************************************************************
    subroutine hdf5_output_slice(lwrite, time, label, suffix, pos, grid_pos, data)
!
!  append to a slice file
!
!  12-nov-02/axel: coded
!  26-jun-06/anders: moved to slices
!  22-sep-07/axel: changed Xy to xy2, to be compatible with Mac
!  28-Oct-2018/PABourdin: cleaned up, moved to IO module
!
      logical, intent(in) :: lwrite
      real, intent(in) :: time
      character (len=*), intent(in) :: label, suffix
      real, intent(in) :: pos
      integer, intent(in) :: grid_pos
      real, dimension (:,:), pointer :: data
!
      if (.not. lwrite .or. .not. associated(data)) return
!
!  files data/procN/slice*.* are distributed and will be synchronized a-posteriori on I/O error
!
      open (lun_output, file=trim(directory)//'/slice_'//trim(label)//'.'//trim(suffix), &
            form='unformatted', position='append')
      write (lun_output) data, time, pos
      close (lun_output)
!
    endsubroutine hdf5_output_slice
!***********************************************************************
    subroutine index_append(varname,ivar,vector,array)
!
! 14-oct-18/PAB: coded
! 09-Jul-2020/PAB: reworked
!
      character (len=*), intent(in) :: varname
      integer, intent(in) :: ivar
      integer, intent(in) :: vector
      integer, intent(in) :: array
!
      character (len=len(varname)) :: component
      integer :: pos, l
!
      if (lroot) then
        open (lun_output, file=trim (datadir)//'/'//trim (index_pro), POSITION='append')
        if ((vector > 0) .and. (array > 0)) then
          ! expand array first: iuud => [iuud1,iuud2,iuud3,...]
          ! expand vector then: iuud# => [iuud#x,iuud#y,iuud#z] => ivar+(#-1)*vector+[0,1,2]
          do pos = 1, array
            write (lun_output,*) trim(varname)//trim(itoa(pos))//'x='//trim(itoa(ivar+(pos-1)*vector))
            write (lun_output,*) trim(varname)//trim(itoa(pos))//'y='//trim(itoa(ivar+(pos-1)*vector+1))
            write (lun_output,*) trim(varname)//trim(itoa(pos))//'z='//trim(itoa(ivar+(pos-1)*vector+2))
          enddo
        elseif ((vector > 0) .and. (array < 0)) then
          write (lun_output,*) trim(varname)//'=indgen('//trim(itoa(-array))//')*'//trim(itoa(vector))//'+'//trim(itoa(ivar))
        elseif (array /= 0) then    ! i.e. vector<0
          ! backwards compatibility: ind => indgen(array)+ivar
          write (lun_output,*) trim (varname)//'=indgen('//trim(itoa(abs(array)))//')+'//trim (itoa (ivar))
          ! expand array: ind => [ind1,ind2,...]
          do pos = 1, array
            write (lun_output,*) trim(varname)//trim(itoa(pos))//'='//trim(itoa(ivar+(pos-1)))
          enddo
        elseif (vector > 0) then    ! i.e. array=0
          ! expand vectors: iuu => [iux,iuy,iuz], iaa => [iax,iay,iaz], etc.
          if (vector == 3) then
            ! backwards compatibility: write original vector
            write (lun_output,*) trim(varname)//'='//trim (itoa(ivar))
            ! apply shortcuts
            component = trim (varname)
            l = len (trim (component))
            if (l == 3) then
              ! double endings: iuu, iaa, etc.
              if (component(2:2) == component(3:3)) l = 2
            endif
            write (lun_output,*) trim(component(1:l))//'x='//trim (itoa(ivar))
            write (lun_output,*) trim(component(1:l))//'y='//trim (itoa(ivar+1))
            write (lun_output,*) trim(component(1:l))//'z='//trim (itoa(ivar+2))
          elseif (vector == 6) then
            ! expand symmetric 3x3 tensor (6 different components)
            write (lun_output,*) trim(varname)//'_xx='//trim (itoa(ivar))
            write (lun_output,*) trim(varname)//'_xy='//trim (itoa(ivar+1))
            write (lun_output,*) trim(varname)//'_xz='//trim (itoa(ivar+2))
            write (lun_output,*) trim(varname)//'_yy='//trim (itoa(ivar+3))
            write (lun_output,*) trim(varname)//'_yz='//trim (itoa(ivar+4))
            write (lun_output,*) trim(varname)//'_zz='//trim (itoa(ivar+5))
          elseif (vector==9) then
            ! expand asymmetric 3x3 tensor (9 different components)
            write (lun_output,*) trim(varname)//'_xx='//trim (itoa(ivar))
            write (lun_output,*) trim(varname)//'_xy='//trim (itoa(ivar+1))
            write (lun_output,*) trim(varname)//'_xz='//trim (itoa(ivar+2))
            write (lun_output,*) trim(varname)//'_yx='//trim (itoa(ivar+3))
            write (lun_output,*) trim(varname)//'_yy='//trim (itoa(ivar+4))
            write (lun_output,*) trim(varname)//'_yz='//trim (itoa(ivar+5))
            write (lun_output,*) trim(varname)//'_zx='//trim (itoa(ivar+6))
            write (lun_output,*) trim(varname)//'_zy='//trim (itoa(ivar+7))
            write (lun_output,*) trim(varname)//'_zz='//trim (itoa(ivar+8))
          else
            call fatal_error('index_append','unknown tensor type')
          endif
        else
          ! scalar: ilnrho => ivar
          write (lun_output,*) trim(varname)//'='//trim(itoa(ivar))
        endif
        close (lun_output)
      endif
!
    endsubroutine index_append
!***********************************************************************
    subroutine particle_index_append(label,ilabel)
!
! 22-Oct-2018/PABourdin: coded
!
      character (len=*), intent(in) :: label
      integer, intent(in) :: ilabel
!
      if (lroot) then
        open(lun_output,file=trim(datadir)//'/'//trim(particle_index_pro), POSITION='append')
        write(lun_output,*) trim(label)//'='//trim(itoa(ilabel))
        close(lun_output)
      endif
!
    endsubroutine particle_index_append
!***********************************************************************
    subroutine pointmass_index_append(label,ilabel)
!
! 13-Apr-2019/PABourdin: copied from 'particle_index_append'
!
      character (len=*), intent(in) :: label
      integer, intent(in) :: ilabel
!
      integer, parameter :: lun_output = 92
!
      if (lroot) then
        open(lun_output,file=trim(datadir)//'/'//trim(pointmass_index_pro), POSITION='append')
        write(lun_output,*) trim(label)//'='//trim(itoa(ilabel))
        close(lun_output)
      endif
!
    endsubroutine pointmass_index_append
!***********************************************************************
    function index_get(ivar,particle,pointmass,quiet)
!
! 13-Nov-2018/PABourdin: coded
!
      character (len=labellen) :: index_get
      integer, intent(in) :: ivar
      logical, optional, intent(in) :: particle, pointmass, quiet
!
      call fatal_error ('index_get', 'You can not use HDF5 without setting an HDF5_IO module.')
      call keep_compiler_quiet(ivar)
      call keep_compiler_quiet(particle)
      call keep_compiler_quiet(pointmass)
      call keep_compiler_quiet(quiet)
!
      index_get=' '
!
    endfunction index_get
!***********************************************************************
    subroutine index_reset
!
! 14-Oct-2018/PABourdin: coded
!
      if (lroot) then
! Axel, please uncomment for debugging on Beskow:
! write (123,*) '============> replace "index.pro" with empty file'
        open(lun_output,file=trim(datadir)//'/'//trim(index_pro),status='replace')
        close(lun_output)
        open(lun_output,file=trim(datadir)//'/'//trim(particle_index_pro),status='replace')
        close(lun_output)
        open(lun_output,file=trim(datadir)//'/'//trim(pointmass_index_pro),status='replace')
        close(lun_output)
      endif
!
    endsubroutine index_reset
!***********************************************************************
    subroutine output_profile(fname, coord, a, type, lsave_name, lhas_ghost)
!
!  Writes a profile to a file.
!
!  10-jul-05/axel: coded
!  07-Nov-2018/PABourdin: moved to IO modules
!
      use General, only: loptest
!
      real, dimension(:) :: coord, a
      character (len=*) :: fname
      character :: type
      logical, optional :: lsave_name, lhas_ghost
!
      integer :: pos, np
!
!  If within a loop, do this only for the first step (indicated by lwrite_prof).
!
      if (lwrite_prof) then
        ! Write profile.
        open(lun_output,file=trim(directory)//'/'//type//'prof_'//trim(fname)//'.dat',position='append')
        np = size(coord)
        do pos=1, np
          write(lun_output,*) coord(pos), a(pos)
        enddo
        close(lun_output)
!
        ! Add file name to list of profiles if lsave_name is true.
        if (loptest(lsave_name)) then
          open(lun_output,file=trim(directory)//'/'//type//'prof_list.dat',position='append')
          write(lun_output,*) fname
          close(lun_output)
        endif
      endif
!
    endsubroutine output_profile
!***********************************************************************
    subroutine input_profile(fname, type, a, np, lhas_ghost)
!
!  Read a profile from a file (taken from stratification, for example).
!
!  15-may-15/axel: coded
!  05-Nov-2018/PABourdin: generalized to any direction and moved to IO module
!
      character (len=*), intent(in) :: fname
      character, intent(in) :: type
      integer, intent(in) :: np
      real, dimension(np), intent(out) :: a
      logical, optional :: lhas_ghost
!
      integer :: pos
      real, dimension(np) :: coord
!
      ! Read profile.
      open(lun_input,file=trim(directory)//type//'/prof_'//trim(fname)//'.dat')
      do pos=1, np
        read(lun_input,*) coord(pos), a(pos)
      enddo
      close(lun_input)
!
!  Should we check that coord == z for type == 'z'?
!
    endsubroutine input_profile
!***********************************************************************
    subroutine output_average_1D(path, label, nc, name, data, time, lbinary, lwrite, header)
!
!   Output 1D average to a file.
!
!   16-Nov-2018/PABourdin: coded
!
      character (len=*), intent(in) :: path, label
      integer, intent(in) :: nc
      character (len=fmtlen), dimension(nc), intent(in) :: name
      real, dimension(:,:), intent(in) :: data
      real, intent(in) :: time
      logical, intent(in) :: lbinary, lwrite
      real, dimension(:), optional, intent(in) :: header
!
      character (len=fnlen) :: filename
!
      if (.not. lwrite .or. (nc <= 0)) return
!
      filename = trim(path) // '/' // trim(label) // 'averages.dat'
      if (trim (label) == 'phi_z') filename = trim(path) // '/phizaverages.dat'
      if (lbinary) then
        open(lun_output, file=filename, form='unformatted', position='append')
        if (present (header)) write(lun_output) header
        write(lun_output) time
        write(lun_output) data(:,1:nc)
        close(lun_output)
      else
        open(lun_output, file=filename, position='append')
        if (present (header)) write(lun_output,'(1p,8'//trim(fmt_avgs)//')') header
        write(lun_output,'(1pe12.5)') time
        write(lun_output,'(1p,8'//trim(fmt_avgs)//')') data(:,1:nc)
        close(lun_output)
      endif
!
    endsubroutine output_average_1D
!***********************************************************************
    subroutine output_average_1D_chunked(path, label, nc, name, data, full, time, lbinary, lwrite, header)
!
!   Output 1D chunked average to a file.
!
!   16-Nov-2018/PABourdin: coded
!
      character (len=*), intent(in) :: path, label
      integer, intent(in) :: nc
      character (len=fmtlen), dimension(nc), intent(in) :: name
      real, dimension(:,:,:), intent(in) :: data
      integer, intent(in) :: full
      real, intent(in) :: time
      logical, intent(in) :: lbinary, lwrite
      real, dimension(:), optional, intent(in) :: header
!
      if (present (header)) then
        call output_average_2D(path, label, nc, name, data, time, lbinary, lwrite, header)
      else
        call output_average_2D(path, label, nc, name, data, time, lbinary, lwrite)
      endif
!
    endsubroutine output_average_1D_chunked
!***********************************************************************
    subroutine output_average_2D(path, label, nc, name, data, time, lbinary, lwrite, header)
!
!   Output average to a file.
!
!   16-Nov-2018/PABourdin: coded
!
      character (len=*), intent(in) :: path, label
      integer, intent(in) :: nc
      character (len=fmtlen), dimension(nc), intent(in) :: name
      real, dimension(:,:,:), intent(in) :: data
      real, intent(in) :: time
      logical, intent(in) :: lbinary, lwrite
      real, dimension(:), optional, intent(in) :: header
!
      character (len=fnlen) :: filename
      integer :: ia
!
      if (.not. lwrite .or. (nc <= 0)) return
!
      filename = trim(path) // '/' // trim(label) // 'averages.dat'
      if (lbinary) then
        open(lun_output, file=filename, form='unformatted', position='append')
        if (present (header)) write(lun_output) header
        write(lun_output) time
        if (label(1:1) == 'z') then
          write(lun_output) ( data(ia,:,:), ia=1, nc )
        else
          write(lun_output) data(:,:,1:nc)
        endif
        close(lun_output)
      else
        open(lun_output, file=filename, position='append')
        if (present (header)) write(lun_output,'(1p,8'//trim(fmt_avgs)//')') header
        write(lun_output,'(1pe12.5)') time
        if (label(1:1) == 'z') then
          write(lun_output,'(1p,8'//trim(fmt_avgs)//')') ( data(ia,:,:), ia=1, nc )
        else
          write(lun_output,'(1p,8'//trim(fmt_avgs)//')') data(:,:,1:nc)
        endif
        close(lun_output)
      endif
!
    endsubroutine output_average_2D
!***********************************************************************
    subroutine output_average_phi(path, number, nr, nc, name, data, time, r, dr)
!       
!   Output phi average to a file with these records:
!   1) nr_phiavg, nz_phiavg, nvars, nprocz
!   2) t, r_phiavg, z_phiavg, dr, dz
!   3) data
!   4) len(labels),labels
!
!   27-Nov-2014/PABourdin: cleaned up code from write_phiaverages
!   25-Nov-2018/PABourdin: coded
!         
      use General, only: safe_character_append
!     
      character (len=*), intent(in) :: path, number
      integer, intent(in) :: nr, nc 
      character (len=fmtlen), dimension(nc), intent(in) :: name
      real, dimension(:,:,:,:), intent(in) :: data
      real, intent(in) :: time
      real, dimension(nr), intent(in) :: r
      real, intent(in) :: dr
!
      character (len=fnlen) :: filename
      character (len=1024) :: labels
      integer :: pos
!     
      if (.not. lroot .or. (nc <= 0)) return
!
      filename = 'PHIAVG' // trim(number)
      open(lun_output, file=trim(path)//'/averages/'//trim(filename), form='unformatted', position='append')
      write(lun_output) nr, nzgrid, nc, nprocz
      write(lun_output) time, r, z(n1)+(/(pos*dz, pos=0, nzgrid-1)/), dr, dz
      ! note: due to passing data as implicit-size array,
      ! the indices (0:nz) are shifted to (1:nz+1),
      ! so that we have to write only the portion (2:nz+1).
      ! data has to be repacked to avoid writing an array temporary
      ! ngrs: writing an array temporary outputs corrupted data on copson
      write(lun_output) pack(data(:,2:nz+1,:,1:nc), .true.)
      labels = trim(name(1))
      if (nc >= 2) then
        do pos = 2, nc
          call safe_character_append (labels, ",", trim(name(pos)))
        enddo
      endif
      write(lun_output) len(labels), labels
      close(lun_output)
!
      ! append filename to list
      open (lun_output, FILE=trim(datadir)//'/averages/phiavg.files', POSITION='append')
      write (lun_output,'(A)') trim(filename)
      close (lun_output)
!
    endsubroutine output_average_phi
!***********************************************************************
    subroutine trim_average(path, plane, ngrid, nname)
!
!  Trim a 1D-average file for times past the current time.
!
!  25-apr-16/ccyang: coded
!  23-Nov-2018/PABourdin: moved to IO module
!
      use File_io, only: file_exists

      character (len=*), intent(in) :: path, plane
      integer, intent(in) :: ngrid, nname
!
      real, dimension(:), allocatable :: tmp
      character(len=fnlen) :: filename
      integer :: pos, num_rec, ioerr, alloc_err
      integer :: lun_input = 84
      real :: time
!
      if (.not. lroot) return
      if ((ngrid <= 0) .or. (nname <= 0)) return
      filename = trim(path) // '/' // trim(plane) // 'averages.dat'
      if (.not. file_exists (filename)) return
!
      allocate (tmp(ngrid * nname), stat = alloc_err)
      if (alloc_err > 0) call fatal_error('trim_average', 'Could not allocate memory for averages')
!
      if (lwrite_avg1d_binary) then
        open(lun_input, file=filename, form='unformatted', action='readwrite')
      else
        open(lun_input, file=filename, action='readwrite')
      endif
!
      ! count number of records
      num_rec = 0
      ioerr = 0
      time = 0.0
      do while ((t > time) .and. (ioerr == 0))
        num_rec = num_rec + 1
        if (lwrite_avg1d_binary) then
          read(lun_input, iostat=ioerr) time
          read(lun_input, iostat=ioerr) tmp
        else
          read(lun_input, *, iostat=ioerr) time
          read(lun_input, *, iostat=ioerr) tmp
        endif
      enddo
      if (time == t) num_rec = num_rec - 1
!
      if ((time >= t) .and. (ioerr == 0)) then
        ! trim excess data at the end of the average file
        rewind(lun_input)
        if (num_rec > 0) then
          do pos = 1, num_rec
            if (lwrite_avg1d_binary) then
              read(lun_input, iostat=ioerr) time
              read(lun_input, iostat=ioerr) tmp
            else
              read(lun_input, *, iostat=ioerr) time
              read(lun_input, *, iostat=ioerr) tmp
            endif
          enddo
        endif
        endfile (lun_input)
        if (ip <= 10) print *, 'trim_average: trimmed '//trim(plane)//'-averages for t >= ', t
      endif
!
      close (lun_input)
      deallocate (tmp)
!
    endsubroutine trim_average
!***********************************************************************
endmodule HDF5_IO