fileio_gnuplot.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: fileio_gnuplot.f90                                                #
!#                                                                           #
!# Copyright (C) 2008-2024                                                   #
!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de>                       #
!# Björn Sperling   <sperling@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_gnuplot_mod
  USE fileio_base_mod
  USE geometry_base_mod
  USE mesh_base_mod
  USE physics_base_mod
  USE timedisc_base_mod
  USE fluxes_base_mod
  USE sources_generic_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
  ! Private Attributes section starts here:
  INTEGER, PARAMETER      :: hlen = 10000
  INTEGER, PARAMETER      :: default_decs = 5
  CHARACTER, PARAMETER    :: sp = achar(32)       
  CHARACTER, PARAMETER    :: lf = achar(10)       
  CHARACTER(LEN=2), PARAMETER  :: recsep = sp // sp    
  CHARACTER(LEN=2), PARAMETER  :: linsep = sp // lf    
  CHARACTER(LEN=2), PARAMETER  :: blksep = lf // lf    
  CHARACTER(LEN=30), PARAMETER :: &
         header_string = "# Data output of fosite" // linsep
  CHARACTER(LEN=HLEN)  :: header_buf
  !--------------------------------------------------------------------------!
  TYPE valptr_typ
    REAL, DIMENSION(:,:,:), POINTER :: val
  END TYPE
 
  TYPE output_typ
    TYPE(valptr_typ), DIMENSION(:), POINTER :: p
    CHARACTER(LEN=MAX_CHAR_LEN)             :: key
    CHARACTER(LEN=1024)                     :: path
    INTEGER                                 :: numbytes
  END TYPE output_typ
 
  TYPE tsoutput_typ
    REAL, POINTER                           :: val
    CHARACTER(LEN=MAX_CHAR_LEN)             :: key
  END TYPE tsoutput_typ
 
  TYPE, EXTENDS(fileio_base) :: fileio_gnuplot
    TYPE(output_typ),DIMENSION(:), POINTER :: &     !< list of output fields
                              output => null()
    TYPE(tsoutput_typ),DIMENSION(:), POINTER :: &   !< list of scalar time step output
                              tsoutput => null()
    CHARACTER(LEN=1), DIMENSION(:,:), POINTER :: &
                              outbuf => null()      
    CHARACTER(LEN=512)     :: heading
    CHARACTER(LEN=512)     :: tsheading
    CHARACTER(LEN=512)     :: linebuf
    CHARACTER(LEN=512)     :: tslinebuf
    CHARACTER(LEN=64)      :: fmtstr
    INTEGER                :: cols
    INTEGER                :: tscols
    INTEGER                :: maxcols
    INTEGER                :: decs
    INTEGER                :: flen
    INTEGER                :: linelen
    INTEGER                :: tslinelen
  CONTAINS
    PROCEDURE :: initfileio_deferred => initfileio_gnuplot
    PROCEDURE :: writeheader
    PROCEDURE :: readheader
    PROCEDURE :: writedataset_deferred => writedataset_gnuplot
!     PROCEDURE :: ReadDataset
    PROCEDURE :: getoutputlist
    final :: finalize
  END TYPE
 
  !--------------------------------------------------------------------------!
  PUBLIC :: &
       ! types
       fileio_gnuplot, output_typ, tsoutput_typ, valptr_typ
  !--------------------------------------------------------------------------!
 
CONTAINS
 
  SUBROUTINE initfileio_gnuplot(this,Mesh,Physics,Timedisc,Sources,config,IO)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fileio_gnuplot),INTENT(INOUT)      :: this
    CLASS(mesh_base),    INTENT(IN)          :: Mesh
    CLASS(physics_base), INTENT(IN)          :: Physics
    CLASS(timedisc_base),INTENT(IN)          :: Timedisc
    CLASS(sources_list), ALLOCATABLE, INTENT(IN) :: Sources
    TYPE(dict_typ),      INTENT(IN), POINTER :: config
    TYPE(dict_typ),      INTENT(IN), POINTER :: IO
    !------------------------------------------------------------------------!
    TYPE(output_typ),DIMENSION(:),POINTER     :: poutput
    TYPE(dict_typ), POINTER                   :: node
    CHARACTER(LEN=MAX_CHAR_LEN),DIMENSION(4)  :: skip
    TYPE(real_t)                              :: dummy1
    REAL, DIMENSION(:,:,:,:), POINTER         :: dummy4
    INTEGER                                   :: cartcoords
    INTEGER                                   :: depth
    INTEGER                                   :: i,j,k,n
