! See copyright notice in the COPYRIGHT file.
! ****************************************************************************** !
!> Module to encapsulate output for various formats.
module hvs_output_module
  use aotus_module,     only: flu_State, aot_get_val, aoterr_Fatal,            &
    &                         aoterr_WrongType, aoterr_NonExistent
  use aot_table_module, only: aot_table_open, aot_table_close

  use env_module,                     only: LabelLen, PathLen, long_k, rk
  use treelmesh_module,               only: treelmesh_type
  use tem_aux_module,                 only: tem_abort
  use tem_comm_env_module,            only: tem_comm_env_type
  use tem_dyn_array_module,           only: PositionOfVal
  use tem_logging_module,             only: logunit
  use tem_subtree_type_module,        only: tem_subtree_type, tem_dump_subTree
  use tem_subTree_module,             only: tem_updatePropertyBits
  use tem_time_module,                only: tem_time_type, tem_time_reset
  use tem_tools_module,               only: upper_to_lower
  use tem_varsys_module,              only: tem_varsys_type
  use tem_varMap_module,              only: tem_varMap_type, tem_create_varMap
  use tem_vrtx_module,                only: tem_vrtx_type, &
    &                                       tem_vrtx_finalize, &
    &                                       tem_calc_vrtx_coord
  use tem_geometry_module,            only: tem_BaryOfId
  ! use tem_transient_reduction_module, only: tem_transient_reduction_type,      &
  !   &                                       tem_transient_reduction_init
  use tem_restart_module,             only: tem_restart_type, tem_init_restart,&
    &                                       tem_restart_dump_data,             &
    &                                       tem_restart_finalize,              &
    &                                       tem_restart_openWrite,             &
    &                                       tem_restart_closeWrite
  use tem_timeControl_module,         only: tem_timeControl_type
  use tem_shape_module,               only: tem_shape_type
  use tem_solveHead_module,           only: tem_solveHead_type


  use hvs_vtk_module, only: hvs_vtk_config_type, hvs_vtk_config_load,          &
    &                       hvs_vtk_file_type, hvs_vtk_open,                   &
    &                       hvs_vtk_write_meshdata, hvs_vtk_close,             &
    &                       hvs_vtk_write_varsys, hvs_vtk_dump_data,           &
    &                       hvs_vtk_init, hvs_vtk_closePVD
  use hvs_ascii_module, only: hvs_ascii_type, hvs_asciiSpatial_type,           &
    &                         hvs_asciiSpatial_open, hvs_ascii_dump_elem_data, &
    &                         hvs_ascii_dump_point_data,                       &
    &                         hvs_asciiSpatial_dump_elem_data,                 &
    &                         hvs_asciiSpatial_dump_point_data,                &
    &                         hvs_ascii_close,                                 &
    &                         hvs_asciiSpatial_close, hvs_ascii_init,          &
    &                         hvs_asciiSpatial_init

  implicit none

  private

  !> parameter for ascii output, writing only one row consisting
  !! of all the elements with all variable values in that row
  !! into an initially opened file
  integer, public, parameter :: hvs_AsciiTransient = 1
  !> Parameter for gnuplot ascii, writing a complete spatial
  !! representation at each time step into a new file
  integer, public, parameter :: hvs_AsciiSpatial   = 2
  !> Parameter for harvester output
  integer, public, parameter :: hvs_Internal       = 3
  !> Parameter for precice spatial coupling
  integer, public, parameter :: hvs_PreciceSpatial = 4
  !> Parameter for VTK output
  integer, public, parameter :: hvs_VTK            = 5

  public :: hvs_output_config_type
  public :: hvs_output_file_type
  public :: hvs_output_load
  public :: hvs_output_init
  public :: hvs_output_open
  public :: hvs_output_write
  public :: hvs_output_close
  public :: hvs_output_finalize

  !> This data type contains data loaded from disk.
  type hvs_output_config_type
    !> Kind of visualization file to use for the output.
    integer :: vis_kind

    !> Description of the vtk configuration (used when vis_kind='vtk')
    type(hvs_vtk_config_type) :: vtk

    !> Logic to decide to use get_point or get_element to dump data
    logical :: useGetPoint

    !> Stores the number of output dofs to be dumped.
    !! @todo KM: 20160425 nDofs is reasonable to configure in output
    !! only for vtk so move to vtk_type.
    integer :: nDofs

  end type hvs_output_config_type

  type hvs_output_file_type
    !> Kind of visualization file to use for the output.
    integer :: vis_kind

    !> Basename to use for the output files.
    character(len=PathLen) :: basename

    !> Process description to use for the output.
    !! Might be only a subset of the global communicator
    type(tem_comm_env_type) :: proc

    !> Description for vtk output
    type(hvs_vtk_file_type) :: vtk

    !> Description for ascii output
    type(hvs_ascii_type) :: ascii

    !> Description for ascii output
    type(hvs_asciiSpatial_type) :: asciiSpatial

    !> Description for harvester output i.e restart format
    type(tem_restart_type) :: restart

    !> Point in time for which this output should be done.
    type(tem_time_type) :: time

    !> Number of variables to write to this file.
    integer :: nVars

    !> List of variable positions to write into this file.
    integer, allocatable :: varPos(:)

    !> Vertex information of elements within the tracking shape.
    !! Required for vtk output
    type(tem_vrtx_type) :: vrtx

    !> Store the barycenters for the linearized tree elements
    !! it has a size of ( nElems, 3 ).
    !! Used in dump AsciiSpatial
    !! It is set  in hvs_output_init
    !! It is used in hvs_output_write
    !! Need to UPDATE after balance
    real(kind=rk), allocatable :: bary(:,:)

    !> is there any transient reduction active?
    !! Then this tracking must be updated every timestep!
    ! logical :: isTransientReduce = .false.

    !> transient reductions which collect and reduce data over several timesteps
    ! type(tem_transient_reduction_type), allocatable :: transientReduce(:)

    !> The number of dofs for each scalar variable of the equation system
    integer :: nDofs

    !> Logic to decide to use get_point or get_element to dump data
    logical :: useGetPoint
  end type hvs_output_file_type


