fileio_vtk.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: fileio_vtk.f90                                                    #
!#                                                                           #
!# Copyright (C) 2010-2017                                                   #
!# Björn Sperling   <sperling@astrophysik.uni-kiel.de>                       #
!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de>                       #
!# Jannes Klee      <jklee@astrophysik.uni-kiel.de>                          #
!#                                                                           #
!# This program is free software; you can redistribute it and/or modify      #
!# it under the terms of the GNU General Public License as published by      #
!# the Free Software Foundation; either version 2 of the License, or (at     #
!# your option) any later version.                                           #
!#                                                                           #
!# This program is distributed in the hope that it will be useful, but       #
!# WITHOUT ANY WARRANTY; without even the implied warranty of                #
!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or        #
!# NON INFRINGEMENT.  See the GNU General Public License for more            #
!# details.                                                                  #
!#                                                                           #
!# You should have received a copy of the GNU General Public License         #
!# along with this program; if not, write to the Free Software               #
!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.                 #
!#                                                                           #
!#############################################################################
!----------------------------------------------------------------------------!
!----------------------------------------------------------------------------!
MODULE fileio_vtk_mod
  USE fileio_base_mod
  USE geometry_base_mod
  USE mesh_base_mod
  USE physics_base_mod
  USE fluxes_base_mod
  USE timedisc_base_mod
  USE sources_base_mod
  USE common_dict
#ifdef PARALLEL
#ifdef HAVE_MPI_MOD
  USE mpi
#endif
#endif
  IMPLICIT NONE
#ifdef PARALLEL
#ifdef HAVE_MPIF_H
  include 'mpif.h'
#endif
#endif
  !--------------------------------------------------------------------------!
  PRIVATE
  !--------------------------------------------------------------------------!
   INTEGER, PARAMETER      :: maxcomp  = 9
   INTEGER, PARAMETER      :: maxcols  = 40
   INTEGER, PARAMETER      :: maxkey   = 64
   CHARACTER, PARAMETER    :: lf = achar(10) 
!   CHARACTER(LEN=16),DIMENSION(6),PARAMETER  :: fluxkey = (/'bflux_WEST  ', &
!                                                            'bflux_EAST  ', &
!                                                            'bflux_SOUTH ', &
!                                                            'bflux_NORTH ', &
!                                                            'bflux_BOTTOM', &
!                                                            'bflux_TOP   ' /)

  !--------------------------------------------------------------------------!
  TYPE, EXTENDS(fileio_base) :: fileio_vtk
    PRIVATE
    INTEGER   :: ip1,jp1,kp1
    CONTAINS
    ! methods
    PROCEDURE :: initfileio_vtk
    PROCEDURE :: openfile
    PROCEDURE :: closefile
    PROCEDURE :: writeheader
    PROCEDURE :: readheader
    PROCEDURE :: writetimestamp
    PROCEDURE :: writeparaviewfile
    PROCEDURE :: readtimestamp
    PROCEDURE :: writedataset
    PROCEDURE :: getprecision
    PROCEDURE :: getoutputlist
    PROCEDURE :: getendianness
    PROCEDURE :: finalize
  END TYPE
  !--------------------------------------------------------------------------!
  PUBLIC :: &
    fileio_vtk

  !--------------------------------------------------------------------------!
CONTAINS
  SUBROUTINE initfileio_vtk(this,Mesh,Physics,Timedisc,Sources,config,IO)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fileio_vtk),   INTENT(INOUT)       :: this
    CLASS(mesh_base),    INTENT(IN)          :: Mesh
    CLASS(physics_base), INTENT(IN)          :: Physics
    CLASS(timedisc_base),INTENT(IN)          :: Timedisc
    CLASS(sources_base), INTENT(IN), POINTER :: Sources
    TYPE(Dict_TYP),      INTENT(IN), POINTER :: config
    TYPE(Dict_TYP),      INTENT(IN), POINTER :: IO
    !------------------------------------------------------------------------!
    TYPE(Dict_TYP), POINTER             :: node
    INTEGER                             :: k,i,j,err
#ifdef PARALLEL
    INTEGER, DIMENSION(:), POINTER      :: sendbuf,recvbuf
    INTEGER                             :: n
#endif
    REAL, DIMENSION(:,:,:,:,:), POINTER :: corners
    !------------------------------------------------------------------------!
#ifdef PARALLEL
    this%multfiles = .true.