#ifdef PARALLEL
    INTEGER                                   :: blocknum
    INTEGER, DIMENSION(Mesh%JMAX-Mesh%JMIN+1,Mesh%KMAX-Mesh%KMIN+1) :: blocklen,indices
#endif
    !------------------------------------------------------------------------!
    CALL this%InitFileio(mesh,physics,timedisc,sources,config,io,"gnuplot","dat",textfile=.false.)
 
    CALL getattr(config, "decimals", this%DECS, default_decs)
    ! compute length of character field for real number output
    ! and check if linebuffer is large enough
    ! flen = 1 (sign) + 1 (one digit) + 1 (decimal point) + decs (decimal places)
    !      + 1 (E character) + 1 (sign of exponent) + 2 (exponent digits) + 2 (spaces)
    this%FLEN = this%DECS + 9
    this%MAXCOLS = len(this%linebuf)/this%FLEN-1
 
    ! this is for registering all arrays which are supposed to be written to the output file
    ALLOCATE(this%output(this%MAXCOLS),this%tsoutput(this%MAXCOLS),stat=this%err)
    IF (this%err.NE.0) &
      CALL this%Error("fileio_gnuplot::InitFileIO","memory allocation failed for this%output, this%tsoutput")
 
    ! get pointer to simulation time to make it the first entry
    ! in the time step data set
    CALL getattr(io,"/timedisc/time",dummy1)
    IF (ASSOCIATED(dummy1%p)) THEN
       this%tsoutput(1)%val => dummy1%p
       this%tsoutput(1)%key = "time"
       this%TSCOLS = 1
    ELSE
       this%TSCOLS = 0
    END IF
 
    ! check if cartesian coordinates are selected for output
    IF (this%cartcoords) THEN
       ! output cartesian coordinates
       CALL getattr(io,"/mesh/bary_centers",dummy4)
       ! requires 3 pointers, one for each cartesian coordinate
       ALLOCATE(this%output(1)%p(3),stat=this%err)
    ELSE
       ! output curvilinear coordinates
       CALL getattr(io,"/mesh/bary_curv",dummy4)
       ! requires NDIMS pointers, one for each curvilinear coordinate, depending on the
       ! dimensionality of the mesh
       ALLOCATE(this%output(1)%p(3),stat=this%err)
      !!!!! only output necessary curvilinear coordinates; needs some addional checks below
!        ALLOCATE(this%output(1)%p(Mesh%NDIMS),STAT=this%err)
    END IF
    IF (this%err.NE.0) &
       CALL this%Error("fileio_gnuplot::InitFileIO","memory allocation failed for this%output(1)%p")
 
    ! register coordinate array pointers for output
    DO n=1,3
      this%output(1)%p(n)%val => dummy4(:,:,:,n)
    END DO
    this%output(1)%key = '# ' // 'x' // repeat(' ',this%FLEN-3) // 'y' // repeat(' ',this%FLEN-1) &
                              // 'z' // repeat(' ',this%FLEN-1)
    
    !!!!!! ATTENTION old code suppresses coordinate output in 1D/2D simulations
    !!!!!! if the coordinate does not vary along a specific dimension