contains


  ! ----------------------------------------------------------------------------!
  !> Read the output configuration from a Lua script.
  subroutine hvs_output_load(me, conf, parent, isReduce)
    ! --------------------------------------------------------------------------!
    !> The output configuration settings to fill.
    type(hvs_output_config_type), intent(out) :: me

    !> Handle of the Lua script to load the configuration from.
    type(flu_state) :: conf

    !> Table handle to the table providing the output settings.
    integer, intent(in), optional :: parent

    !> true if reduction is defined
    logical, intent(in) :: isReduce
    ! ----------------------------------------------------------------------!
    character(len=labelLen) :: vkind
    integer :: thandle
    integer :: iError
    ! ----------------------------------------------------------------------!
    write(logUnit(3), *) '  Loading output table ...'

    call aot_table_open( L       = conf,    &
      &                  parent  = parent,  &
      &                  thandle = thandle, &
      &                  key     = 'output' )

    if (thandle > 0) then
      ! If reduction is active set output%vis_kind to hvs_asciiTransient
      ! else load vis_kind from output table
      if (isReduce) then
        write(logUnit(7),*) 'Spatial reduction is active, set output format to ascii'
        vkind = 'ascii'
      else
        call aot_get_val( L       = conf,     &
          &               thandle = thandle,  &
          &               key     = 'format', &
          &               val     = vkind,    &
          &               ErrCode = iError    )

        if ( btest(iError, aoterr_Fatal) ) then
          write(logUnit(0),*) 'Fatal Error: In reading format for output'
          if ( btest(iError, aoterr_NonExistent) ) then
            write(logUnit(0),*) 'NonExistent: No format for output provided.'
          end if
          if( btest(iError, aoterr_WrongType) ) then
            write(logUnit(0),*) 'WrongType: No format for output provided.'
          end if
          write(logUnit(0),*) 'STOPPING'
          call tem_abort()
        end if
      end if

      vkind = adjustl(vkind)
      vkind = upper_to_lower(vkind)
      select case(trim(vkind))
      case('vtk')
        ! Output data in paraview unstructured vtk format (.vtu)
        me%vis_kind = hvs_VTK
        call hvs_vtk_config_load( me      = me%vtk, &
          &                       conf    = conf,   &
          &                       thandle = thandle )
      case('ascii')
        me%vis_kind = hvs_AsciiTransient
      case('asciispatial')
        ! write an ascii file for each time step and write all three barycenters
        ! each line has corrdinates and variables of only ONE element
        ! #       coordX  coordY  coordZ    var1       var2        var3
        me%vis_kind = hvs_AsciiSpatial
      case('harvester')
        ! ... in the harvester format, meaning
        ! there is a mesh in treelm format and corresponding to it
        ! there are restart files which hold the state information for the
        ! elements in the mesh for the requested time step
        me%vis_kind = hvs_Internal
      case('precice')
        ! ... for preCICE
        ! there is a mesh in treelm format and corresponding to it
        ! there are restart files which hold the state information for the
        ! elements in the mesh for the requested time step
        me%vis_kind = hvs_PreciceSpatial
      case default
        write(logunit(0),*) 'ERROR in hvs_output_load: '     &
          &                 // 'unknown visualization kind ' &
          &                 // trim(vkind)
        write(logunit(0),*) 'format has to be one of:'
        write(logunit(0),*) '* vtk'
        write(logunit(0),*) '* ascii'
        write(logunit(0),*) '* asciiSpatial'
        write(logunit(0),*) '* harvester'
        write(logunit(0),*) '* precice'
        write(logunit(0),*)
        call tem_abort()
      end select
    else
      write(logUnit(0),*) 'FATAL Error: output table is not defined.'
      write(logunit(0),*) 'Please provide an output table with at least the'
      write(logunit(0),*) 'format defined in it!'
      write(logUnit(0),*) 'STOPPING'
      call tem_abort()
    end if
    write(logUnit(5),*) '  Output format: '//trim(vkind)

    ! To decide whether to use get_point or get_element
    call aot_get_val( L       = conf,           &
      &               thandle = thandle,        &
      &               val     = me%useGetPoint, &
      &               ErrCode = iError,         &
      &               default = .false.,        &
      &               key     = 'use_get_point' )
    write(logUnit(7),*) '  Use get_point: ', me%useGetPoint

    ! Get the number of Dofs to be written in the output
    ! The default is set to -1. If the dofs are not specified,
    ! all the dofs should be dumped
    call aot_get_val( L       = conf,     &
      &               thandle = thandle,  &
      &               val     = me%nDofs, &
      &               ErrCode = iError,   &
      &               default = -1,       &
      &               key     = 'ndofs'   )

    call aot_table_close( L       = conf,   &
      &                   thandle = thandle )

  end subroutine hvs_output_load
  ! ----------------------------------------------------------------------------!
  ! ----------------------------------------------------------------------------!


  ! ----------------------------------------------------------------------------!
  !> Initialize the output for a given mesh.
  !!
  !! This creates vertex for a mesh and fill hvs_output_file_type.
  subroutine hvs_output_init(out_file, out_config, tree, varsys, subtree,      &
    &                        varPos, basename, timeControl, nDofs, globProc,   &
    &                        solver, geometry, solSpec_unit)
    ! --------------------------------------------------------------------------!
    !> Output file settings
    !! It must be intent inout since ascii%reduction and trasient%reduction
    !! are loaded in tem_load_trackingHeader
    type(hvs_output_file_type), intent(inout) :: out_file

    !> The output configuration settings to use.
    type(hvs_output_config_type), intent(in) :: out_config

    !> Mesh of the data to visualize.
    type(treelmesh_type), intent(in) :: tree

    !> Description of the available variable system to get the given varnames
    !! from.
    type(tem_varSys_type), intent(in) :: varsys

    !> Optional restriction of the elements to output.
    type(tem_subtree_type), optional, intent(in) :: subtree

    !> List of variable positions that should be written in the output.
    !!
    !! If this is not provided, all variables from the varsys will be written
    !! to the vtk file.
    integer, optional, intent(in) :: varPos(:)

    !> An extension to the output basename.
    !!
    !! The filename will be constructed by tracking%header%prefix//
    !! tracking%header%label
    character(len=*), intent(in) :: basename

    !> output timeControl
    type(tem_timeControl_type), optional, intent(in) :: timeControl

    !> The number of dofs for each scalar variable of the equation system
    integer, intent(in), optional :: nDofs

    !> Global communicator type for global rank information
    type(tem_comm_env_type ), intent(in) :: globProc

    !> Global solver information
    type(tem_solveHead_type ),intent(in) :: solver

    !> shape defined for this ascii output
    type(tem_shape_type), optional, intent(in) :: geometry(:)

    !> Solver specific unit for restart header
    integer, optional, intent(in) :: solSpec_unit
    ! ----------------------------------------------------------------------!
    integer :: iVar, iElem
    integer(kind=long_k) :: glob_nElems, glob_nPoints
    integer :: nElems
    integer :: nPoints
    integer(kind=long_k) :: tTreeID
    ! local varMap for restart init
    type(tem_varMap_type) :: varMap
    ! ----------------------------------------------------------------------!
    ! Copy visualization kind
    out_file%vis_kind = out_config%vis_kind

    out_file%useGetPoint = out_config%useGetPoint

    ! copy basename
    out_file%basename = trim(basename)

    ! nDofs is valid only for get_element
    if (out_file%useGetPoint) then
      out_file%nDofs = 1
    else
      if (present(nDofs)) then
        ! out_config%nDofs is set to -1 if unspecied
        ! in the config file. In this case all the dof's
        ! should be dumped
        if (out_config%nDofs < 0) then
          out_file%nDofs = nDofs
        else
          ! Otherwise the number of dofs dumped should
          ! be what's specified in the config
          out_file%nDofs = out_config%nDofs
        end if
      else
        out_file%nDofs = 1
      end if
    end if

    ! Gather global nElems, local nElems and communicator environment
    if ( present(subTree) ) then
      if (out_file%useGetPoint) then
        nPoints = subTree%nPoints
        glob_nPoints = subTree%glob_nPoints
      else
        nPoints = subTree%nElems
        glob_nPoints = subTree%global%nElems
      end if

      nElems = subTree%nElems
      glob_nElems = subTree%global%nElems
      ! set output communicator
      out_file%proc%comm      = subTree%global%comm
      out_file%proc%rank      = subTree%global%myPart
      out_file%proc%comm_size = subTree%global%nParts
      out_file%proc%root      = 0
    else
      nElems = tree%nElems
      glob_nElems = tree%global%nElems

      ! @todo KM: use nDofs to convert nElems to nPoints
      nPoints = tree%nElems
      glob_nPoints = tree%global%nElems
      ! set output communicator
      out_file%proc = globProc
    end if

    if (allocated(out_file%varpos)) deallocate(out_file%varpos)

    if (present(varPos)) then

      out_file%nVars = size(varPos)
      allocate(out_file%varpos(out_file%nVars))
      out_file%varpos = varpos

    else

      out_file%nVars = varsys%varname%nVals
      allocate(out_file%varpos(out_file%nVars))
      do iVar=1,out_file%nVars
        out_file%varpos(iVar) = iVar
      end do

    end if

    ! ! Init transient reduction
    ! if (out_file%isTransientReduce) then
    !   ! initialize transient reductions
    !   call tem_transient_reduction_init( me     = out_file%transientReduce,    &
    !     &                                nElems = nElems,                      &
    !     &                                varSys = varSys,                      &
    !     &                                varPos = out_file%varPos,             &
    !     &                                nDofs  = out_file%nDofs,              &
    !     &                                time   = timeControl%min              )
    ! end if

    select case(out_file%vis_kind)
    case(hvs_AsciiTransient)
      call hvs_ascii_init( ascii        = out_file%ascii,       &
        &                  varSys       = varSys,               &
        &                  varPos       = out_file%varPos,      &
        &                  basename     = trim(basename),       &
        &                  globProc     = globProc,             &
        &                  nDofs        = out_file%nDofs,       &
        &                  outproc      = out_file%proc,        &
        &                  nElems       = nElems,               &
        &                  glob_nElems  = glob_nElems,          &
        &                  timeControl  = timeControl,          &
        &                  solver       = solver,               &
        &                  useGetPoint  = out_file%useGetPoint, &
        &                  nPoints      = nPoints,              &
        &                  glob_nPoints = glob_nPoints,         &
        &                  geometry     = geometry              )
    case(hvs_AsciiSpatial)
      ! Store barycenter to dump at every time step
      if ( present(subTree) ) then

        if (out_file%useGetPoint) then
          allocate( out_file%bary( nPoints, 3 ) )
          out_file%bary = subTree%points
        else
          allocate( out_file%bary( nElems, 3 ) )
          do iElem = 1, nElems
            tTreeID = tree%treeID( subTree%map2global(iElem) )
            out_file%bary(iElem, :) = tem_BaryOfId(tree, tTreeID)
          end do
        end if

      else

        allocate( out_file%bary( nElems, 3 ) )
        do iElem = 1, nElems
          tTreeID = tree%treeID( iElem )
          out_file%bary(iElem, :) = tem_BaryOfId(tree, tTreeID)
        end do

      end if

      call hvs_asciiSpatial_init( asciiSpatial = out_file%asciiSpatial, &
        &                         varSys       = varSys,                &
        &                         varPos       = out_file%varPos,       &
        &                         basename     = trim(basename),        &
        &                         globProc     = globProc,              &
        &                         outproc      = out_file%proc,         &
        &                         nDofs        = out_file%nDofs,        &
        &                         nElems       = nElems,                &
        &                         glob_nElems  = glob_nElems,           &
        &                         useGetPoint  = out_file%useGetPoint,  &
        &                         nPoints      = nPoints,               &
        &                         glob_nPoints = glob_nPoints,          &
        &                         timeControl  = timeControl,           &
        &                         solver       = solver,                &
        &                         geometry     = geometry               )
    case(hvs_Internal)
      ! -----------------------------------------------------------------------
      ! Initialize restart header, communicator and chunk info
      ! -----------------------------------------------------------------------
      out_file%restart%controller%writePrefix =  trim(basename) // '_'
      out_file%restart%controller%writeRestart = .true.

      ! this restart object is not meant to read data!
      out_file%restart%controller%readRestart = .false.

      ! create varMap for restart
      call tem_create_varMap(varName = varSys%varName%val(out_file%varPos), &
        &                    varSys  = varSys,                              &
        &                    varMap  = varMap                               )

      ! init the restart typed file format
      ! should be called when the first time the mesh is dumped
      call tem_init_restart( me           = out_file%restart, &
        &                    solver       = solver,           &
        ! &                    varSys       = varsys,           &
        &                    varMap       = varMap,           &
        &                    tree         = tree,             &
        &                    subTree      = subTree,          &
        &                    solSpec_unit = solSpec_unit,     &
        &                    nDofs_write  = out_file%nDofs    )

      ! Dump tree if output is harvester format and shape is not global
      ! i.e subTree is present
      if (present(subTree)) then
        call tem_dump_subTree( subTree, tree )
      end if
    case(hvs_VTK)
      ! Calculate vertex for vtk output
      call tem_calc_vrtx_coord( tree    = tree,          &
        &                       vrtx    = out_file%vrtx, &
        &                       subtree = subtree        )

      call hvs_vtk_init( vtk_file   = out_file%vtk,   &
        &                vtk_config = out_config%vtk, &
        &                basename   = trim(basename), &
        &                proc       = out_file%proc   )
    end select

  end subroutine hvs_output_init
  ! ----------------------------------------------------------------------------!
  ! ----------------------------------------------------------------------------!


  ! ----------------------------------------------------------------------------!
  !> Open the output for a given mesh.
  !!
  !! This writes the mesh data to the output files.
  !! Subsequently it can then be enriched by restart information.
  subroutine hvs_output_open(out_file, use_iter, mesh, varsys, time, subtree )
    ! --------------------------------------------------------------------------!
    type(hvs_output_file_type), intent(inout) :: out_file

    !> The output configuration settings to use.
    ! type(hvs_output_config_type), intent(in) :: out_config
    ! Use iteration in filename?
    logical, intent(in) :: use_iter

    !> Mesh of the data to visualize.
    type(treelmesh_type), intent(in) :: mesh

    !> Description of the available variable system to get the given varnames
    !! from.
    type(tem_varSys_type), intent(in) :: varsys

    !> Optional restriction of the elements to output.
    type(tem_subtree_type), optional, intent(inout) :: subtree

    !> Time information.
    !!
    !! If this is present, the filename will be built with a time stamp and
    !! the time point information is written into the vtu file.
    type(tem_time_type), intent(in), optional :: time
    ! ----------------------------------------------------------------------!
    integer :: nElems
    ! ----------------------------------------------------------------------!
    if (present(subtree)) then
      nElems = subtree%nElems
    else
      nElems = mesh%nElems
    end if

    if (present(time)) then
      out_file%time = time
    else
      call tem_time_reset(out_file%time)
    end if

    select case(out_file%vis_kind)
    case(hvs_AsciiSpatial)
      call hvs_asciiSpatial_open( asciiSpatial = out_file%asciiSpatial, &
        &                         outProc      = out_file%proc,         &
        &                         time         = out_file%time,         &
        &                         varsys       = varsys,                &
        &                         varpos       = out_file%varpos,       &
        &                         nDofs        = out_file%nDofs         )
    case(hvs_Internal)
      ! prepare the header differently for using global tree and tracking
      ! subTree. For harvester format, prepare the header.
      if( present(subTree) ) then
        ! For harvester format, prepare the header
        ! Here we pass the tracking subTree and the global tree to the header
        ! In tem_restart_openWrite the subTree will be written to disk

        ! update the property bits according to the global mesh if this has
        ! changed
        call tem_updatePropertyBits( mesh, subTree )

        call tem_restart_openWrite(me      = out_file%restart,             &
          &                        tree    = mesh,                         &
          &                        timing  = out_file%time,                &
          &                        varSys  = varSys,                     &
          &                        subTree = subTree,                      &
          &                        label   = trim(out_file%basename)//'_'  )
      else
        call tem_restart_openWrite( me     = out_file%restart, &
          &                         tree   = mesh,             &
          &                         timing = out_file%time,    &
          &                         varSys = varSys          )
      end if ! present subTree
    case(hvs_VTK)
      call hvs_vtk_open( vtk_file = out_file%vtk,   &
        &                use_iter = use_iter,       &
        &                proc     = out_file%proc,  &
        &                time     = out_file%time   )

      call hvs_vtk_write_meshdata( vtk_file = out_file%vtk,  &
        &                          vrtx     = out_file%vrtx, &
        &                          nElems   = nElems         )

      call hvs_vtk_write_varSys( vtk_file = out_file%vtk,   &
        &                        varsys   = varsys,         &
        &                        varpos   = out_file%varpos )
    end select

  end subroutine hvs_output_open
  ! ----------------------------------------------------------------------------!
  ! ----------------------------------------------------------------------------!


  ! ----------------------------------------------------------------------------!
  subroutine hvs_output_write(out_file, varsys, mesh, subtree)
    ! --------------------------------------------------------------------------!
    type(hvs_output_file_type), intent(inout) :: out_file

    !> Description of the available variable system to get the given varnames
    !! from.
    type(tem_varSys_type), intent(in) :: varsys

    !> Mesh to write the data on.
    type(treelmesh_type), intent(in) :: mesh

    !> Optional restriction of the elements to output.
    type(tem_subtree_type), optional, intent(in) :: subtree
    ! --------------------------------------------------------------------------!
    integer :: iVar
    ! --------------------------------------------------------------------------!

    select case(out_file%vis_kind)
    case(hvs_AsciiTransient)
      if (out_file%useGetPoint) then
        ! Call the get point routine and dump the data at the points
        call hvs_ascii_dump_point_data(                                      &
        &                       ascii           = out_file%ascii,            &
        &                       outProc         = out_file%proc,             &
        &                       varPos          = out_file%varPos,           &
        &                       varSys          = varSys,                    &
        &                       mesh            = mesh,                      &
        &                       time            = out_file%time,             &
        &                       subTree         = subtree )
        ! &                       transientReduce = out_file%transientReduce   )
      else
        ! call the get_element routine and dump all dofs
        call hvs_ascii_dump_elem_data(                                        &
          &                     ascii           = out_file%ascii,             &
          &                     outProc         = out_file%proc,              &
          &                     varPos          = out_file%varPos,            &
          &                     varSys          = varSys,                     &
          &                     mesh            = mesh,                       &
          &                     time            = out_file%time,              &
          &                     nDofs           = out_file%nDofs,             &
          &                     subTree         = subTree )
          ! &                     transientReduce = out_file%transientReduce    )

      endif
    case(hvs_AsciiSpatial)
      if (out_file%useGetPoint) then
        call hvs_asciiSpatial_dump_point_data(                            &
          &                   asciiSpatial    = out_file%asciiSpatial,    &
          &                   varPos          = out_file%varPos,          &
          &                   varSys          = varSys,                   &
          &                   mesh            = mesh,                     &
          &                   time            = out_file%time,            &
          &                   subTree         = subTree,                  &
          ! &                   transientReduce = out_file%transientReduce, &
          &                   bary            = out_file%bary             )
      else
        call hvs_asciiSpatial_dump_elem_data(                             &
          &                   asciiSpatial    = out_file%asciiSpatial,    &
          &                   varPos          = out_file%varPos,          &
          &                   varSys          = varSys,                   &
          &                   mesh            = mesh,                     &
          &                   time            = out_file%time,            &
          &                   nDofs           = out_file%nDofs,           &
          &                   subTree         = subTree,                  &
          ! &                   transientReduce = out_file%transientReduce, &
          &                   bary            = out_file%bary             )
      end if
    case(hvs_Internal)
       call tem_restart_dump_data( restart = out_file%restart,        &
         &                         varSys  = varSys,                  &
         &                         tree    = mesh,                    &
         &                         time    = out_file%time,           &
         &                         subTree = subTree )
         ! &                         transientReduce = out_file%transientReduce )
    case(hvs_VTK)
      do iVar=1,out_file%nVars
        call hvs_vtk_dump_data( vtk_file = out_file%vtk,          &
          &                     varpos   = out_file%varpos(iVar), &
          &                     varSys   = varSys,                &
          &                     mesh     = mesh,                  &
          &                     subtree  = subtree,               &
          &                     time     = out_file%time          )
      end do
    end select

  end subroutine hvs_output_write
  ! ----------------------------------------------------------------------------!
  ! ----------------------------------------------------------------------------!


  ! ----------------------------------------------------------------------------!
  !> Close all files open for current time step.
  subroutine hvs_output_close(out_file, mesh, varSys, subTree)
    ! --------------------------------------------------------------------------!
    type(hvs_output_file_type), intent(inout) :: out_file

    !> Description of the available variable system to get the given varnames
    !! from.
    type(tem_varSys_type), intent(in) :: varsys

    !> Mesh to write the data on.
    type(treelmesh_type), intent(in) :: mesh

    !> Optional restriction of the elements to output.
    type(tem_subtree_type), optional, intent(in) :: subtree
    ! --------------------------------------------------------------------------!
    select case(out_file%vis_kind)
    case(hvs_AsciiSpatial)
      call hvs_asciiSpatial_close( asciiSpatial = out_file%asciiSpatial )
    case(hvs_Internal)
      call tem_restart_closeWrite( me      = out_file%restart, &
        &                          tree    = mesh,             &
        &                          timing  = out_file%time,    &
        &                          varSys  = varSys,         &
        &                          subTree = subTree           )
    case(hvs_VTK)
      call hvs_vtk_close( vtk_file = out_file%vtk, &
        &                 proc     = out_file%proc )
    end select

  end subroutine hvs_output_close
  ! ----------------------------------------------------------------------------!
  ! ----------------------------------------------------------------------------!


  ! ----------------------------------------------------------------------------!
  !> Finalize the output
  subroutine hvs_output_finalize( out_file )
    ! --------------------------------------------------------------------------!
    type(hvs_output_file_type), intent(inout) :: out_file
    ! --------------------------------------------------------------------------!
    select case(out_file%vis_kind)
    case(hvs_AsciiTransient)
      call hvs_ascii_close( ascii = out_file%ascii)
    case(hvs_VTK)
      call hvs_vtk_closePVD(vtk_file = out_file%vtk)
      call tem_vrtx_finalize(out_file%vrtx)
    case(hvs_Internal)
      call tem_restart_finalize( me = out_file%restart )
    end select
  end subroutine hvs_output_finalize
  ! ----------------------------------------------------------------------------!
  ! ----------------------------------------------------------------------------!

end module hvs_output_module