#endif
    ! start init form base class in beginning
    CALL this%InitFileIO(mesh,physics,timedisc,sources,config,io,"VTK","vts")

    ! some sanity checks
    IF (this%cycles .NE. this%count+1) CALL this%Error('InitFileIO','VTK need filecycles = count+1')

    ! determine dimensions for the coordinate array
    IF (abs(mesh%xmax-mesh%xmin).LT.2*epsilon(mesh%xmin) &
       .OR.(mesh%INUM.EQ.1.AND.mesh%ROTSYM.EQ.1)) THEN
      ! suppress x-direction in output
      this%IP1 = 0
    ELSE
      this%IP1 = 1
    END IF
    IF (abs(mesh%ymax-mesh%ymin).LT.2*epsilon(mesh%ymin) &
       .OR.(mesh%JNUM.EQ.1.AND.mesh%ROTSYM.EQ.2)) THEN
      ! suppress y-direction in output
      this%JP1  = 0
    ELSE
      this%JP1  = 1
    END IF
    IF (abs(mesh%zmax-mesh%zmin).LT.2*epsilon(mesh%zmin) &
       .OR.(mesh%KNUM.EQ.1.AND.mesh%ROTSYM.EQ.3)) THEN
      ! suppress z-direction in output
      this%KP1  = 0
    ELSE
      this%KP1  = 1
    END IF

    ! allocate memory for auxilliary data fields
    ALLOCATE(this%binout(maxcomp,                &
                         mesh%IMIN:mesh%IMAX,    &
                         mesh%JMIN:mesh%JMAX,    &
                         mesh%KMIN:mesh%KMAX),   &
             this%vtkcoords(3,                   &
                         mesh%IMIN:mesh%IMAX+this%IP1,  &
                         mesh%JMIN:mesh%JMAX+this%JP1,  &
                         mesh%KMIN:mesh%KMAX+this%KP1), &
             this%output(maxcols),               &
             this%tsoutput(maxcols),             &
             this%bflux(physics%VNUM,6),         &
             stat = err)
    IF (err.NE.0) &
         CALL this%Error("InitFileio_vtk", "Unable to allocate memory.")

    ! determine extend and endianness of real numbers
    CALL this%GetPrecision(this%realsize)
    ! save size of floats to string: realfmt
    WRITE(this%linebuf,fmt='(I4)',iostat=this%error_io) 8*this%realsize
    WRITE(this%realfmt,fmt='(A)',iostat=this%error_io) '"Float' &
          // trim(adjustl(this%linebuf)) // '"'
    CALL this%GetEndianness(this%endianness,'"LittleEndian"','"BigEndian"')

    ! sanity check in case someone ignores the comment in
    ! fileio_common.f90 regarding the KIND type of this variable
    IF (bit_size(this%output(1)%numbytes) .NE. 4*8) &
       CALL this%Error("InitFileio_vtk", "output%numbytes must be a 4 byte integer !!!")

#ifdef PARALLEL
    ! allocate memory for this%extent, send & receive buffers;
    ! recvbuf and this%extent are only needed on the receiver node (rank 0)
    ! but parallel profiling using scalasca with mpich3 crashes
    ! if recvbuf is not allocated on each node hence we allocate a minimal
    ! array of size 1 on all nodes except for the receiver node;
    ! remark: according to the MPI-2 standard recvbuf must only be allocated
    !         on the receiver node
    IF (this%GetRank() .EQ. 0 ) THEN
       n = this%GetNumProcs()
       ALLOCATE(this%extent(0:(n-1)),sendbuf(6),recvbuf(6*n),stat = this%error_io)
    ELSE
       ! see comment above
       ALLOCATE(sendbuf(6),recvbuf(1),stat = this%error_io)
    END IF
    IF (this%error_io.NE.0) &
       CALL this%Error( "InitFileio_vtk", "Unable to allocate memory.")

    ! store information about grid extent on each node
    sendbuf(1) = mesh%IMIN
    sendbuf(2) = mesh%IMAX+this%IP1
    sendbuf(3) = mesh%JMIN
    sendbuf(4) = mesh%JMAX+this%JP1
    sendbuf(5) = mesh%KMIN
    sendbuf(6) = mesh%KMAX+this%KP1

    ! gather information about grid extent on each node at the rank 0 node
    CALL mpi_gather(sendbuf, 6, mpi_integer, recvbuf, 6, mpi_integer, &
               0, mpi_comm_world, this%error_io)

    ! store information about grid extent at the rank 0 node
    IF (this%GetRank() .EQ. 0 ) THEN
       DO i=0,n-1
          WRITE(this%extent(i),fmt='(6(I7))',iostat=this%error_io)&
          recvbuf(i*6+1),recvbuf(i*6+2),recvbuf(i*6+3),recvbuf(i*6+4),recvbuf(i*6+5),recvbuf(i*6+6)
       END DO
    END IF
    ! free buffer memory
    DEALLOCATE(sendbuf,recvbuf,stat=this%error_io)
#endif

    ! check whether mesh coordinates should be cartesian or curvilinear
    IF (this%cartcoords) THEN
       corners => mesh%RemapBounds(mesh%cart%corners)
    ELSE
       corners => mesh%RemapBounds(mesh%curv%corners)
    END IF

    ! store mesh corners in appropriate array for VTK coordinate output
    DO k=mesh%KMIN,mesh%KMAX+mesh%KP1
      DO j=mesh%JMIN,mesh%JMAX+mesh%JP1
         DO i=mesh%IMIN,mesh%IMAX+mesh%IP1
            this%vtkcoords(1:3,i,j,k) = corners(i,j,k,1,1:3)
         END DO
      END DO
    END DO

    ! add coordinates in 1D/2D simulations with non-vanishing extent
    IF (mesh%INUM.EQ.1.AND.this%IP1.EQ.1) THEN
      ! no transport in x-direction, but x-extent > 0
      i = mesh%IMIN
      DO k=mesh%KMIN,mesh%KMAX+mesh%KP1
         DO j=mesh%JMIN,mesh%JMAX+mesh%JP1
            this%vtkcoords(1:3,i+1,j,k) = corners(i,j,k,2,1:3)
         END DO
      END DO
    END IF

    IF (mesh%JNUM.EQ.1.AND.this%JP1.EQ.1) THEN
      ! no transport in y-direction, but y-extent > 0
      j = mesh%JMIN
      DO k=mesh%KMIN,mesh%KMAX+mesh%KP1
         DO i=mesh%IMIN,mesh%IMAX+mesh%IP1
            this%vtkcoords(1:3,i,j+1,k) = corners(i,j,k,3,1:3)
         END DO
      END DO
    END IF

    IF (mesh%KNUM.EQ.1.AND.this%KP1.EQ.1) THEN
      ! no transport in z-direction, but z-extent > 0
      k = mesh%KMIN
      DO j=mesh%JMIN,mesh%JMAX+mesh%JP1
         DO i=mesh%IMIN,mesh%IMAX+mesh%IP1
            this%vtkcoords(1:3,i,j,k+1) = corners(i,j,k,5,1:3)
         END DO
      END DO
    END IF

    ! byte offset, i.e. extent of coordinate data + 4 bytes for storing
    ! the size information which must be a 4 byte integer number
    this%offset = SIZE(this%vtkcoords)*this%realsize + 4

    ! create list of output data arrays
    node => io
    this%cols = 0
    this%tscols = 0
    CALL this%GetOutputlist(mesh,node,this%cols,this%tscols)

    IF (this%GetRank().EQ.0) CALL this%WriteParaviewFile()
  END SUBROUTINE initfileio_vtk

  SUBROUTINE writeparaviewfile(this)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fileio_vtk), INTENT(INOUT) :: this
    !------------------------------------------------------------------------!
    INTEGER                          :: k
    REAL                             :: ftime
    CHARACTER(LEN=256)               :: basename