!     IF (Mesh%INUM.EQ.1) THEN
!       ! use format string as temp
!       WRITE (this%fmtstr,'(A5,I2,A1)')'(A1,A',this%FLEN-3,')'
!       WRITE(this%linebuf,TRIM(this%fmtstr))'#','y'
!       this%COLS = 1
!     ELSE IF (Mesh%JNUM.EQ.1) THEN
!       WRITE (this%fmtstr,'(A5,I2,A1)')'(A1,A',this%FLEN-3,')'
!       WRITE(this%linebuf,TRIM(this%fmtstr))'#','x'
!       this%COLS = 1
!     ELSE
!       WRITE (this%fmtstr,'(A5,I2,A2,I2,A1)')'(A1,A',this%FLEN-3,',A',this%FLEN-1,')'
!       WRITE(this%linebuf,TRIM(this%fmtstr))'#','x','y'
!       this%COLS = 2
!     END IF
 
    ! register scalars and arrays for output
    node => io
    skip(1:4) = [CHARACTER(LEN=MAX_CHAR_LEN) :: "bary_centers", "bary_curv", "corners", "time"]
    k = 1
    CALL this%GetOutputList(mesh,node,k,this%TSCOLS,skip)
 
    ! shrink this%output
    poutput => this%output
    ALLOCATE(this%output,source=poutput(1:k),stat=this%err)
    IF (this%err.NE.0) &
      CALL this%Error("fileio_gnuplot::InitFileIO","memory allocation failed for this%output")
    DEALLOCATE(poutput)
 
    ! count number of output columns for array data
    this%COLS = 0
    DO n=1,SIZE(this%output)
       this%COLS = this%COLS + SIZE(this%output(n)%p)
    END DO
 
    ! length of one output line
    this%linelen = this%COLS * this%FLEN
    IF (this%linelen.GT.len(this%linebuf)) &
       CALL this%Error("fileio_gnuplot::InitFileIO", &
          "linebuffer to small; reducing decimals or number of output fields may help")
 
    ! create heading for field data
    this%linebuf = ""
    ! copy heading for coordinates first
    n = SIZE(this%output(1)%p)*this%FLEN
    this%linebuf(1:n-1) = this%output(1)%key(1:n-1)
    ! format string for writing one string of width FLEN
    WRITE(this%fmtstr,'(A,I2,A1)') '(A',this%FLEN,')'
    DO k=2,SIZE(this%output)  ! skip entry for coordinates, i.e. this%output(1)
       i = index(this%output(k)%key,"/",back=.true.) ! find last occurance of slash in key
       ! write key (without leading path)
       WRITE(this%linebuf(n:),trim(this%fmtstr)) this%output(k)%key(i+1:i+this%FLEN-1)
       ! append spaces
       DO i=2,SIZE(this%output(k)%p)
         WRITE(this%linebuf(n+(i-1)*this%FLEN:),trim(this%fmtstr)) repeat(' ',this%FLEN)
       END DO
       n = n + (i - 1)*this%FLEN
    END DO
    this%heading = trim(this%linebuf) // lf
 
    ! create heading for time step data
    ! prepend hash -> comment line in GNUPLOT
    this%linebuf = ""
    DO k=1,this%TSCOLS
       WRITE(this%linebuf(1+(k-1)*this%FLEN:),trim(this%fmtstr)) trim(this%tsoutput(k)%key(1:this%FLEN-1))
    END DO
    this%tsheading = repeat("#",this%FLEN-1) // trim(this%linebuf(1:)) // lf
 
    ! local domain size
    this%INUM = mesh%IMAX - mesh%IMIN + 1
    this%JNUM = mesh%JMAX - mesh%JMIN + 1
    this%KNUM = mesh%KMAX - mesh%KMIN + 1
 
   ! size of the output buffer (on this process)
   this%bufsize  = this%INUM
 
   ! allocate memory for output buffer and displacement records
   this%err = 0
   ALLOCATE(this%outbuf(this%linelen,this%bufsize),stat=this%err)
   IF (this%err.NE.0) &
      CALL this%Error("fileio_gnuplot::InitFileIO_gnuplot","memory allocation failed for this%outbuf")
 
#ifdef PARALLEL
   ! number of output blocks (on this process)
   ! = number of cells (excluding ghost cells)
   blocknum = this%JNUM * this%KNUM
   blocklen(:,:) = this%bufsize
   DO k=1,this%KNUM
      DO j=1,this%JNUM
          indices(j,k) = ((k+mesh%KMIN-2)*mesh%JNUM + (j+mesh%JMIN-2) )*mesh%INUM + mesh%IMIN - 1
      END DO
   END DO
 
   ! create new MPI data types
   SELECT TYPE(df=>this%datafile)
   CLASS IS(filehandle_mpi)
      ! basic type is one output line corresponding to one data point (coordinates + data)
      CALL mpi_type_contiguous(this%linelen,mpi_character,df%basictype,this%err)
      IF (this%err.EQ.0) CALL mpi_type_commit(df%basictype,this%err)
 
      ! file type for blocknum indexed blocks of basic type
      IF (this%err.EQ.0) CALL mpi_type_indexed(blocknum,reshape(blocklen,(/blocknum/)), &
            reshape(indices,(/blocknum/)),df%basictype,df%filetype,this%err)
      IF (this%err.EQ.0) CALL mpi_type_commit(df%filetype,this%err)
   END SELECT
   IF (this%err.NE.0) &
      CALL this%Error("fileio_gnuplot::InitFileIO_gnuplot","creation of MPI file types failed")
 
