WriteXDMF Subroutine

  subroutine WriteXDMF(time,nout)
    use, intrinsic :: iso_fortran_env, only: int64,real64
    implicit none
    real(8),intent(in)::time
    integer, intent(in) :: nout
    integer:: itot, jtot, ktot
    integer:: n

    character(*), parameter :: dirname = datadir
    character(256) :: xmfname, fgridx, fgridy, fgridz, fdata
    integer :: u, ncell
    integer(int64) :: bytes_per_real, bytes_per_field
    integer(int64) :: off_base


    itot = ntotal(1)
    jtot = ntotal(2)
    ktot = ntotal(3)
    
    ! ---- file names (match your ReadData) ----
    write(xmfname,'(a,i5.5,a)') trim(filefield), nout, ".xmf"
    xmfname = trim(dirname)//trim(xmfname)

    fgridx = filegrid1D
    fgridy = filegrid2D
    fgridz = filegrid3D
    write(fdata,'(a,i5.5,a)') trim(filefield), nout,trim(binaryextension)

    ! ---- sizes & offsets ----
    ! stream/unformatted wrote raw reals; assume real64 (8 bytes) because iso_fortran_env real64
    bytes_per_real = int(storage_size(0.0_real64)/8, int64)
    
    ncell = itot*jtot*ktot
    bytes_per_field = int(ncell, int64) * bytes_per_real

    ! ---- write XDMF (XML) ----
    open(newunit=u, file=xmfname, status="replace", action="write", form="formatted")
    
    write(u,'(a)') '<?xml version="1.0" ?>'
    write(u,'(a)') '<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>'
    write(u,'(a)') '<Xdmf Version="2.0">'
    write(u,'(a)') '  <Domain>'
    write(u,'(a)') '    <Grid Name="Grid" GridType="Uniform">'
    write(u,'(a,es24.16,a)') '      <Time Value="', time, '"/>'

    ! Topology: Rectilinear mesh uses node dimensions = (in+1, jn+1, kn+1)
    ! Use Order="Fortran" consistently so Dimensions list matches (i,j,k).
    write(u,'(a,i0,1x,i0,1x,i0,a)') '      <Topology TopologyType="3DRectMesh" Dimensions="', &
         ktot+1, jtot+1, itot+1, '"/>'

    write(u,'(a)') '      <Geometry GeometryType="VXVYVZ">'
    call write_axis(u, fgridx, itot+1, int(itot+1,int64)*bytes_per_real, bytes_per_real)  ! x1a(:)
    call write_axis(u, fgridy, jtot+1, int(jtot+1,int64)*bytes_per_real, bytes_per_real)  ! x2a(:)
    call write_axis(u, fgridz, ktot+1, int(ktot+1,int64)*bytes_per_real, bytes_per_real)  ! x3a(:)

    write(u,'(a)') '      </Geometry>'

    off_base = 0_int64
    ! ---- Cell-centered attributes (itot,jtot,ktot) ----
    call write_attr(u, "d" , fdata, itot, jtot, ktot, off_base, bytes_per_real); off_base=off_base + 1_int64 * bytes_per_field
    call write_attr(u, "v1", fdata, itot, jtot, ktot, off_base, bytes_per_real); off_base=off_base + 1_int64 * bytes_per_field
    call write_attr(u, "v2", fdata, itot, jtot, ktot, off_base, bytes_per_real); off_base=off_base + 1_int64 * bytes_per_field
    call write_attr(u, "v3", fdata, itot, jtot, ktot, off_base, bytes_per_real); off_base=off_base + 1_int64 * bytes_per_field
    call write_attr(u, "b1", fdata, itot, jtot, ktot, off_base, bytes_per_real); off_base=off_base + 1_int64 * bytes_per_field
    call write_attr(u, "b2", fdata, itot, jtot, ktot, off_base, bytes_per_real); off_base=off_base + 1_int64 * bytes_per_field
    call write_attr(u, "b3", fdata, itot, jtot, ktot, off_base, bytes_per_real); off_base=off_base + 1_int64 * bytes_per_field
    call write_attr(u, "bp", fdata, itot, jtot, ktot, off_base, bytes_per_real); off_base=off_base + 1_int64 * bytes_per_field
    call write_attr(u, "p" , fdata, itot, jtot, ktot, off_base, bytes_per_real); off_base=off_base + 1_int64 * bytes_per_field
    call write_attr(u, "gp", fdata, itot, jtot, ktot, off_base, bytes_per_real); off_base=off_base + 1_int64 * bytes_per_field
    do n=1,ncomp
       call write_attr(u, "X"//trim(adjustl(to_str(n))) , fdata, itot, jtot, ktot, off_base , bytes_per_real); off_base=off_base + 1_int64 * bytes_per_field
    enddo
    write(u,'(a)') '    </Grid>'
    write(u,'(a)') '  </Domain>'
    write(u,'(a)') '</Xdmf>'
    
    close(u)

  contains

    subroutine write_axis(u, fname, n, seek_bytes, bpr)
      use, intrinsic :: iso_fortran_env, only: int64
      implicit none
      integer, intent(in) :: u, n
      character(*), intent(in) :: fname
      integer(int64), intent(in) :: seek_bytes, bpr

      ! XDMF Binary DataItem:
      ! - Format="Binary"
      ! - Seek="bytes" to skip within file
      ! - Precision="8" etc.
      write(u,'(a)') '        <DataItem Dimensions="'//trim(itoa(n))// &
           '" NumberType="Float" Precision="'//trim(itoa(int(bpr)))// &
           '" Format="Binary" Endian="Little" Seek="'//trim(i64toa(seek_bytes))// &
           '"  >'//trim(fname)//'</DataItem>'
    end subroutine write_axis

    subroutine write_attr(u, name, fname, nx, ny, nz, seek_bytes, bpr)
      use, intrinsic :: iso_fortran_env, only: int64
      implicit none
      integer, intent(in) :: u, nx, ny, nz
      character(*), intent(in) :: name, fname
      integer(int64), intent(in) :: seek_bytes, bpr

      write(u,'(a)') '      <Attribute Name="'//trim(name)//'" AttributeType="Scalar" Center="Cell">'
      write(u,'(a)') '        <DataItem Dimensions="'//trim(itoa(nz))//' '//trim(itoa(ny))//' '//trim(itoa(nx))// &
           '" NumberType="Float" Precision="'//trim(itoa(int(bpr)))// &
           '" Format="Binary" Endian="Little" Seek="'//trim(i64toa(seek_bytes))// &
           '" >'//trim(fname)//'</DataItem>'
      write(u,'(a)') '      </Attribute>'
    end subroutine write_attr
    
    function itoa(i) result(s)
      implicit none
      integer, intent(in) :: i
      character(32) :: s
      write(s,'(i0)') i
    end function itoa

    function i64toa(i) result(s)
      use, intrinsic :: iso_fortran_env, only: int64
      implicit none
      integer(int64), intent(in) :: i
      character(64) :: s
      write(s,'(i0)') i
    end function i64toa

  end subroutine WriteXDMF