#ifdef PARALLEL
    INTEGER                          :: i
#endif
    !------------------------------------------------------------------------!
    ! write a pvd-file: this is a "master" file of all timesteps of all ranks
    CALL this%OpenFile(replace,'pvd')

    ! write vtk header
    WRITE(this%unit+200,fmt='(A)',iostat=this%error_io) &
          '<?xml version="1.0"?>'//lf &
          //'<VTKFile type="Collection" version="0.1" byte_order=' &
          //this%endianness//'>'//lf &
          //repeat(' ',2)//'<Collection>'

    ! write entries for each time step
    DO k=0,this%cycles-1
       this%step=k  ! we need this to get the correct file name with time step appended
                    ! ATTENTION: remember to reset the this%step (see below)
       ftime = this%stoptime*k/(this%cycles-1)
#ifdef PARALLEL
       ! in parallel mode each node generates its own file
       DO i=0,this%GetNumProcs()-1
          basename=this%GetBasename(i)
          IF (this%error_io.EQ. 0) WRITE(this%unit+200,fmt='(A,E11.5,A,I4.4,A)',iostat=this%error_io) &
             repeat(' ',4)//'<DataSet timestep="',&
             ftime,'" part="', i ,'" file="'//trim(basename)//'"/>'
       END DO
#else
       basename=this%GetBasename()
       IF (this%error_io.EQ. 0) WRITE(this%unit+200,fmt='(A,E11.5,A)',iostat=this%error_io) &
          repeat(' ',4)//'<DataSet timestep="',ftime,'" part="0" file="' &
                       //trim(basename)//'"/>'