#endif
    ! write the format string for one entry in the data file:
    ! FLEN-2 characters for the number and 2 for the separators
    WRITE (this%fmtstr,'(A3,I2,A,I2.2,A5)') '(ES', this%FLEN-2, '.', this%DECS,',A2)'
 
    ! print some information
    ! ...
  END SUBROUTINE initfileio_gnuplot
 
 
  RECURSIVE SUBROUTINE getheaderstring(string,root,k,prefix)
  IMPLICIT NONE
  !------------------------------------------------------------------------!
    TYPE(dict_typ), POINTER    :: root,node,subnode
    CHARACTER(LEN=*)           :: string
    CHARACTER(LEN=*), OPTIONAL :: prefix
    CHARACTER(LEN=128)         :: buf
    !------------------------------------------------------------------------!
    INTEGER                    :: idummy, k
    LOGICAL                    :: ldummy
    CHARACTER(LEN=MAX_CHAR_LEN):: cdummy
    REAL                       :: rdummy
    !------------------------------------------------------------------------!
    INTENT(INOUT)              :: string,k
    !------------------------------------------------------------------------!
 
    node => root
    DO WHILE(ASSOCIATED(node))
       SELECT CASE(getdatatype(node))
       CASE(dict_int)
          CALL getattr(node,getkey(node),idummy)
          WRITE(buf,'(A1,A25,I14,A)')'#',trim(getkey(node))//": ",idummy, linsep
          WRITE(string(k:),'(A)')buf
          k = k + len(trim(buf))
       CASE(dict_real)
          CALL getattr(node,getkey(node),rdummy)
          WRITE(buf,'(A1,A25,ES14.4E3,A)')'#',trim(getkey(node))//": ",rdummy, linsep
          WRITE(string(k:),'(A)')buf
          k = k + len(trim(buf))
       CASE(dict_char)
          CALL getattr(node,getkey(node),cdummy)
          WRITE(buf,'(A1,A25,A,A)')'#',trim(getkey(node))//": ",trim(cdummy), linsep
          WRITE(string(k:),'(A)')buf
          k = k + len(trim(buf))
       CASE(dict_bool)
          CALL getattr(node,getkey(node),ldummy)
          WRITE(buf,'(A1,A25,L14,A)')'#',trim(getkey(node))//": ",ldummy, linsep
          WRITE(string(k:),'(A)')buf
          k = k + len(trim(buf))
       END SELECT
       IF(haschild(node)) THEN
          IF (present(prefix)) THEN
             WRITE(buf,'(A)')'#  ['//trim(prefix)//'/'//trim(getkey(node))//']' // linsep
          ELSE
             WRITE(buf,'(A)')'#  ['//trim(getkey(node))//']' // linsep
          END IF
          WRITE(string(k:),'(A)')buf
          k = k + len(trim(buf))
          IF (present(prefix)) THEN
             buf = trim(prefix)//'/'//trim(getkey(node))
          ELSE
             buf = trim(getkey(node))
          END IF
          subnode => getchild(node)
          CALL getheaderstring(string,subnode,k,trim(buf))
       END IF
       node => getnext(node)
    END DO
  END SUBROUTINE getheaderstring
 
  RECURSIVE SUBROUTINE getoutputlist(this,Mesh,node,oarr,onum,skip,prefix)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fileio_gnuplot), INTENT(INOUT) :: this
    CLASS(mesh_base), INTENT(IN)         :: mesh
    TYPE(dict_typ), POINTER              :: node
    INTEGER, INTENT(INOUT)               :: oarr
    INTEGER, INTENT(INOUT)               :: onum
    CHARACTER(LEN=MAX_CHAR_LEN), DIMENSION(:), OPTIONAL, INTENT(IN) &
                                         :: skip
    CHARACTER(LEN=*), OPTIONAL, INTENT(INOUT) &
                                         :: prefix
    !------------------------------------------------------------------------!
    TYPE(dict_typ), POINTER              :: dir
    CHARACTER(LEN=MAX_CHAR_LEN)          :: key
    TYPE(real_t)                         :: dummy1
    REAL, DIMENSION(:,:), POINTER        :: dummy2
    REAL, DIMENSION(:,:,:), POINTER      :: dummy3
    REAL, DIMENSION(:,:,:,:), POINTER    :: dummy4
    REAL, DIMENSION(:,:,:,:,:), POINTER  :: dummy5
    INTEGER, DIMENSION(5)                :: dims
    INTEGER                              :: m,n
    !------------------------------------------------------------------------!
    ! reset error code
    DO WHILE(ASSOCIATED(node))
       ! check if node is directory
       IF(haschild(node)) THEN
          ! recursion
          IF (PRESENT(prefix)) THEN
             ! add prefix to key
             key = trim(prefix)//'/'//trim(getkey(node))
          ELSE
             key = '/'//trim(getkey(node))
          END IF
          dir => getchild(node)
          IF (PRESENT(skip)) THEN
             CALL getoutputlist(this,mesh,dir,oarr,onum,skip,key)
          ELSE
             CALL getoutputlist(this,mesh,dir,oarr,onum,prefix=key)
          END IF
       ELSE IF (hasdata(node)) THEN
          ! node contains data
          IF (PRESENT(skip)) THEN
             IF (any(skip(:) == getkey(node))) THEN
                ! skip entry and continue with next node
                node=>getnext(node)
                cycle
             END IF
          END IF
          ! check shape of data
          dims = 0
          SELECT CASE(getdatatype(node))