#endif
    END DO
    ! reset the step (see comment above)
    this%step=0

    IF (this%error_io.EQ. 0) WRITE(this%unit+200,fmt='(A)',iostat=this%error_io) &
             repeat(' ',2)//'</Collection>'//lf//'</VTKFile>'
    IF (this%error_io.GT.0) CALL this%Error("InitFileIO_vtk","cannot write pvd-file")

    CALL this%CloseFile('pvd')
  END SUBROUTINE writeparaviewfile


  SUBROUTINE getprecision(this,realsize)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fileio_vtk), INTENT(INOUT) :: this
    INTEGER                          :: realsize
    !------------------------------------------------------------------------!
    CHARACTER(LEN=4)                 :: cTIPO4
    CHARACTER(LEN=8)                 :: cTIPO8
    CHARACTER(LEN=16)                :: cTIPO16
    REAL                             :: rTIPO1, rTIPO2
    INTEGER                          :: k,iTIPO
    !------------------------------------------------------------------------!
    INTENT(OUT)                      :: realsize
    !------------------------------------------------------------------------!

    !endianness
    k = bit_size(itipo)/8
    !HOW big is a REAL?
    rtipo1 = acos(0.0)
    ctipo4 = transfer(rtipo1,ctipo4)
    rtipo2 = transfer(ctipo4,rtipo2)
    IF (rtipo2 == rtipo1) THEN
       realsize = 4
    ELSE
       ctipo8 = transfer(rtipo1,ctipo8)
       rtipo2 = transfer(ctipo8,rtipo2)
       IF (rtipo2 == rtipo1) THEN
          realsize = 8
       ELSE
          ctipo16 = transfer(rtipo1,ctipo16)
          rtipo2 = transfer(ctipo16,rtipo2)
          IF (rtipo2 == rtipo1) THEN
             realsize = 16
          ELSE
             CALL this%Error("GetPrecision_vtk", "Could not estimate size of float type")
          END IF
       END IF
    END IF

  END SUBROUTINE getprecision

  SUBROUTINE getendianness(this, res, littlestr, bigstr)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fileio_vtk), INTENT(INOUT) :: this
    CHARACTER(LEN=*)                 :: res
    CHARACTER(LEN=*)                 :: littlestr
    CHARACTER(LEN=*)                 :: bigstr
    !------------------------------------------------------------------------!
    INTEGER                          :: k,iTIPO
    CHARACTER, POINTER               :: cTIPO(:)
    !------------------------------------------------------------------------!
    INTENT(IN)                       :: littlestr, bigstr
    INTENT(OUT)                      :: res
    !------------------------------------------------------------------------!

    !endianness
    k = bit_size(itipo)/8
    ALLOCATE(ctipo(k),stat = this%error_io)
       IF (this%error_io.NE.0) &
         CALL this%Error("GetEndianness", "Unable to allocate memory.")
    ctipo(1)='A'
    !cTIPO(2:k-1) = That's of no importance.
    ctipo(k)='B'

    itipo = transfer(ctipo, itipo)
    DEALLOCATE(ctipo)
    !Test of 'B'=b'01000010' ('A'=b'01000001')
    IF (btest(itipo,1)) THEN
       write(res,'(A)',iostat=this%error_io)bigstr
    ELSE
       write(res,'(A)',iostat=this%error_io)littlestr
    END IF
  END SUBROUTINE getendianness


  RECURSIVE SUBROUTINE getoutputlist(this,Mesh,node,k,l,prefix)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fileio_vtk), INTENT(INOUT) :: this
    CLASS(mesh_base),  INTENT(IN)    :: mesh
    TYPE(dict_typ),    POINTER       :: node
    INTEGER,           INTENT(INOUT) :: k
    INTEGER,           INTENT(INOUT) :: l
    CHARACTER(LEN=*), OPTIONAL, INTENT(INOUT) :: prefix
    !------------------------------------------------------------------------!
    TYPE(dict_typ), POINTER           :: dir
    CHARACTER(LEN=MAXKEY)             :: key
    TYPE(real_t)                      :: dummy1
    REAL, DIMENSION(:,:), POINTER     :: dummy2
    REAL, DIMENSION(:,:,:), POINTER   :: dummy3
    REAL, DIMENSION(:,:,:,:), POINTER :: dummy4
    INTEGER, DIMENSION(4)             :: dims
    INTEGER                           :: n
    !------------------------------------------------------------------------!
    DO WHILE(ASSOCIATED(node))
      IF(haschild(node)) THEN
      ! recursion
        IF(PRESENT(prefix)) THEN
           key = trim(prefix)//'/'//trim(getkey(node))
        ELSE
           key = '/'//trim(getkey(node))
        END IF
        dir => getchild(node)
        CALL this%GetOutputlist(mesh,dir,k,l,key)
      ELSE IF (hasdata(node).AND.&
               .NOT.(getkey(node).EQ."time")) THEN ! skip time, this is handled in WriteTimestep
         IF (k+1 .GT. maxcols) CALL this%Error("GetOutputlist","reached MAXCOLS")
         ! check shape of output
         dims = 0
         SELECT CASE(getdatatype(node))
         CASE(dict_real_twod)
            CALL getattr(node,trim(getkey(node)),dummy2)
            dims(1:2) = shape(dummy2)
            dims(3:4) = 1
         CASE(dict_real_threed)
            CALL getattr(node,trim(getkey(node)),dummy3)
            dims(1:3) = shape(dummy3)
            dims(4) = 1
         CASE(dict_real_fourd)
            CALL getattr(node,trim(getkey(node)),dummy4)
            dims(1:4) = shape(dummy4)
         CASE(dict_real_p)
            CALL getattr(node,getkey(node),dummy1)
            IF (ASSOCIATED(dummy1%p)) THEN
               l=l+1
               this%tsoutput(l)%val => dummy1%p
               this%tsoutput(l)%key = trim(getkey(node))
               IF (PRESENT(prefix)) THEN
                  this%tsoutput(l)%key = trim(prefix) // '/' // this%tsoutput(l)%key
               END IF
            END IF
         END SELECT
         ! only register output if it contains mesh data
         ! if not => reset k
         ! TODO: HACKATHON 3D - Dritte Meshdimension abfragen.
         IF ((dims(1).EQ.(mesh%IMAX-mesh%IMIN+1)).AND.&
             (dims(2).EQ.(mesh%JMAX-mesh%JMIN+1)).AND.&
             (dims(3).EQ.(mesh%KMAX-mesh%KMIN+1))) THEN
            ! increase output index
            k = k + 1
            ! store name of the output array
            this%output(k)%key = trim(getkey(node))
            ! prepend prefix if present
            IF (PRESENT(prefix)) THEN
               this%output(k)%key = trim(prefix) // '/' // this%output(k)%key
            END IF
            ! allocate memory for pointer to output arrays
            ALLOCATE(this%output(k)%p(dims(4)),stat=this%error_io)
            IF (this%error_io.NE.0) &
               CALL this%Error( "GetOutputlist", "Unable to allocate memory.")
            ! set pointer to output arrays
            IF (dims(4).EQ.1) THEN
               ! 3D data
               this%output(k)%p(1)%val => dummy3
               this%output(k)%numbytes = SIZE(dummy3)*this%realsize
            ELSE
               ! 4D data (for example pvars - density,xvel,yvel)
               DO n=1,dims(4)
                  this%output(k)%p(n)%val => dummy4(:,:,:,n)
                  this%output(k)%numbytes = SIZE(dummy4)*this%realsize
               END DO
            END IF
         END IF
      END IF
      node=>getnext(node)
    END DO
  END SUBROUTINE getoutputlist

  SUBROUTINE openfile(this,action,ftype)
     IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fileio_vtk), INTENT(INOUT)        :: this
    INTEGER,           INTENT(IN)           :: action
    CHARACTER(LEN=*),  INTENT(IN), OPTIONAL :: ftype
    !------------------------------------------------------------------------!
    CHARACTER(LEN=32)  :: sta,act,pos,fext
    !------------------------------------------------------------------------!
    this%error_io=1
    SELECT CASE(action)
    CASE(readonly)
       sta = 'OLD'
       act = 'READ'
       pos = 'REWIND'
    CASE(readend)
       sta = 'OLD'
       act = 'READ'
       pos = 'APPEND'
    CASE(replace)
       sta = 'REPLACE'
       act = 'WRITE'
       pos = 'REWIND'
    CASE(append)
       sta = 'OLD'
       act = 'READWRITE'
       pos = 'APPEND'
    CASE DEFAULT
       CALL this%Error("OpenFile","Unknown access mode.")
    END SELECT

    ! check file type to open
    IF (PRESENT(ftype)) THEN
       fext=trim(ftype)
    ELSE
       ! default
       fext='vts'
    END IF

    SELECT CASE(trim(fext))
    CASE('vts')
       ! open the VTS file
       OPEN(this%unit, file=this%GetFilename(),status=trim(sta), &
            access = 'STREAM' ,   &
            action=trim(act),position=trim(pos),iostat=this%error_io)
#ifdef PARALLEL
    CASE('pvts')
       ! open PVTS file (only parallel mode)
       IF (this%GetRank().EQ.0) THEN
          this%extension='pvts' ! change file name extension
          OPEN(this%unit+100, file=this%GetFilename(-1),status=trim(sta), &
               form='FORMATTED',action=trim(act),position=trim(pos),iostat=this%error_io)
          this%extension='vts' ! reset default extension
       END IF
#endif
    CASE('pvd')
       OPEN(this%unit+200, file=trim(this%filename)//'.'//trim(fext),status=trim(sta), &
            form='FORMATTED',action=trim(act),position=trim(pos),iostat=this%error_io)
    CASE DEFAULT
       CALL this%Error("OpenFile","Unknown file type")
    END SELECT

    IF (this%error_io.GT.0) &
       CALL this%Error("OpenFile","cannot open " // trim(this%extension) // "-file")
  END SUBROUTINE openfile

  SUBROUTINE closefile(this,ftype)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fileio_vtk), INTENT(INOUT) :: this
    CHARACTER(LEN=*), OPTIONAL       :: ftype
    !------------------------------------------------------------------------!
    CHARACTER(LEN=4)                 :: fext
    !------------------------------------------------------------------------!
    INTENT(IN)                       :: ftype
    !------------------------------------------------------------------------!
    this%error_io=1
    ! check file type to open
    IF (PRESENT(ftype)) THEN
       fext=trim(ftype)
    ELSE
       ! default
       fext='vts'
    END IF

    ! close file depending on the file type
    SELECT CASE(trim(fext))
    CASE('vts')
       CLOSE(this%unit,iostat=this%error_io)
#ifdef PARALLEL
    CASE('pvts')
        IF (this%GetRank() .EQ. 0 ) THEN
           CLOSE(this%unit+100,iostat=this%error_io)
        END IF
#endif
    CASE('pvd')
        CLOSE(this%unit+200,iostat=this%error_io)
    CASE DEFAULT
       CALL this%Error("OpenFile","Unknown file type")
    END SELECT
    ! set default extension

    IF(this%error_io .NE. 0) &
       CALL this%Error( "CloseFileIO", "cannot close "//trim(fext)//"-file")
  END SUBROUTINE closefile

!  !> Extract a subkey of a string (obsolete)
!  !!
!  !! Extract a subkey from a string (key) of type "abcd/SUBKEY/efgh"
!  FUNCTION GetSubKey(key) RESULT(subkey)
!    IMPLICIT NONE
!    !------------------------------------------------------------------------!
!    CHARACTER(LEN=*)      :: key
!    !------------------------------------------------------------------------!
!    CHARACTER(LEN=MAXKEY) :: subkey
!    !------------------------------------------------------------------------!
!    INTENT(IN)            :: key
!    !------------------------------------------------------------------------!
!    ! format of key like: "abcd/SUBKEY/efgh"
!    subkey = key(SCAN(key,"/",.TRUE.)+1:)
!  END FUNCTION GetSubKey

  SUBROUTINE writeheader(this,Mesh,Physics,Header,IO)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fileio_vtk),   INTENT(INOUT) :: this
    CLASS(mesh_base),    INTENT(IN)    :: Mesh
    CLASS(physics_base), INTENT(IN)    :: Physics
    TYPE(Dict_TYP),      POINTER       :: IO
    TYPE(Dict_TYP),      POINTER       :: Header
    !------------------------------------------------------------------------!
    CALL this%OpenFile(replace)

    this%error_io = 1
    ! write header in vts files
    WRITE(this%linebuf,fmt='(A,6(I7),A)') '<?xml version="1.0"?>' //lf// &
        '<VTKFile type="StructuredGrid" version="0.1" byte_order='&
                   //this%endianness//'>'//lf//&
        '  <StructuredGrid WholeExtent="',&
                  mesh%IMIN,mesh%IMAX+this%IP1,mesh%JMIN,mesh%JMAX+this%JP1, &
                  mesh%KMIN,mesh%KMAX+this%KP1,'">'//lf// &
!                   Mesh%IMIN,Mesh%IMAX+1,Mesh%JMIN,Mesh%JMAX+1,Mesh%KMIN,Mesh%KMAX+1,'">'//LF// &
        '    <FieldData>'//lf// &
        '      <InformationKey type="String" Name="fosite version" format="ascii">'//lf// &
        '        '//trim(version)//lf// &
        '      </InformationKey>'//lf
    WRITE(this%unit,iostat=this%error_io) trim(this%linebuf)