!           CASE(DICT_REAL_TWOD)
!               CALL GetAttr(node,TRIM(GetKey(node)),dummy2)
!               dims(1:2) = SHAPE(dummy2)
!               dims(3:5) = 1
          CASE(dict_real_threed)
              CALL getattr(node,trim(getkey(node)),dummy3)
              dims(1:3) = shape(dummy3)
              dims(4:5) = 1
          CASE(dict_real_fourd)
              CALL getattr(node,trim(getkey(node)),dummy4)
              dims(1:4) = shape(dummy4)
              dims(5)   = 1
          CASE(dict_real_fived)
              CALL getattr(node,trim(getkey(node)),dummy5)
              dims(1:5) = shape(dummy5)
          CASE(dict_real_p)
              CALL getattr(node,getkey(node),dummy1)
              IF (ASSOCIATED(dummy1%p)) THEN
                onum=onum+1
                IF (onum.GT.this%MAXCOLS-1) &
                   CALL this%Error("fileio_gnuplot::GetOutputList", &
                            "number of scalar output fields exceeds upper limit")
                ! register number data for output
                this%tsoutput(onum)%val => dummy1%p
                this%tsoutput(onum)%key = trim(getkey(node))
                IF (PRESENT(prefix)) THEN
                    this%tsoutput(onum)%key = trim(prefix) // '/' // this%tsoutput(onum)%key
                END IF
              END IF
          CASE DEFAULT
              CALL this%Warning("fileio_gnuplot::GetOutputList", &
                       "'" // getkey(node) // "'" // " registered for output," &
                       // " but data type is currently not supported")
          END SELECT
          ! only register array data for output if it contains
          ! mesh data
          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
             ! count number of output columns currently registered
             n=0
             DO m=1,oarr
                n = n + SIZE(this%output(m)%p)
             END DO
             ! check limit for output of array data
             IF (n+dims(4)*dims(5).GT.this%MAXCOLS) &
                CALL this%Error("fileio_gnuplot::GetOutputList", &
                         "number of array output fields exceeds upper limit")
             ! increase array output index
             oarr = oarr + 1
             ! store name of the output array
             this%output(oarr)%key = trim(getkey(node))
             ! prepend prefix if present
             IF (PRESENT(prefix)) THEN
                this%output(oarr)%key = trim(prefix) // '/' // this%output(oarr)%key
             END IF
             ! allocate memory for pointer to output arrays
             ALLOCATE(this%output(oarr)%p(dims(4)*dims(5)),stat=this%err)
             IF (this%err.NE.0) &
                CALL this%Error( "fileio_gnuplot::GetOutputList", "Unable to allocate memory.")
             ! set pointer to output arrays
             IF (dims(4).EQ.1) THEN
                ! 3D data
                this%output(oarr)%p(1)%val => dummy3
             ELSE IF (dims(5).EQ.1) THEN
                ! 4D data (for example pvars - density,xvel,yvel)
                DO n=1,dims(4)
                   this%output(oarr)%p(n)%val => dummy4(:,:,:,n)
                END DO
             ELSE
                ! 5D data, e.g., stress tensor components
                DO n=1,dims(5)
                   DO m=1,dims(4)
                      this%output(oarr)%p(n+(m-1)*dims(5))%val => dummy5(:,:,:,m,n)
                   END DO
                END DO
             END IF
          END IF
 
       END IF
       node=>getnext(node)
    END DO
  END SUBROUTINE getoutputlist
 
 
  SUBROUTINE writeheader(this,Mesh,Physics,Header,IO)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fileio_gnuplot), INTENT(INOUT) :: this
    CLASS(mesh_base),      INTENT(IN)    :: Mesh
    CLASS(physics_base),   INTENT(IN)    :: Physics
    TYPE(dict_typ), POINTER              :: Header,IO
    !------------------------------------------------------------------------!
    TYPE(dict_typ), POINTER              :: node
    INTEGER                              :: depth
    !------------------------------------------------------------------------!
    IF (this%GetRank().EQ.0) THEN
       ! get settings from config dictionary and store in header_buf
       depth = 1
       node => header
       CALL getheaderstring(header_buf,node,depth)
       header_buf = trim(header_string) // trim(header_buf)
 
       SELECT TYPE(df=>this%datafile)
#ifndef PARALLEL
       CLASS IS(filehandle_fortran)
          WRITE (unit=df%GetUnitNumber(),iostat=this%err) trim(header_buf) !(1:HLEN-1)
#else
       CLASS IS(filehandle_mpi)
          CALL mpi_file_write(df%GetUnitNumber(),trim(header_buf),len(trim(header_buf)), &
                mpi_character,df%status,this%err)
#endif
       END SELECT
    END IF
  END SUBROUTINE writeheader
 
  SUBROUTINE readheader(this,success)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fileio_gnuplot) :: this
    LOGICAL            :: success
    !------------------------------------------------------------------------!
    INTENT(OUT)        :: success
    !------------------------------------------------------------------------!
#ifdef PARALLEL
    IF (this%GetRank().EQ.0) THEN
    END IF
#else
#endif
    CALL this%Warning("fileio_gnuplot::ReadHeader","reading file header not implemented yet")
    success = .false.
  END SUBROUTINE readheader
 
  SUBROUTINE writedataset_gnuplot(this,Mesh,Physics,Fluxes,Timedisc,Header,IO)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fileio_gnuplot), 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,IO
    INTEGER             :: i,j,k,m,l,n
#ifdef PARALLEL
    INTEGER(KIND=MPI_OFFSET_KIND) :: offset
    INTEGER                       :: request, status(MPI_STATUS_SIZE)