#ifdef PARALLEL
    ! write header in pvts file (only rank 0)
    IF (this%GetRank() .EQ. 0 ) THEN
       CALL this%OpenFile(replace,'pvts')
       WRITE(this%unit+100,fmt='(A,6(I7),A)',iostat=this%error_io) '<?xml version="1.0"?>' //lf// &
           '<VTKFile type="PStructuredGrid" version="0.1" byte_order='&
                   //this%endianness//'>'//lf//&
           '  <PStructuredGrid WholeExtent="',&
                   1,mesh%INUM+this%IP1,1,mesh%JNUM+this%JP1,1,mesh%KNUM+this%KP1,'">'//lf//&
           '    <PFieldData>'//lf// &
           '      <InformationKey type="String" Name="fosite version" format="ascii">'//lf// &
           '        '//trim(version)//lf//&
           '      </InformationKey>'
       CALL this%CloseFile('pvts')
    END IF
#endif
    IF (this%error_io.GT.0) CALL this%Error("WriteHeader","cannot write (P)VTS file header")

    CALL this%CloseFile()
  END SUBROUTINE writeheader

  SUBROUTINE readheader(this,success)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fileio_vtk), INTENT(INOUT) :: this
    LOGICAL                          :: success
    !------------------------------------------------------------------------!
    INTENT(OUT)                      :: success
    !------------------------------------------------------------------------!
    CALL this%OpenFile(readonly)
    success = .false.
    CALL this%CloseFile()
  END SUBROUTINE readheader

  SUBROUTINE writetimestamp(this,time)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fileio_vtk) :: this
    REAL              :: time
    !------------------------------------------------------------------------!
    INTENT(IN)        :: time
    !------------------------------------------------------------------------!
    CALL this%OpenFile(append)

    WRITE(this%linebuf,fmt='(A,ES14.7,A)') &
        '      <DataArray type="Float64" Name="TIME" NumberOfTuples="1" format="ascii">'//lf// &
                   repeat(' ',8),time,lf//&
        '      </DataArray>'//lf
    WRITE(this%unit,iostat=this%error_io) trim(this%linebuf)
#ifdef PARALLEL
    ! write time stamp in PVTS file in parallel mode
    IF (this%GetRank() .EQ. 0 ) THEN
       CALL this%OpenFile(append,'pvts')
       WRITE(this%unit+100,fmt='(A)',advance='NO',iostat=this%error_io) trim(this%linebuf)
       CALL this%CloseFile('pvts')
    END IF

    CALL this%CloseFile()
#endif

  END SUBROUTINE writetimestamp

  SUBROUTINE readtimestamp(this,time)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fileio_vtk), INTENT(INOUT) :: this
    REAL                             :: time
    !------------------------------------------------------------------------!
    INTENT(OUT)                      :: time
    !------------------------------------------------------------------------!
    CALL this%OpenFile(readonly)
    time = 0.0
    CALL this%CloseFile()
  END SUBROUTINE readtimestamp

  SUBROUTINE writedataset(this,Mesh,Physics,Fluxes,Timedisc,Header,IO)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fileio_vtk),    INTENT(INOUT) :: this
    CLASS(mesh_base),     INTENT(IN)    :: Mesh
    CLASS(physics_base),  INTENT(INOUT) :: Physics
    CLASS(fluxes_base),   INTENT(IN)    :: Fluxes
    CLASS(timedisc_base), INTENT(IN)    :: Timedisc
    TYPE(Dict_TYP),       POINTER       :: Header
    TYPE(Dict_TYP),       POINTER       :: IO
    !------------------------------------------------------------------------!
    INTEGER                             :: k,m
    INTEGER                             :: n, offset
#ifdef PARALLEL
    CHARACTER(LEN=256)                  :: basename
    INTEGER                             :: i
#endif
    !------------------------------------------------------------------------!
    ! this part is formerly from fileio_generic.f90
    ! calculate boundary fluxes, if they were requested for write
!    IF(ASSOCIATED(Timedisc%bflux)) THEN
!        print *, "test"
!        DO k=1,4
!            Timedisc%bflux(:,k) = Fluxes%GetBoundaryFlux(Mesh,Physics,k)
!        END DO
!#ifdef PARALLEL
!        CALL MPI_BCAST(Timedisc%bflux, Physics%VNUM*4, DEFAULT_MPI_REAL, 0, &
!                Mesh%comm_cart, ierror)
!#endif
!    END IF

    IF (ASSOCIATED(timedisc%w)) THEN
      IF (mesh%FARGO.EQ.3.AND.mesh%shear_dir.EQ.1) THEN
        CALL physics%AddBackgroundVelocityX(mesh,timedisc%w,timedisc%pvar,timedisc%cvar)
      ELSE
        CALL physics%AddBackgroundVelocityY(mesh,timedisc%w,timedisc%pvar,timedisc%cvar)
      END IF
    END IF


    ! write the header if either this is the first data set we write or
    ! each data set is written into a new file
    IF ((this%step.EQ.0).OR.(this%cycles.GT.0)) THEN
       CALL this%WriteHeader(mesh,physics,header,io)
    END IF
    ! write the time stamp
    CALL this%WriteTimestamp(timedisc%time)
    CALL this%OpenFile(append)

    !------------------------------------------------------------------------!
    ! REMARK: VTS/PVTS data files are opened in unformated or stream mode,
    !         hence we always write the formated data into a line buffer
    !         and then write the line buffer to the data file
    ! -----------------------------------------------------------------------!
    ! 1. META data
    ! -----------------------------------------------------------------------!
    ! a) global information; currently only scalar real numbers supported
#ifdef PARALLEL
    IF (this%GetRank() .EQ. 0 ) THEN
       CALL this%OpenFile(append,'pvts')
       IF(this%error_io .NE. 0) CALL this%Error( "WriteDataset", "Can't open pvts file")
    END IF
#endif
    ! loop over all time step data registred in GetOutputlist
    DO k=1,this%TSCOLS
       WRITE(this%linebuf,fmt='(A,ES14.7,A)') &
        '      <DataArray type="Float64" Name="'//trim(this%tsoutput(k)%key)// &
                       '" NumberOfTuples="1" format="ascii">'//lf// &
                   repeat(' ',8),this%tsoutput(k)%val,lf//&
        '      </DataArray>'//lf
       WRITE(this%unit,iostat=this%error_io) trim(this%linebuf)
#ifdef PARALLEL
       ! write time step data in PVTS file in parallel mode
       IF (this%GetRank() .EQ. 0 ) THEN
          WRITE(this%unit+100,fmt='(A)',advance='NO',iostat=this%error_io) trim(this%linebuf)
       END IF
#endif
    END DO


! simple hack to write positions of the binary star system
!     IF (HasKey(IO, "/sources/grav/binary/binpos/value"))THEN
!       CALL GetAttr(IO, "/sources/grav/binary/binpos", dummy2)
!
!       WRITE(this%tslinebuf,fmt='(ES14.7,A,ES14.7)',IOSTAT=this%error_io)&
!             dummy2(1,1), repeat(' ',4),dummy2(2,1)
!
!       WRITE(this%unit, IOSTAT=this%error_io)&
!              LF//repeat(' ',6)//'<DataArray type="Float64" Name="position primary"' &
!                            //' NumberOfTuples="2" format="ascii">' &
!            //LF//repeat(' ',8)//TRIM(this%tslinebuf) &
!            //LF//repeat(' ',6)//'</DataArray>'
!
! #ifdef PARALLEL
!       IF (GetRank(this) .EQ. 0 ) &
!          WRITE(this%unit+100, IOSTAT=this%error_io)&
!              LF//repeat(' ',6)//'<DataArray type="Float64" Name="position primary"' &
!                            //' NumberOfTuples="2" format="ascii">' &
!            //LF//repeat(' ',8)//TRIM(this%tslinebuf) &
!            //LF//repeat(' ',6)//'</DataArray>'
!
! #endif
!
!       WRITE(this%tslinebuf,fmt='(ES14.7,A,ES14.7)',IOSTAT=this%error_io)&
!             dummy2(1,2), repeat(' ',4),dummy2(2,2)
!
!       WRITE(this%unit, IOSTAT=this%error_io)&
!              LF//repeat(' ',6)//'<DataArray type="Float64" Name="position secondary"' &
!                            //' NumberOfTuples="2" format="ascii">' &
!            //LF//repeat(' ',8)//TRIM(this%tslinebuf) &
!            //LF//repeat(' ',6)//'</DataArray>'
!
! #ifdef PARALLEL
!       IF (GetRank(this) .EQ. 0 ) &
!          WRITE(this%unit+100, IOSTAT=this%error_io)&
!              LF//repeat(' ',6)//'<DataArray type="Float64" Name="position secondary"' &
!                            //' NumberOfTuples="2" format="ascii">' &
!            //LF//repeat(' ',8)//TRIM(this%tslinebuf) &
!            //LF//repeat(' ',6)//'</DataArray>'
! #endif
!
!     END IF

    ! -----------------------------------------------------------------------!
    ! b) coordinates
    WRITE(this%linebuf,fmt='(A,(6(I7)),A)')&
         '    </FieldData>'//lf// &
         '    <Piece Extent="', &
                    mesh%IMIN,mesh%IMAX+this%IP1,mesh%JMIN,mesh%JMAX+this%JP1,&
                    mesh%KMIN,mesh%KMAX+this%KP1,'">'//lf// &
         '    <Points>'//lf//&
         '      <DataArray type='//this%realfmt//' NumberOfComponents="3"'//&
                         ' Name="Point" format="appended" offset="0">'//lf// &
         '      </DataArray>'//lf// &
         '    </Points>'//lf// &
         '    <CellData>'//lf
    WRITE(this%unit,iostat=this%error_io) trim(this%linebuf)

#ifdef PARALLEL
    IF (this%GetRank() .EQ. 0 ) THEN
       WRITE(this%unit+100,fmt='(A)',iostat=this%error_io)&
              repeat(' ',4)//'</PFieldData>' &
        //lf//repeat(' ',4)//'<PPoints>' &
        //lf//repeat(' ',6)//'<DataArray type='//this%realfmt &
                //' NumberOfComponents="3" Name="Point"/>' &
        //lf//repeat(' ',4)//'</PPoints>' &
        //lf//repeat(' ',4)//'<PCellData>'
    END IF
#endif

    ! -----------------------------------------------------------------------!
    ! c) output data fields
    offset = 0
    DO k = 1, this%cols
       WRITE(this%linebuf,fmt='(A,I3,A,I8,A)')&
           repeat(' ',6)//'<DataArray type='//this%realfmt//&
                          ' NumberOfComponents="', SIZE(this%output(k)%p), &
              '" Name="'//trim(this%output(k)%key)//&
              '" format="appended" offset="',this%offset+offset,'"/>'//lf
       WRITE(this%unit,iostat=this%error_io) trim(this%linebuf)