#endif
    !------------------------------------------------------------------------!
    ! some sanity checks:
    ! 1. no output data
    IF (.NOT.ASSOCIATED(this%output)) RETURN
 
    ! 2. uninitilized output data structure
    DO l=1,SIZE(this%output)
       IF (.NOT.ASSOCIATED(this%output(l)%p)) THEN
          CALL this%Error("fileio_gnuplot::WriteDataset", &
                          "this should not happen: output data pointer for " &
                          // trim(this%output(l)%key) // " not associated")
       ELSE
          DO m=1,SIZE(this%output(l)%p)
             IF (.NOT.ASSOCIATED(this%output(l)%p(m)%val)) &
                CALL this%Error("fileio_gnuplot::WriteDataset", &
                                "this should not happen: one of the data array pointers for " &
                                // trim(this%output(l)%key) // " not associated")
          END DO
       END IF
    END DO
 
    this%err = 0
    IF (this%GetRank().EQ.0) THEN
       ! generate string with time step data
       DO k=1,this%TSCOLS
          WRITE(this%tslinebuf(1+(k-1)*this%FLEN:),trim(this%fmtstr)) this%tsoutput(k)%val
       END DO
       this%tslinebuf = repeat("#",this%FLEN-1) // sp // trim(this%tslinebuf) // lf
 
       ! write time step data to file
       SELECT TYPE(df=>this%datafile)
#ifndef PARALLEL
       CLASS IS(filehandle_fortran)
          WRITE (df%GetUnitNumber(),iostat=this%err) trim(this%tsheading) &
              // trim(this%tslinebuf) // trim(this%heading)
#else
       CLASS IS(filehandle_mpi)
          CALL mpi_file_write(df%GetUnitNumber(),trim(this%tsheading),len(trim(this%tsheading)), &
                mpi_character,df%status,this%err)
          IF (this%err.EQ.0) CALL mpi_file_write(df%GetUnitNumber(),trim(this%tslinebuf), &
                len(trim(this%tslinebuf)),mpi_character,df%status,this%err)
          IF (this%err.EQ.0) CALL mpi_file_write(df%GetUnitNumber(),trim(this%heading), &
                len(trim(this%heading)),mpi_character,df%status,this%err)
#endif
       END SELECT
    END IF
 
    IF (this%err.NE.0) &
       CALL this%Error("fileio_gnuplot::WriteDataset","writing time step data to file failed")
 
    SELECT TYPE(df=>this%datafile)
#ifndef PARALLEL
    CLASS IS(filehandle_fortran)
      ! write _one_ line feed at the beginning of each time step
      WRITE (df%GetUnitNumber(),iostat=this%err) lf
#else
    CLASS IS(filehandle_mpi)
      ! very important: wait for the header write command to finish
      CALL mpi_barrier(mpi_comm_world,this%err)
      ! be sure to write at the end, i.e. after the time step data, by getting the offset from the file's size
      IF (this%err.EQ.0) CALL mpi_file_get_size(df%GetUnitNumber(),offset,this%err)
      ! write _one_ line feed at the beginning of each time step
      IF (this%GetRank().EQ.0) THEN
        IF (this%err.EQ.0) CALL mpi_file_write_at(df%GetUnitNumber(), offset, lf, 1, &
             mpi_character, df%status, this%err)
      END IF
#endif
    END SELECT
 
    IF (this%err.NE.0) &
        CALL this%Error("fileio_gnuplot::WriteDataset","writing preceeding line feed to file failed")
 
#ifdef PARALLEL
    SELECT TYPE(df=>this%datafile)
    CLASS IS(filehandle_mpi)
      ! add the initial line feed and the general offset (depends on Mesh%IMIN)
      offset = offset + 1
      ! create the file view
      IF (this%err.EQ.0) CALL mpi_file_set_view(df%GetUnitNumber(),offset,df%basictype,df%filetype, &
          'native',mpi_info_null,this%err)
    END SELECT
 
    IF (this%err.NE.0) &
         CALL this%Error("fileio_gnuplot::WriteDataset","creating MPI file view failed")
#endif
 
    ! trim the floating point values for gnuplot output, i.e. set small numbers to 0
    DO l=1,SIZE(this%output)
      DO m=1,SIZE(this%output(l)%p)
        WHERE (abs(this%output(l)%p(m)%val(:,:,:)).LT.max(tiny(this%output(l)%p(m)%val),1.0d-99))
          this%output(l)%p(m)%val(:,:,:) = 0.0e+00
        END WHERE
      END DO
    END DO
 
    ! write array data to file
    ! do not use Mesh%IMIN/MAX etc., because the indices of this%output()%p()%val start at 1
    this%err = 0
    DO k=mesh%KMIN,mesh%KMAX
      DO j=mesh%JMIN,mesh%JMAX
        DO i=mesh%IMIN,mesh%IMAX
          ! write array data to line buffer
          n = 1
          DO l=1,SIZE(this%output)
            DO m=1,SIZE(this%output(l)%p)
               ! fill output line buffer with data
               WRITE (this%linebuf((n-1)*this%FLEN+1:n*this%FLEN),trim(this%fmtstr)) &
                     this%output(l)%p(m)%val(i-mesh%IMIN+1,j-mesh%JMIN+1,k-mesh%KMIN+1), recsep
               n = n + 1
            END DO
          END DO
 
          IF (mesh%INUM.GT.1) THEN
             IF (i.EQ.mesh%INUM) THEN
                ! finish the block
                this%linebuf(this%linelen-1:this%linelen) = blksep
             ELSE
                ! finish the line
                this%linebuf(this%linelen-1:this%linelen) = linsep
             END IF
          ELSE IF (mesh%JNUM.GT.1) THEN
             IF (j.EQ.mesh%JNUM) THEN
                ! finish the block
                this%linebuf(this%linelen-1:this%linelen) = blksep
             ELSE
                ! finish the line
                this%linebuf(this%linelen-1:this%linelen) = linsep
             END IF
          ELSE
             IF (k.EQ.mesh%KNUM) THEN
                ! finish the block
                this%linebuf(this%linelen-1:this%linelen) = blksep
             ELSE
                ! finish the line
                this%linebuf(this%linelen-1:this%linelen) = linsep
             END IF
          END IF
 
          ! write line buffer to output buffer
          this%outbuf(:,i-mesh%IMIN+1) = transfer(this%linebuf,this%outbuf(:,1),size=this%linelen)
        END DO ! i-loop
 
        ! write output buffer to file
        SELECT TYPE(df=>this%datafile)
#ifndef PARALLEL
        CLASS IS(filehandle_fortran)
           ! write line buffer to file
           IF (this%err.EQ.0) &
              WRITE (df%GetUnitNumber(),iostat=this%err) this%outbuf
#else
        CLASS IS(filehandle_mpi)
          !*****************************************************************!
          ! This collective call didn't work for pvfs2 -> bug in ROMIO ?
          IF (this%err.EQ.0) CALL mpi_file_write_all(df%GetUnitNumber(),this%outbuf, &
                                      this%bufsize,df%basictype,status,this%err)
          !*****************************************************************!
          ! If the collective call above fails try this:
!           IF (this%err.EQ.0) CALL MPI_File_iwrite(df%GetUnitNumber(),this%outbuf, &
!                                      this%bufsize,df%basictype,request,this%err)
!           IF (this%err.EQ.0) CALL MPI_Wait(request,df%status,this%err)
#endif
        END SELECT
 
      END DO ! j-loop
    END DO ! k-loop
 
    IF (this%err.NE.0) &
         CALL this%Error("fileio_gnuplot::WriteDataset","creating MPI file view failed")
 
  END SUBROUTINE writedataset_gnuplot
 
!   !> \public Reads the data arrays from file (not yet implemented)
!   !!
!   SUBROUTINE ReadDataset(this,Mesh,Physics,Timedisc)
!     IMPLICIT NONE
!     !------------------------------------------------------------------------!
!     CLASS(fileio_gnuplot)   :: this     !< \param [in,out] this fileio type
!     CLASS(mesh_base)     :: Mesh     !< \param [in] mesh mesh type
!     CLASS(physics_base)  :: Physics  !< \param [in] physics physics type
!     CLASS(timedisc_base) :: Timedisc !< \param [in] timedisc timedisc type
!     !------------------------------------------------------------------------!
!     INTENT(IN)           :: Mesh,Physics,Timedisc
!     INTENT(INOUT)        :: this
!     !------------------------------------------------------------------------!
!   END SUBROUTINE ReadDataset
 
 
  SUBROUTINE error(this,modproc,msg)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fileio_gnuplot), INTENT(INOUT) :: this
    CHARACTER(LEN=*),  INTENT(IN)     :: modproc
    CHARACTER(LEN=*),  INTENT(IN)     :: msg
    !------------------------------------------------------------------------!
    IF (this%Initialized()) &
      CALL this%datafile%CloseFile(this%step)
    CALL this%Error(modproc,msg)
  END SUBROUTINE error
 
  SUBROUTINE finalize(this)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    TYPE(fileio_gnuplot),INTENT(INOUT) :: this
    !------------------------------------------------------------------------!
    INTEGER                          :: k
    !------------------------------------------------------------------------!
    IF (ASSOCIATED(this%outbuf)) DEALLOCATE(this%outbuf)
    IF (ASSOCIATED(this%output)) THEN
      DO k=1,SIZE(this%output)
        IF (ASSOCIATED(this%output(k)%p)) DEALLOCATE(this%output(k)%p)
      END DO
      DEALLOCATE(this%output)
    END IF
    IF (ASSOCIATED(this%tsoutput)) DEALLOCATE(this%tsoutput)
 
    NULLIFY(this%outbuf,this%output,this%tsoutput)
 
    CALL this%Finalize_base()
  END SUBROUTINE finalize
 
END MODULE fileio_gnuplot_mod