#ifdef PARALLEL
       IF (this%GetRank() .EQ. 0 ) THEN
           WRITE(this%unit+100,fmt='(A,I3,A)',iostat=this%error_io)&
               repeat(' ',6)//'<DataArray type='//this%realfmt//&
                              ' NumberOfComponents="', SIZE(this%output(k)%p),&
               '" Name="'//trim(this%output(k)%key)//'"/>'
       END IF
#endif
       ! add size of data set (in bytes) + 4 (4 byte integer for size information)
       ! -> "jump address" for next data set in the file in bytes
       offset = offset + this%output(k)%numbytes + 4
    END DO

    WRITE(this%unit,iostat=this%error_io)&
             repeat(' ',6)//'</CellData>' &
       //lf//repeat(' ',4)//'</Piece>' &
       //lf//repeat(' ',2)//'</StructuredGrid>' &
       //lf//repeat(' ',2)//'<GlobalData>'

#ifdef PARALLEL
    IF (this%GetRank() .EQ. 0 ) THEN
      WRITE(this%unit+100,fmt='(A)',iostat=this%error_io) repeat(' ',4)//'</PCellData>'

      DO i=0,this%GetNumProcs()-1
         basename=this%GetBasename(i)
         WRITE(this%unit+100,fmt='(A)',iostat=this%error_io) &
            repeat(' ',4)//'<Piece Extent="'//trim(this%extent(i))&
            //'" Source="'//trim(basename)//'"/>'
      END DO

      WRITE(this%unit+100,fmt='(A)',iostat=this%error_io)&
        repeat(' ',2)//'</PStructuredGrid>' &
        //lf//'</VTKFile>'

    END IF
#endif

    ! -----------------------------------------------------------------------!
    ! d) boundary fluxes
!     DO i=1,4
!        ! in PARALLEL mode the result is sent to process with rank 0
!        this%bflux(:,i) = GetBoundaryFlux(Fluxes,Mesh,Physics,i)
!
!        WRITE(this%unit,IOSTAT=this%error_io)&
!               LF//repeat(' ',4)//'<DataArray type='//this%realfmt//&
!               'Name="'//TRIM(fluxkey(i))//'" format="ascii" >' &
!               //LF//repeat(' ',6)
!        DO k=1,Physics%VNUM
!            WRITE(this%linebuf((k-1)*20+1:k*20),fmt='(E16.9,A)',IOSTAT=this%error_io)&
!                 this%bflux(k,i),repeat(' ',4)
!        END DO
!
!        WRITE(this%unit,IOSTAT=this%error_io)TRIM(this%linebuf(1:Physics%VNUM*20)) &
!               //LF//repeat(' ',4)//'</DataArray>'
!     END DO

    ! end META data
    ! -----------------------------------------------------------------------------------!

    WRITE(this%unit,iostat=this%error_io)&
       lf//repeat(' ',2)//'</GlobalData>' &
       //lf//'<AppendedData encoding="raw">' &
       //lf//'_'

#ifdef PARALLEL
    IF (this%GetRank() .EQ. 0 ) THEN
       CALL this%CloseFile('pvts')
       IF(this%error_io .NE. 0) CALL this%Error( "WriteDataset", "Can't close pvts file")
    END IF
#endif

    ! -----------------------------------------------------------------------!
    ! 2. REAL data
    ! -----------------------------------------------------------------------!
    ! TODO: Hier nochmal rübergucken für 3D
    ! a) coordinates:
    !    (i)  size of data array in bytes (must be a 4 byte integer)
    !    (ii) coordinate data array (generated in InitFileio)
    WRITE(this%unit,iostat=this%error_io) int(this%offset,4),this%vtkcoords(:,:,:,:)
    IF (this%error_io.GT.0) &
       CALL this%Error( "WriteDataset", "cannot write coordinates")

    ! -----------------------------------------------------------------------!
    ! b) data fields
    !    (i)  size of data array in bytes (must be a 4 byte integer);
    !         this is determined in GetOutputlist and storen in
    !         this%output(:)%numbytes
    !    (ii) the data given by a list of pointers which is also
    !         generated in GetOutputlist
    ! loop over all output data fields
    DO k = 1, this%cols
         ! reorganize the data to comply with VTK style ordering
         n = SIZE(this%output(k)%p)  ! get first dimension
         DO m=1,n
            this%binout(m,:,:,:)=this%output(k)%p(m)%val(:,:,:)
         END DO

         WRITE(this%unit,iostat=this%error_io) this%output(k)%numbytes, &
            this%binout(1:n,:,:,:)

         IF (this%error_io.GT.0) &
            CALL this%Error( "WriteDataset", "cannot write data")
    END DO
    WRITE(this%unit,iostat=this%error_io)lf//repeat(' ',2)//&
         '</AppendedData>'//lf//'</VTKFile>'//lf

    CALL this%CloseFile()
    CALL this%IncTime()
  END SUBROUTINE writedataset

  SUBROUTINE finalize(this)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fileio_vtk), INTENT(INOUT) :: this
    !------------------------------------------------------------------------!
    INTEGER                          :: k
    !------------------------------------------------------------------------!
    !------------------------------------------------------------------------!
    DO k=1,this%cols
       IF (ASSOCIATED(this%output(k)%p)) DEALLOCATE(this%output(k)%p)
    END DO
    DEALLOCATE(this%binout,this%vtkcoords,this%output,this%tsoutput,this%bflux)
    CALL this%Finalize_base()
  END SUBROUTINE finalize

END MODULE fileio_vtk_mod