! See copyright notice in the COPYRIGHT file.
! **************************************************************************** !
!> summary: Module providing datatypes and routines regarding face information.
!! author: Jens Zudrop
!!
!! This module provides information about the existence and properties of
!! faces in an Octree based mesh. It is able to deliver info about the
!! horizontal (i.e. on the same refinement level) and vertical (i.e. between
!! the different levels) dependencies.
!!
!! If you want to make use of face descriptions in your code, you just have to
!! call [[tem_face_module:tem_build_face_info]] to buid up the face
!! descriptions in your solver. You have to make use of the created iterators
!! and face descriptions.
!! The `faces` structure created by [[tem_face_module:tem_build_face_info]]
!! is an array of [[tem_face_type]] with the face information for each level.
!! On each level there are three [[tem_face_descriptor_type]] for the three
!! spatial dimensions `faces(iLevel)%faces(1:3)`. Each of these descriptors
!! holds the complete list of faces on that level for this direction in
!! the `faceList`.
!!
!! In general, the construction of the face description follows the following
!! procedure:
!!
!! 1. Generate level descriptor in a dimension-by-dimension approach.
!! 2. Collect all the faces and attach properties to the faces.
!! 3. Extend the collected properties of the remote elements by from finer or
!!    from coarser properties (to indicate if a halo element is from coarser or
!!    from finer).
!! 4. Create vertical dependencies between the faces.
!! 5. Create list of faces which are relevant for computation and interpolation.
!! 6. Create communication buffers for send and receive of face information.
!!
!! Each face has two neighboring elements. To identify each face uniquely, we
!! use the tree ID of its left element.
!!```
!!       ---------------------------------------
!!      /                  /                  /
!!     /                  /                  /
!!    /   left elem.     /   right elem.    /
!!   /                  /                  /
!!  /                  /                  /
!!  --------------------------------------
!!```
!!
!! This module provides a basic description of the faces per spatial direction.
!! Building this description is achieved in several steps:
!!
!! 1. Collect all the faces (each level independently).
!! 2. Determine the properties of the faces.
!!    a. Check if halo elements are refined on the remote partition.
!! 3. Build the vertical dependencies between the faces of the different levels.
!! 4. Build buffers and lists of faces from the previous steps.
!!
!! To define a correct list of faces for the solvers, we attach left and right
!! properties to each face (in part 1 of the upper algorithm). We differentiate
!! between a left and right property of each face. The left one is determined
!! by it's left element and the right one by it's right element respectively.
!!```
!!              tem_left / / / tem_right
!!       -----------------------------------------------------
!!      /                      / / /                        /
!!     /                      / / /                        /
!!    /   left elem. --->    / / /    <--- right elem.    /
!!   /                      / / /                        /
!!  /                      / / /                        /
!!  ----------------------------------------------------
!!                       / / /
!!```
!!
!! We use one or more of the following properties to the left and right
!! properties of the faces:
!!
!! a. [[tem_faceData_module:tem_fluidface_prp]]
!! b. [[tem_faceData_module:tem_fromfinerface_prp]]
!! c. [[tem_faceData_module:tem_fromcoarserface_prp]]
!! d. [[tem_faceData_module:tem_remoteface_prp]]
!!
!! These properties are the basis for building the correct list of faces in step
!! 4 of the above algorithm.
!!
!! See [[Face Details]] for examples on face configurations.

module tem_face_module

  ! include treelm modules
  use env_module,              only: long_k
  use tem_dyn_array_module,    only: dyn_longArray_type, dyn_intArray_type, &
    &                                init, append, PositionOfVal
  use tem_grow_array_module,   only: grw_intArray_type, grw_longArray_type, &
    &                                init, append
  use tem_construction_module, only: tem_find_allElements,             &
    &                                tem_init_elemLevels,              &
    &                                tem_build_VerticalDependencies,   &
    &                                tem_build_HorizontalDependencies, &
    &                                tem_cleanupDependencyArrays,      &
    &                                tem_levelDesc_type, tem_treeIdInTotal
  use tem_param_module,        only: q__E, q__N, q__T, q__W, q__S, q__B, &
    &                                qOffset, qInvDir
  use tem_topology_module,     only: tem_directChildren, tem_ParentOf, &
    &                                tem_coordOfID, tem_IDofCoord
  use tem_logging_module,      only: logUnit
  use tem_aux_module,          only: tem_abort
  use tem_geometry_module,     only: tem_eligibleChildren
  use tem_stencil_module,      only: tem_stencilHeader_type, &
    &                                tem_stencil_dump
  use treelmesh_module,        only: treelmesh_type
  use tem_bc_prop_module,      only: tem_bc_prop_type
  use tem_comm_module,         only: tem_commpattern_type,tem_communication_type
  use tem_comm_env_module,     only: tem_comm_env_type
  use tem_faceData_module,     only: tem_face_type, tem_face_descriptor_type, &
    &                                tem_left, tem_right, &
    &                                tem_faceList_type,                   &
    &                                tem_notExist_prp, tem_fluidFace_prp, &
    &                                tem_fromFinerFace_prp,               &
    &                                tem_fromCoarserFace_prp,             &
    &                                tem_remoteFace_prp, tem_bndFace_prp
  use tem_compteFaceRules_module, only: tem_isComputeFace
  use tem_element_module,    only: eT_fluid, eT_halo,  &
    &                              eT_ghostFromCoarser, eT_ghostFromFiner
  use tem_spacetime_fun_module,  only: tem_st_fun_listElem_type

  implicit none

  private

  public :: tem_build_face_info


contains


  ! ************************************************************************** !
  !> summary: Routine to get face information for your whole mesh.
  !!
  !! This routine builds a description of the faces on each level. The
  !! resulting array of [[tem_face_type]] ranges from the minimal level
  !! of the mesh to the maximal level of the mesh.
  subroutine tem_build_face_info( tree, boundary, commPattern, proc, faces, &
                                & nEligibleChildren )
    ! --------------------------------------------------------------------------
    !> Tree representation of the mesh to buid the faces for.
    type(treelmesh_type), intent(inout) :: tree

    !> The boundaries of your simulation domain.
    type(tem_bc_prop_type), intent(in) :: boundary

    !> The communication pattern you use for the buffer.
    type(tem_commpattern_type), intent(in) :: commPattern

    !> Process description to use.
    type(tem_comm_env_type), intent(in) :: proc

    !> The number of eligible children for the vertical face dependency
    integer, intent(in) :: nEligibleChildren

    !> The created face descriptors (one for each level).
    type(tem_face_type), allocatable, intent(out) :: faces(:)
    ! --------------------------------------------------------------------------
    ! Level descriptor for each spatial direction and each level of your mesh.
    ! The level descriptor have to be constructed with the dimension by
    ! dimension stencils (+1, 0, -1) for each spatial direction.
    type(tem_levelDesc_type) :: levelDesc( 1:3, tree%global%minLevel  &
      &                                         :tree%global%maxLevel )
    type(tem_levelDesc_type), allocatable :: levelDescX(:)
    type(tem_levelDesc_type), allocatable :: levelDescY(:)
    type(tem_levelDesc_type), allocatable :: levelDescZ(:)
    integer :: iLevel
    ! --------------------------------------------------------------------------


    allocate(faces(tree%global%minLevel:tree%global%maxLevel))

    ! create the level desriptor in a dim-by-dim manner.
    write(logUnit(3),*)'Building dim-by-dim level descriptors for faces ...'
    call tem_dimByDim_construction( tree, boundary, commPattern, proc, &
      &                             levelDescX, levelDescY, levelDescZ )
    do iLevel = tree%global%minLevel,tree%global%maxLevel
      levelDesc(1,iLevel) = levelDescX(iLevel)
      levelDesc(2,iLevel) = levelDescY(iLevel)
      levelDesc(3,iLevel) = levelDescZ(iLevel)
    end do
    deallocate(levelDescX, levelDescY, levelDescZ)

    ! First we collect all the faces. A face is always identified by its left
    ! element tree id.
    ! We collect all the faces, even if only one element next to the face
    ! exists.
    ! Therefore, we run over fluids, ghostFromCoarser, ghostFromFiner and halo
    ! elements.
    ! We attach two-sided properties to the faces, too (depending on the element
    ! left and right to the face).
    call tem_collect_faces( levelDesc, tree%global%minLevel, &
      &                     tree%global%maxLevel, faces      )

    ! Now, we have to check all the faces which have the communication property
    ! again.
    ! We have to check if the halo element next to it is refined/from coarser on
    ! the remote rank. So, we extend the faces' properties.
    call tem_extend_remotePrp( levelDesc, tree%global%minLevel, &
      &                        tree%global%maxLevel, faces      )

    ! Now, we build the vertical dependency of the faces.
    call tem_faceDep_vertical( tree%global%minLevel, tree%global%maxLevel, &
      &                        faces, nEligibleChildren                    )

    !> \todo Remove all the unnecessary faces from the face descriptor

    ! Build the buffers for sent/received faces.
    call tem_build_faceBuffers( tree%global%minLevel, tree%global%maxLevel, &
      &                         faces, levelDesc                            )

    ! Build the iterators for the faces. These ones will be used by the solvers.
    call tem_build_faceLists(tree%global%minLevel, tree%global%maxLevel, &
      &                       nEligibleChildren, faces                   )

    !> \todo JZ: this is very memory consuming since we store them in different
    !! arrays several times.
    ! Buffer the dimension by dimension descriptor
    do iLevel = tree%global%minLevel,tree%global%maxLevel
      faces(iLevel)%dimByDimDesc(1) = levelDesc(1,iLevel)
      faces(iLevel)%dimByDimDesc(2) = levelDesc(2,iLevel)
      faces(iLevel)%dimByDimDesc(3) = levelDesc(3,iLevel)
    end do

  end subroutine tem_build_face_info
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> summary: Subroutine to build communication buffers for the faces.
  !!
  subroutine tem_build_faceBuffers( minLevel, maxLevel, faces, levelDesc )
    ! --------------------------------------------------------------------------
    !> The minimum refinement level of your mesh.
    integer, intent(in) :: minLevel
    !> The maximum refinement level of your mesh.
    integer, intent(in) :: maxLevel
    !> The communication pattern you want use for the buffer.
    ! type(tem_commpattern_type), intent(in) :: commPattern
    !> The created face descriptor.
    type(tem_face_type),intent(inout)  :: faces(minLevel:maxLevel)
    !> Dimension-by-dimension level descriptors
    type(tem_levelDesc_type), intent(in) :: levelDesc(1:3,minLevel:maxLevel)
    ! --------------------------------------------------------------------------
    integer :: iLevel, iDir
    ! --------------------------------------------------------------------------

    ! Create the face buffers for which I will receive information before
    ! my rank can do the timestep. Attention: this is consistent with the
    ! definition of a compute face.
    do iLevel = minLevel, maxLevel
      do iDir = 1, 3
        ! buffer for the state
        call tem_build_faceRecvBuffers(                                 &
          &    faces       = faces(iLevel)%faces(iDir),                 &
          &    levelDesc   = levelDesc(iDir, iLevel),                   &
          &    buf         = faces(iLevel)%faces(iDir)%recvBuffer_state )

        ! buffer for the flux
        call tem_build_faceRecvBuffers(                                &
          &    faces       = faces(iLevel)%faces(iDir),                &
          &    levelDesc   = levelDesc(iDir, iLevel),                  &
          &    buf         = faces(iLevel)%faces(iDir)%recvBuffer_flux )
      end do
    end do

    ! Create the face buffers for which I will send information before
    ! the other ranks can the timestep.
    do iLevel = minLevel, maxLevel
      do iDir = 1, 3
      ! buffer for the state
      call tem_build_faceSendBuffers(                                 &
        &    faces       = faces(iLevel)%faces(iDir),                 &
        &    levelDesc   = levelDesc(iDir, iLevel),                   &
        &    buf         = faces(iLevel)%faces(iDir)%sendBuffer_state )
      ! buffer for the flux
      call tem_build_faceSendBuffers(                                &
        &    faces       = faces(iLevel)%faces(iDir),                &
        &    levelDesc   = levelDesc(iDir, iLevel),                  &
        &    buf         = faces(iLevel)%faces(iDir)%sendBuffer_flux )
      end do
    end do

  end subroutine tem_build_faceBuffers
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> summary: Subroutine to build communication buffers for the faces the
  !! process will receive information for (before it can make the compute step).
  subroutine tem_build_faceRecvBuffers( levelDesc, faces, buf )
    ! -------------------------------------------------------------------------
    !> Dimension-by-dimension level descriptor for the current level and
    !! direction.
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> The communication pattern you want use for the buffer.
    ! type(tem_commpattern_type), intent(in) :: commPattern
    !> The created face descriptor.
    type(tem_face_descriptor_type),intent(inout)  :: faces
    !> The created receive buffer. Size is two, due to left and right limes of
    !! the face values. To access it use tem_left and tem_right.
    type(tem_communication_type), intent(out) :: buf(2)
    ! --------------------------------------------------------------------------
    integer :: iFace, iProc, iSide, elemPos, faceSide, tIdPos, rank, rankPos,  &
      &        elemAddPos, faceIndex
    logical :: wasAdded
    integer(kind=long_k) :: elemId
    ! The elements ids I will receive data for. We have to differentiate the
    ! left and right face of a face here.
    type(dyn_longArray_type) :: elemIds(2)
    ! The elements positions I will receive data for. We have to differentiate
    ! the left and right face of a face here.
    type(grw_intArray_type) :: elemPositions(2)
    ! The source ranks.
    type(dyn_intArray_type) :: srcRank(2)
    ! The position of the source ranks.
    type(grw_intArray_type) :: srcRankPos(2)
    ! --------------------------------------------------------------------------

    ! Init our growing arrays.
    call init( me = elemIds(1) , length = 16 )
    call init( me = elemIds(2) , length = 16 )
    call init( me = elemPositions(1) , length = 16 )
    call init( me = elemPositions(2) , length = 16 )
    call init( me = srcRankPos(1) , length = 16)
    call init( me = srcRankPos(2) , length = 16)
    call init( me = srcRank(1) , length = 16 )
    call init( me = srcRank(2) , length = 16 )

    ! Iterate over all the faces and count the communicate faces.
    do iFace = 1, faces%faceList%faceId%nVals

      ! Check whether we will receive info of one of the adjacent elements
      ! of this face.
      faceSide = tem_isRecvFace(iFace, faces)
      if (faceSide /= 0) then

        ! Check for which adjacent element (left or right) the current rank
        ! will receive information about and count it.
        if (faceSide == tem_left) then
          elemPos = faces%faceList%leftElemPos%val(iFace)
          elemId = faces%faceList%faceId%val(iFace)
          ! If we receive information about the left element, we receive
          ! information about it's right face!
          faceSide = tem_right
        else
          elemPos = faces%faceList%rightElemPos%val(iFace)
          elemId = faces%faceList%rightElemId%val(iFace)
          ! If we receive information about the right element, we receive
          ! information about it's left face!
          faceSide = tem_left
        end if

        ! Get the source rank. Therefore we have a look at the level
        ! descriptor of the current spatial direction and level.
        tIdPos = PositionOfVal( me = levelDesc%elem%tId, val = elemId )
        if (tIdPos <= 0) then
          write(*,*) 'ERROR in tem_build_faceRecvBuffers: not able to '//  &
            &        'identify source rank for a certain face, stopping ...'
          call tem_abort()
        else
          rank = levelDesc%elem%sourceProc%val(tIdPos)
        end if

        ! Store element position and source proc in a dynamic array.
        call append( me       = elemIds(faceSide), &
          &          val      = elemId,            &
          &          pos      = elemAddPos,        &
          &          wasAdded = wasAdded           )
        if (wasAdded) then
          call append( me = elemPositions(faceSide), val = elemPos )
          call append( me = srcRank(faceSide), val = rank, pos = rankPos )
          call append( me = srcRankPos(faceSide), val = rankPos )
        end if

      end if
    end do

    ! Now, we build the buffer for the faces (left and right values).
    do iSide = 1, 2
      ! Init the datastructures in the buffer
      buf(iSide)%nProcs = srcRank(iSide)%nVals
      allocate(buf(iSide)%proc( srcRank(iSide)%nVals ))
      allocate(buf(iSide)%nElemsProc( srcRank(iSide)%nVals ))
      allocate(buf(iSide)%elemPos( srcRank(iSide)%nVals ))
      do iProc = 1, buf(iSide)%nProcs
        call init( me = buf(iSide)%elemPos(iProc), length = 16 )
      end do

      ! Set the positions of the elements I will receive.
      do iFace = 1, elemPositions(iSide)%nVals
        faceIndex = elemIds(iSide)%sorted(iFace)
        elemPos = elemPositions(iSide)%val(faceIndex)
        rankPos = srcRankPos(iSide)%val(faceIndex)
        call append( me = buf(iSide)%elemPos(rankPos), val = elemPos )
      end do

      ! Count the number of elements we recv from each proc.
      do iProc = 1, buf(iSide)%nProcs
        buf(iSide)%proc(iProc) = srcRank(iSide)%val(iProc) - 1
        buf(iSide)%nElemsProc(iProc) = buf(iSide)%elemPos(iProc)%nVals
      end do
    end do

  end subroutine tem_build_faceRecvBuffers
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> summary: Subroutine to check if the current rank will receive information
  !! for the a given face before the compute step. Since face data
  !!
  pure function tem_isRecvFace(facePos, faces ) result( isRecvFace )
    ! --------------------------------------------------------------------------
    !> The position of the face in faces to be checked.
    integer , intent(in) :: facePos
    !> The description of the faces.
    type(tem_face_descriptor_type),intent(in)  :: faces
    !> Returns 0 if no information has to be received from remote ranks
    !! before a compute step can be done. In case that information has to be
    !! received it will return a non-zero positive value. Furhtermore,
    !! it determines whether information for the left face value
    !! (left adjacent element)
    !! or the right face value (right adjacent element) has to be received.
    !! In the first case it will return
    !! [[tem_faceData_module:tem_left]] and in the second case
    !! it will return
    !! [[tem_faceData_module:tem_right]].
    !! If something went wrong, e.g. an unknown combination occurs,
    !! it will return -1.
    integer :: isRecvFace
    ! --------------------------------------------------------------------------
    integer :: leftPrp, rightPrp
    ! --------------------------------------------------------------------------

    ! get the left and right property of the faces.
    leftPrp = faces%faceList%leftPrp%val(facePos)
    rightPrp = faces%faceList%rightPrp%val(facePos)

    isrecvface = -1 ! default value (error)

    ! So, check for all conditions under which face info has to be received.
    ! Fluid-fluid face (purely local)
    ! -> nothing to be sent
    if ( (leftPrp == tem_fluidFace_prp)        &
      &  .and. (rightPrp == tem_fluidFace_prp) ) then
      isRecvFace = 0

    ! Fluid-remote face (the halo is not refined)
    ! -> left element is on my rank, so I will compute this face. Therefore,
    !    I will also receive the remote information. In this case I will
    !    receive info from the right adjacent element.
    elseif ( (leftPrp == tem_fluidFace_prp)        &
      &     .and. (rightPrp == tem_remoteFace_prp) ) then
      isRecvFace = tem_right

    ! Remote-fluid face (the halo is not refined)
    ! -> left element is NOT on my rank, so another rank will compute this face.
    !    So, I have to send info for the left adjacent element, but I do not
    !    not receive information.
    elseif( (leftPrp == tem_remoteFace_prp)       &
      &     .and. (rightPrp == tem_fluidFace_prp) ) then
      isRecvFace = 0

    ! Remote-notExisting face (the remote is not refined)
    ! -> remote is not located on my rank and the other element does not
    !    exist. This face does not have to be computed on any rank, so my rank
    !    is not recv any info about any element adjacent to this face.
    elseif( leftPrp.eq.tem_remoteFace_prp .and.                                &
      &     rightPrp.eq.tem_notExist_prp ) then
      isRecvFace = 0

    ! notExisting-remote face (the remote is not refined)
    ! -> remote is not located on my rank and the other element does not
    !    exist. This face does not have to be computed on any rank, so my rank
    !    is not recv any info about any element adjacent to this face.
    elseif( leftPrp.eq.tem_notExist_prp .and.                                  &
      &     rightPrp.eq.tem_remoteFace_prp ) then
      isRecvFace = 0

    ! Fluid-fromFiner face (purely local)
    ! -> everything is local, so there is no need to send or recv information
    !    from another rank.
    elseif( leftPrp.eq.tem_fluidFace_prp .and.                                 &
      &     rightPrp.eq.tem_fromFinerFace_prp ) then
      isRecvFace = 0

    ! FromFiner-fluid face (purely local)
    ! -> everything is local, so there is no need to send or recv information
    !    from another rank.
    elseif( leftPrp.eq.tem_fromFinerFace_prp .and.                             &
      &     rightPrp.eq.tem_fluidFace_prp ) then
      isRecvFace = 0

    ! FromFiner-FromFiner face (purely local)
    ! -> everything is local. Somehow two ghost elements meet each other, but
    !    we do not have to recv or send any data for these faces.
    elseif( leftPrp.eq.tem_fromFinerFace_prp .and.                             &
      &     rightPrp.eq.tem_fromFinerFace_prp ) then
      isRecvFace = 0

    ! FromCoarser-fluid face (purely local)
    ! -> everything is local, so no data has to be send or received.
    elseif( leftPrp.eq.tem_fromCoarserFace_prp .and.                           &
      &     rightPrp.eq.tem_fluidFace_prp ) then
      isRecvFace = 0

    ! Fluid-fromCoarser face (purely local)
    ! -> everything is local, so no data has to be send or received.
    elseif( leftPrp.eq.tem_fluidFace_prp .and.                                 &
      &     rightPrp.eq.tem_fromCoarserFace_prp ) then
      isRecvFace = 0

    ! FromCoarser-fromCoarser face (purely local)
    ! -> everything is local, so no data has to be send or received.
    elseif( leftPrp.eq.tem_fromCoarserFace_prp .and.                           &
      &     rightPrp.eq.tem_fromCoarserFace_prp ) then
      isRecvFace = 0

    ! Remote-fromCoarser face
    ! -> My rank holds the right elements. Therefore, my element is on a coarser
    !    level.
    !    However, the remote elements are refined, but not located on a single
    !    rank (otherwise it would have the notExist property instead of the
    !    remoteFace property). Therefore, my rank will send data to the remote
    !    ranks (as they are on a finer level), but we do not receive any info.
    elseif( leftPrp.eq.tem_remoteFace_prp .and.                                &
      &     rightPrp.eq.tem_fromCoarserFace_prp ) then
      isRecvFace = 0

    ! FormCoarser-Remote face
    ! -> My rank holds the left element. Therefore, my element is on a coarser
    !    level. However, the remote elements are refined, but not located on a
    !    single rank (otherwise it would have the notExist property instead of
    !    the remoteFace property). Therefore, my rank will send data to the
    !    remote ranks (as they are on a finer level), but we do not receive any
    !    info.
    elseif( leftPrp.eq.tem_fromCoarserFace_prp .and.                           &
      &     rightPrp.eq.tem_remoteFace_prp ) then
      isRecvFace = 0

    ! Remote-Remote face
    ! -> My rank holds none of the elements located next to this face. So we do
    !    not receive any information for this face.
    elseif( leftPrp.eq.tem_remoteFace_prp .and.                                &
      &     rightPrp.eq.tem_remoteFace_prp ) then
      isRecvFace = 0

    ! notExist - remote face
    elseif( leftPrp.eq.tem_notExist_prp .and.                                  &
      &     rightPrp.eq.tem_remoteFace_prp) then
      isRecvFace = 0

    ! remote - notExist face
    elseif( leftPrp.eq.tem_remoteFace_prp .and.                                &
      &     rightPrp.eq.tem_notExist_prp ) then
      isRecvFace = 0

    ! fluid-remote+fromCoaser face
    elseif( leftPrp.eq.tem_fluidFace_prp .and.                                 &
      &     rightPrp.eq.ior( tem_remoteFace_prp, tem_fromCoarserFace_prp )) then
      isRecvFace = tem_right

    ! remote+fromCoarser-fluid face
    elseif( leftPrp.eq.ior(tem_remoteFace_prp,tem_fromCoarserFace_prp) .and.   &
      &     rightPrp.eq.tem_fluidFace_prp ) then
      isRecvFace = tem_left

    ! notExist-remote+fromCoaser face (nothing to be received, this rank does
    ! not have any element shared by this face.)
    elseif( leftPrp.eq.tem_notExist_prp .and.                                  &
      &     rightPrp.eq.ior(tem_remoteFace_prp, tem_fromCoarserFace_prp)) then
      isRecvFace = 0

    ! remote+fromCoarser-remote face (nothing to be received, everything is
    ! remote)
    elseif( leftPrp.eq.ior(tem_remoteFace_prp, tem_fromCoarserFace_prp) .and.  &
      &     rightPrp.eq.tem_notExist_prp ) then
      isRecvFace = 0

    ! remote+fromCoarser-remote+fromCoaser face (nothing to be received,
    ! everything is remote)
    elseif( leftPrp.eq.ior(tem_remoteFace_prp, tem_fromCoarserFace_prp) .and.  &
      &     rightPrp.eq.ior(tem_remoteFace_prp, tem_fromCoarserFace_prp)) then
      isRecvFace = 0

    ! fromFiner - not exist face (nothing to be done)
    elseif( leftPrp.eq.tem_fromFinerFace_prp .and.                             &
      &     rightPrp.eq.tem_notExist_prp ) then
      isRecvFace = 0

    ! notExist - fromFiner face (nothing to be done)
    elseif( leftPrp.eq.tem_notExist_prp .and.                                  &
      &     rightPrp.eq.tem_fromFinerFace_prp ) then
      isRecvFace = 0

    ! fromFiner+remote - fluid face (we do not have to recv info here, but we
    ! send it)
    elseif( leftPrp.eq.ior(tem_fromFinerFace_prp,tem_remoteFace_prp) .and.     &
      &     rightPrp.eq.tem_fluidFace_prp ) then
      isRecvFace = 0

    ! fluid - remote+fromFiner face (we do not recv info here, but we send it)
    elseif( leftPrp.eq.tem_fluidFace_prp .and.                                 &
      &     rightPrp.eq.ior(tem_fromFinerFace_prp,tem_remoteFace_prp)) then
      isRecvFace = 0

    ! fromCoarser - notExist face (nothing to receive)
    elseif( leftPrp.eq.tem_fromCoarserFace_prp .and.                           &
      &     rightPrp.eq.tem_notExist_prp ) then
      isRecvFace = 0

    ! notExist - fromCoarser face (nothing to receive)
    elseif( leftPrp.eq.tem_notExist_prp .and.                                  &
      &     rightPrp.eq.tem_fromCoarserFace_prp ) then
      isRecvFace = 0

    ! remote+fromFiner-remote face (everything is remote)
    elseif( leftPrp.eq.ior(tem_fromFinerFace_prp, tem_remoteFace_prp) .and.    &
      &     rightPrp.eq.tem_remoteFace_prp) then
      isRecvFace = 0

    ! remote-remote+fromFiner face (everything is remote)
    elseif( leftPrp.eq. tem_remoteFace_prp .and.                               &
      &     rightPrp.eq.ior(tem_remoteFace_prp, tem_fromFinerFace_prp)) then
      isRecvFace = 0

    ! notExist-remote+fromFiner face (nothing is on my rank)
    elseif( leftPrp.eq. tem_notExist_prp .and.                                 &
      &     rightPrp.eq.ior(tem_remoteFace_prp, tem_fromFinerFace_prp)) then
      isRecvFace = 0

    ! remote+fromFiner-notExist face (nothing is on my rank)
    elseif( leftPrp.eq.ior(tem_remoteFace_prp, tem_fromFinerFace_prp) .and.    &
      &     rightPrp.eq.tem_notExist_prp ) then
      isRecvFace = 0

    ! fromFiner-fromCoarser face (everything is local, this is an intermediate
    ! face level (i.e. leveljumps of difference larger than 1 occure here) )
    elseif( leftPrp.eq.tem_fromFinerFace_prp .and.                             &
      &     rightPrp.eq.tem_fromCoarserFace_prp ) then
      isRecvFace = 0

    ! fromCoarser-fromFiner face (everything is local, this is an intermediate
    ! face level (i.e. leveljumps of difference larger than 1 occure here) )
    elseif( leftPrp.eq.tem_fromCoarserFace_prp .and.                           &
      &     rightPrp.eq.tem_fromFinerFace_prp ) then
      isRecvFace = 0

    ! fromFiner+remote - fromCoarser face (we do not have to recv info here,
    ! but we send it)
    elseif( leftPrp.eq.ior(tem_fromFinerFace_prp,tem_remoteFace_prp) .and.     &
      &     rightPrp.eq.tem_fromCoarserFace_prp ) then
      isRecvFace = 0

    ! fromCoarser - remote+fromFiner face (we do not recv info here, but we
    ! send it)
    elseif( leftPrp.eq.tem_fromCoarserFace_prp .and.                           &
      &     rightPrp.eq.ior(tem_fromFinerFace_prp,tem_remoteFace_prp)) then
      isRecvFace = 0

    ! remote+fromFiner - remote+fromFiner face (everything is remote, nothing
    ! to be done)
    elseif( leftPrp.eq.ior(tem_fromFinerFace_prp,tem_remoteFace_prp) .and.     &
      &     rightPrp.eq.ior(tem_fromFinerFace_prp,tem_remoteFace_prp)) then
      isRecvFace = 0

    ! remote+fromFiner - fromFiner face (everything is from finer, nothing to
    ! be done)
    elseif( leftPrp.eq.ior(tem_fromFinerFace_prp,tem_remoteFace_prp) .and.     &
      &     rightPrp.eq.tem_fromFinerFace_prp ) then
      isRecvFace = 0

    ! fromFiner - remote+fromFiner face (everything is from finer, nothing to
    ! be done)
    elseif( leftPrp.eq.tem_fromFinerFace_prp .and.                             &
      &     rightPrp.eq.ior(tem_fromFinerFace_prp,tem_remoteFace_prp)) then
      isRecvFace = 0

    ! fromFiner - remote face (nothing to be done)
    elseif( leftPrp.eq.tem_fromFinerFace_prp .and.                             &
      &     rightPrp.eq.tem_remoteFace_prp ) then
      isRecvFace = 0

    ! remote - fromFiner face (nothing to be done)
    elseif( leftPrp.eq.tem_remoteFace_prp .and.                                &
      &     rightPrp.eq.tem_fromFinerFace_prp ) then
      isRecvFace = 0

    ! remote+fromCoarser - fromCoarser face (elements are on coarser
    ! level, here and remote so nothing to be done)
    elseif( leftPrp.eq.ior(tem_fromCoarserFace_prp,tem_remoteFace_prp) .and.   &
      &     rightPrp.eq.tem_fromCoarserFace_prp ) then
      isRecvFace = 0

    ! fromCoarser - remote+fromCoarser face (elements are on coarser
    ! level, here and remote so nothing to be done)
    elseif( leftPrp.eq.tem_fromCoarserFace_prp  .and.                          &
      &     rightPrp.eq.ior(tem_fromCoarserFace_prp,tem_remoteFace_prp)) then
      isRecvFace = 0

    ! fluid - boundary face (nothing to be done)
    elseif(leftPrp.eq.tem_fluidFace_prp .and. rightPrp.eq.tem_bndFace_prp) then
      isRecvFace = 0

    ! boundary - fluid face (nothing to be done)
    elseif(leftPrp.eq.tem_bndFace_prp .and. rightPrp.eq.tem_fluidFace_prp) then
      isRecvFace = 0

    ! remote - boundary face (nothing to be done)
    elseif(leftPrp.eq.tem_remoteFace_prp .and. rightPrp.eq.tem_bndFace_prp) then
      isRecvFace = 0

    ! boundary - remote face (nothing to be done)
    elseif(leftPrp.eq.tem_bndFace_prp .and. rightPrp.eq.tem_remoteFace_prp) then
      isRecvFace = 0

    ! fromFiner - boundary face (nothing to be done)
    elseif( leftPrp.eq.tem_fromFinerFace_prp .and.                             &
      &     rightPrp.eq.tem_bndFace_prp ) then
      isRecvFace = 0

    ! boundary - fromFiner face (nothing to be done)
    elseif( leftPrp.eq.tem_bndFace_prp .and.                                   &
      &     rightPrp.eq.tem_fromFinerFace_prp ) then
      isRecvFace = 0

    ! fromFiner+remote - boundary face (nothing to be done)
    elseif( leftPrp.eq.ior(tem_fromFinerFace_prp,tem_remoteFace_prp) .and.     &
      &     rightPrp.eq.tem_bndFace_prp ) then
      isRecvFace = 0

    ! boundary - fromFiner+remote face (nothing to be done)
    elseif( leftPrp.eq.tem_bndFace_prp .and.                                   &
      &     rightPrp.eq.ior(tem_fromFinerFace_prp,tem_remoteFace_prp)) then
      isRecvFace = 0

    ! fromCoarser - boundary face (nothing to be done)
    elseif( leftPrp.eq.tem_fromCoarserFace_prp .and.                           &
      &     rightPrp.eq.tem_bndFace_prp ) then
      isRecvFace = 0

    ! boundary - fromCoarser face (nothing to be done)
    elseif( leftPrp.eq.tem_bndFace_prp .and.                                   &
      &     rightPrp.eq.tem_fromCoarserFace_prp ) then
      isRecvFace = 0

    end if

  end function tem_isRecvFace
! ****************************************************************************** !


  ! ************************************************************************** !
  !> summary: Subroutine to build communication buffers for the faces the
  !! current rank has to send information for (before the other ranks can make
  !! the compute step).
  subroutine tem_build_faceSendBuffers( levelDesc, faces, buf )
    ! --------------------------------------------------------------------------
    !> Dimension-by-dimension level descriptor for the current level and
    !! direction.
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> The communication pattern you want use for the buffer.
    ! type(tem_commpattern_type), intent(in) :: commPattern
    !> The created face descriptor.
    type(tem_face_descriptor_type),intent(inout)  :: faces
    !> The created send buffer. Size is two, due to left and right limes of
    !! the face values. To access it use tem_left and tem_right.
    type(tem_communication_type), intent(out) :: buf(2)
    ! --------------------------------------------------------------------------
    integer :: iFace, iProc, iSide, elemPos, faceSide, tIdPos, rank, rankPos,  &
      &        faceIndex, elemAddPos
    logical :: wasAdded
    integer(kind=long_k) :: elemId, trgElemId
    ! The element treeids I will send data for. We have to differentiate the
    ! left and right face of a face here.
    type(dyn_longArray_type) :: elemIds(2)
    ! The elements positions I will send data for. We have to differentiate the
    ! left and right face of a face here.
    type(grw_intArray_type) :: elemPositions(2)
    ! The source ranks.
    type(dyn_intArray_type) :: trgRank(2)
    ! The position of the source ranks.
    type(grw_intArray_type) :: trgRankPos(2)
    ! --------------------------------------------------------------------------

    ! Init our growing arrays.
    call init( me = elemIds(1) , length = 16 )
    call init( me = elemIds(2) , length = 16 )
    call init( me = elemPositions(1) , length = 16)
    call init( me = elemPositions(2) , length = 16)
    call init( me = trgRankPos(1) , length = 16)
    call init( me = trgRankPos(2) , length = 16)
    call init( me = trgRank(1) , length = 16 )
    call init( me = trgRank(2) , length = 16 )

    ! Iterate over all the faces and count the communicate faces.
    do iFace = 1, faces%faceList%faceId%nVals

      ! Check whether we will send info of one of the adjacent elements
      ! of this face.
      faceSide = tem_isSendFace(iFace, faces)
      if ( faceSide /= 0 ) then

        ! Check for which adjacent element (left or right) the current rank
        ! will send information about and count it.
        if (faceSide == tem_left) then
          ! The position and id of the elment I will send data for
          elemPos = faces%faceList%leftElemPos%val(iFace)
          elemId = faces%faceList%faceId%val(iFace)
          ! I will send information about my left element, so I will send
          ! the right face info of this element.
          faceSide = tem_right
          ! Get the element id in the opposite direction of the face. This
          ! will be used to determine the target rank.
          trgElemId = faces%faceList%rightElemId%val(iFace)
        else
          ! The position and id of the elment I will send data for
          elemPos = faces%faceList%rightElemPos%val(iFace)
          elemId = faces%faceList%rightElemId%val(iFace)
          ! I will send info about the right element, so I send info about
          ! it's left face.
          faceSide = tem_left
          ! Get the element id in the opposite direction of the face. This
          ! will be used to determine the target rank.
          trgElemId = faces%faceList%faceId%val(iFace)
        end if

        ! Get the source rank of the target element id. This will give us
        ! the target rank of the face.
        ! Therefore we have a look at the level
        ! descriptor of the current spatial direction and level.
        tIdPos = PositionOfVal( me = levelDesc%elem%tId, val = trgElemId )
        if (tIdPos <= 0) then
          write(*,*) 'ERROR in tem_build_faceSendBuffers: not able to '//  &
            &        'identify target rank for a certain face, stopping ...'
          call tem_abort()
        end if
        rank = levelDesc%elem%sourceProc%val(tIdPos)

        ! Store element position and soruce proc in a dynamic array.
        call append( me       = elemIds(faceSide), &
          &          val      = elemId,            &
          &          pos      = elemAddPos,        &
          &          wasAdded = wasAdded           )
        if(wasAdded) then
          call append( me = elemPositions(faceSide), val = elemPos )
          call append( me = trgRank(faceSide), val = rank, pos = rankPos )
          call append( me = trgRankPos(faceSide), val = rankPos )
        end if

      end if
    end do

    ! Now, we build the buffer for the faces (left and right values).
    do iSide = 1, 2
      ! Init the datastructures in the buffer
      buf(iSide)%nProcs = trgRank(iSide)%nVals
      allocate(buf(iSide)%proc( trgRank(iSide)%nVals ))
      allocate(buf(iSide)%nElemsProc( trgRank(iSide)%nVals ))
      allocate(buf(iSide)%elemPos( trgRank(iSide)%nVals ))
      do iProc = 1, buf(iSide)%nProcs
        call init( me = buf(iSide)%elemPos(iProc), length = 16 )
      end do

      ! Set the positions of the elements I will receive.
      do iFace = 1, elemPositions(iSide)%nVals
        faceIndex = elemIds(iSide)%sorted(iFace)
        elemPos = elemPositions(iSide)%val(faceIndex)
        rankPos = trgRankPos(iSide)%val(faceIndex)
        call append( me = buf(iSide)%elemPos(rankPos), val = elemPos )
      end do

      ! Count the number of elements we recv from each proc.
      do iProc = 1, buf(iSide)%nProcs
        buf(iSide)%proc(iProc) = trgRank(iSide)%val(iProc) - 1
        buf(iSide)%nElemsProc(iProc) = buf(iSide)%elemPos(iProc)%nVals
      end do
    end do

  end subroutine tem_build_faceSendBuffers
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> summary: Subroutine to check if the current rank will send information
  !! for the a given face before the compute step.
  !!
  pure function tem_isSendFace(facePos, faces) result( isSendFace )
    ! --------------------------------------------------------------------------
    !> The position of the face in faces to be checked.
    integer , intent(in) :: facePos
    !> The description of the faces.
    type(tem_face_descriptor_type),intent(in)  :: faces
    !> Returns 0 if no information has to be received from remote ranks
    !! before a compute step can be done. In case that information has to be
    !! received it will return a non-zero positive value. Furhtermore,
    !! it determines whether information for the left face value
    !! (left adjacent element)
    !! or the right face value (right adjacent element) has to be received.
    !! In the first case it will return
    !! [[tem_faceData_module:tem_left]] and in the second case
    !! it will return
    !! [[tem_faceData_module:tem_right]].
    !! If something went wrong, e.g. an unknown combination occurs,
    !! it will return -1.
    integer :: isSendFace
    ! --------------------------------------------------------------------------
    integer :: leftPrp, rightPrp
    ! --------------------------------------------------------------------------

    leftPrp = faces%faceList%leftPrp%val(facePos)
    rightPrp = faces%faceList%rightPrp%val(facePos)

    issendface = -1 ! default value (error)

    ! So, check for all conditions under which face info has to be send.

    ! Fluid-fluid face (purely local)
    ! -> nothing to send
    if( leftPrp.eq.tem_fluidFace_prp .and. rightPrp.eq.tem_fluidFace_prp ) then
      isSendFace = 0

    ! Fluid-remote face (the halo is not refined)
    ! -> left element is on my rank, so I will compute it. Therefore, I will
    !    receive information about the right element, but I do not have to send
    !    anything.
    elseif( leftPrp.eq.tem_fluidFace_prp .and.                                 &
      &     rightPrp.eq.tem_remoteFace_prp ) then
      isSendFace = 0

    ! Remote-fluid face (the halo is not refined)
    ! -> left element is on a remote rank, so I will not compute it. Therfore
    !    I will send the element to the remote rank.
    elseif( leftPrp.eq.tem_remoteFace_prp .and.                                &
      &     rightPrp.eq.tem_fluidFace_prp ) then
      isSendFace = tem_right

    ! Remote-notExisting face (the remote is not refined)
    ! -> remote is not located on my rank and the other element does not
    !    exist. This face does not have to be computed on any rank, so my rank
    !    is not sending any info about any element adjacent to this face.
    elseif( leftPrp.eq.tem_remoteFace_prp .and.                                &
      &     rightPrp.eq.tem_notExist_prp ) then
      isSendFace = 0

    ! notExisting-remote face (the remote is not refined)
    ! -> remote is not located on my rank and the other element does not
    !    exist. This face does not have to be computed on any rank, so my rank
    !    is not send any info about any element adjacent to this face.
    elseif( leftPrp.eq.tem_notExist_prp .and.                                  &
      &     rightPrp.eq.tem_remoteFace_prp ) then
      isSendFace = 0

    ! Fluid-fromFiner face (purely local)
    ! -> everything is local, so there is no need to send or recv information
    !    from another rank.
    elseif( leftPrp.eq.tem_fluidFace_prp .and.                                 &
      &     rightPrp.eq.tem_fromFinerFace_prp ) then
      isSendFace = 0

    ! FromFiner-fluid face (purely local)
    ! -> everything is local, so there is no need to send or recv information
    !    from another rank.
    elseif( leftPrp.eq.tem_fromFinerFace_prp .and.                             &
      &     rightPrp.eq.tem_fluidFace_prp ) then
      isSendFace = 0

    ! FromFiner-FromFiner face (purely local)
    ! -> everything is local. Somehow two ghost elements meet each other, but
    !    we do not have to recv or send any data for these faces.
    elseif( leftPrp.eq.tem_fromFinerFace_prp .and.                             &
      &     rightPrp.eq.tem_fromFinerFace_prp ) then
      isSendFace = 0

    ! FromCoarser-fluid face (purely local)
    ! -> everything is local, so no data has to be send or received.
    elseif( leftPrp.eq.tem_fromCoarserFace_prp .and.                           &
      &     rightPrp.eq.tem_fluidFace_prp ) then
      isSendFace = 0

    ! Fluid-fromCoarser face (purely local)
    ! -> everything is local, so no data has to be send or received.
    elseif( leftPrp.eq.tem_fluidFace_prp .and.                                 &
      &     rightPrp.eq.tem_fromCoarserFace_prp ) then
      isSendFace = 0

    ! FromCoarser-fromCoarser face (purely local)
    ! -> everything is local, so no data has to be send or received.
    elseif( leftPrp.eq.tem_fromCoarserFace_prp .and.                           &
      &     rightPrp.eq.tem_fromCoarserFace_prp ) then
      isSendFace = 0

    ! Remote-fromCoarser face
    ! -> The remote fluid elements are located on a finer level
    elseif( leftPrp.eq.tem_remoteFace_prp .and.                                &
      &     rightPrp.eq.tem_fromCoarserFace_prp ) then
      isSendFace = tem_right

    ! FromCoarser-remote face
    ! -> The remote fluid elements are located on a finer level
    elseif( leftPrp.eq.tem_fromCoarserFace_prp .and.                           &
      &     rightPrp.eq.tem_remoteFace_prp ) then
      isSendFace = tem_left

    ! Remote-remote face
    ! -> The halo elements meet each other. This can happen, but we do not have
    !    to send or receive information for this face, because we do not
    !    compute anything here.
    elseif( leftPrp.eq.tem_remoteFace_prp .and.                                &
      &     rightPrp.eq.tem_remoteFace_prp ) then
      isSendFace = 0

    ! notExist - remote face (nothing to be send, we do not have info about
    ! adjacent elements)
    elseif( leftPrp.eq.tem_notExist_prp .and.                                  &
      &     rightPrp.eq.tem_remoteFace_prp ) then
      isSendFace = 0

    ! remote - notExist face (nothing to be send, we do not have any info about
    ! adjacent elements)
    elseif( leftPrp.eq.tem_remoteFace_prp .and.                                &
      &     rightPrp.eq.tem_notExist_prp ) then
      isSendFace = 0

    ! fluid-remote+fromCoaser face (here, we receive data in some situations,
    ! but we never send info in this case)
    elseif( leftPrp.eq.tem_fluidFace_prp .and.                                 &
      &     rightPrp.eq.ior(tem_remoteFace_prp, tem_fromCoarserFace_prp)) then
      isSendFace = 0

    ! remote+fromCoarser-fluid face (here, we receive data in some situations,
    ! but we never send info for this faces)
    elseif( leftPrp.eq.ior(tem_remoteFace_prp,tem_fromCoarserFace_prp) .and.   &
      &     rightPrp.eq.tem_fluidFace_prp ) then
      isSendFace = 0

    ! notExist-remote+fromCoaser face (nothing to be send here)
    elseif( leftPrp.eq.tem_notExist_prp .and.                                  &
      &     rightPrp.eq.ior(tem_remoteFace_prp, tem_fromCoarserFace_prp)) then
      isSendFace = 0

    ! remote+fromCoarser-remote face (everything is remote, so nothing to be
    ! sent)
    elseif( leftPrp.eq.ior(tem_remoteFace_prp, tem_fromCoarserFace_prp) .and.  &
      &     rightPrp.eq.tem_notExist_prp ) then
      isSendFace = 0

    ! remote+fromCoarser-remote+fromCoaser face (everything is remote, so
    ! nothing to be sent)
    elseif( leftPrp.eq.ior(tem_remoteFace_prp, tem_fromCoarserFace_prp) .and.  &
      &     rightPrp.eq.ior(tem_remoteFace_prp, tem_fromCoarserFace_prp)) then
      isSendFace = 0

    ! fromFiner - not exist face (nothing to be send here).
    elseif( leftPrp.eq.tem_fromFinerFace_prp .and.                             &
      &     rightPrp.eq.tem_notExist_prp ) then
      isSendFace = 0

    ! notExist - fromFiner face (nothing to be send here)
    elseif( leftPrp.eq.tem_notExist_prp .and.                                  &
      &     rightPrp.eq.tem_fromFinerFace_prp ) then
      isSendFace = 0

    ! fromFiner+remote - fluid face (will be send on the finer level)
    elseif( leftPrp.eq.ior(tem_fromFinerFace_prp,tem_remoteFace_prp) .and.     &
      &     rightPrp.eq.tem_fluidFace_prp ) then
      isSendFace = 0

    ! fluid - remote+fromFiner face (will be send on the finer level)
    elseif( leftPrp.eq.tem_fluidFace_prp .and.                                 &
      &     rightPrp.eq.ior(tem_fromFinerFace_prp,tem_remoteFace_prp)) then
      isSendFace = 0

    ! fromCoarser - notExist face (nothing to send)
    elseif( leftPrp.eq.tem_fromCoarserFace_prp .and.                           &
      &     rightPrp.eq.tem_notExist_prp ) then
      isSendFace = 0

    ! notExist - fromCoarser face (nothing to send)
    elseif( leftPrp.eq.tem_notExist_prp .and.                                  &
      &     rightPrp.eq.tem_fromCoarserFace_prp ) then
      isSendFace = 0

    ! remote+fromFiner-remote face (everything is remote)
    elseif( leftPrp.eq.ior(tem_fromFinerFace_prp, tem_remoteFace_prp) .and.    &
      &     rightPrp.eq.tem_remoteFace_prp ) then
      isSendFace = 0

    ! remote-remote+fromFiner face (everything is remote)
    elseif( leftPrp.eq. tem_remoteFace_prp .and.                               &
      &     rightPrp.eq.ior(tem_remoteFace_prp, tem_fromFinerFace_prp)) then
      isSendFace = 0

    ! notExist-remote+fromFiner face (nothing is on my rank)
    elseif( leftPrp.eq. tem_notExist_prp .and.                                 &
      &     rightPrp.eq.ior(tem_remoteFace_prp, tem_fromFinerFace_prp)) then
      isSendFace = 0

    ! remote+fromFiner-notExist face (nothing is on my rank)
    elseif( leftPrp.eq.ior(tem_remoteFace_prp, tem_fromFinerFace_prp) .and.    &
      &     rightPrp.eq.tem_notExist_prp ) then
      isSendFace = 0

    ! fromFiner-fromCoarser face (everything is local, this is an intermediate
    ! face level (i.e. leveljumps of difference larger than 1 occure here) )
    elseif( leftPrp.eq.tem_fromFinerFace_prp .and.                             &
      &     rightPrp.eq.tem_fromCoarserFace_prp ) then
      isSendFace = 0

    ! fromCoarser-fromFiner face (everything is local, this is an intermediate
    ! face level (i.e. leveljumps of difference larger than 1 occure here) )
    elseif( leftPrp.eq.tem_fromCoarserFace_prp .and.                           &
      &     rightPrp.eq.tem_fromFinerFace_prp ) then
      isSendFace = 0

    ! fromFiner+remote - fromCoarser face (send on a finer level)
    elseif( leftPrp.eq.ior(tem_fromFinerFace_prp,tem_remoteFace_prp) .and.     &
      &     rightPrp.eq.tem_fromCoarserFace_prp ) then
      isSendFace = 0

    ! fromCoarser - remote+fromFiner face (send on a finer level)
    elseif( leftPrp.eq.tem_fromCoarserFace_prp .and.                           &
      &     rightPrp.eq.ior(tem_fromFinerFace_prp,tem_remoteFace_prp)) then
      isSendFace = 0

    ! remote+fromFiner - remote+fromFiner face (everything is remote, nothing
    ! to be done)
    elseif( leftPrp.eq.ior(tem_fromFinerFace_prp,tem_remoteFace_prp) .and.     &
      &     rightPrp.eq.ior(tem_fromFinerFace_prp,tem_remoteFace_prp)) then
      isSendFace = 0

    ! remote+fromFiner - fromFiner face (everything is from finer, nothing to
    ! be done)
    elseif( leftPrp.eq.ior(tem_fromFinerFace_prp,tem_remoteFace_prp) .and.     &
      &     rightPrp.eq.tem_fromFinerFace_prp ) then
      isSendFace = 0

    ! fromFiner - remote+fromFiner face (everything is from finer, nothing to
    ! be done)
    elseif( leftPrp.eq.tem_fromFinerFace_prp .and.                             &
      &     rightPrp.eq.ior(tem_fromFinerFace_prp,tem_remoteFace_prp)) then
      isSendFace = 0

    ! fromFiner - remote face (nothing to be done)
    elseif( leftPrp.eq.tem_fromFinerFace_prp .and.                             &
      &     rightPrp.eq.tem_remoteFace_prp ) then
      isSendFace = 0

    ! remote - fromFiner face (nothing to be done)
    elseif( leftPrp.eq.tem_remoteFace_prp .and.                                &
      &     rightPrp.eq.tem_fromFinerFace_prp ) then
      isSendFace = 0

    ! remote+fromCoarser - fromCoarser face (elements are on coarser
    ! level, here and remote so nothing to be done)
    elseif( leftPrp.eq.ior(tem_fromCoarserFace_prp,tem_remoteFace_prp) .and.   &
      &     rightPrp.eq.tem_fromCoarserFace_prp ) then
      isSendFace = 0

    ! fromCoarser - remote+fromCoarser face (elements are on coarser
    ! level, here and remote so nothing to be done)
    elseif( leftPrp.eq.tem_fromCoarserFace_prp .and.                           &
      &     rightPrp.eq.ior(tem_fromCoarserFace_prp,tem_remoteFace_prp)) then
      isSendFace = 0

    ! fluid - boundary face (nothing to be done)
    elseif(leftPrp.eq.tem_fluidFace_prp .and. rightPrp.eq.tem_bndFace_prp) then
      isSendFace = 0

    ! boundary - fluid face (nothing to be done)
    elseif(leftPrp.eq.tem_bndFace_prp .and. rightPrp.eq.tem_fluidFace_prp) then
      isSendFace = 0

    ! remote - boundary face (nothing to be done)
    elseif(leftPrp.eq.tem_remoteFace_prp .and. rightPrp.eq.tem_bndFace_prp) then
      isSendFace = 0

    ! boundary - remote face (nothing to be done)
    elseif(leftPrp.eq.tem_bndFace_prp .and. rightPrp.eq.tem_remoteFace_prp) then
      isSendFace = 0

    ! fromFiner - boundary face (nothing to be done)
    elseif( leftPrp.eq.tem_fromFinerFace_prp .and.                             &
      &     rightPrp.eq.tem_bndFace_prp ) then
      isSendFace = 0

    ! boundary - fromFiner face (nothing to be done)
    elseif( leftPrp.eq.tem_bndFace_prp .and.                                   &
      &     rightPrp.eq.tem_fromFinerFace_prp ) then
      isSendFace = 0

    ! fromFiner+remote - boundary face (nothing to be done)
    elseif( leftPrp.eq.ior(tem_fromFinerFace_prp,tem_remoteFace_prp) .and.     &
      &     rightPrp.eq.tem_bndFace_prp ) then
      isSendFace = 0

    ! boundary - fromFiner+remote face (nothing to be done)
    elseif( leftPrp.eq.tem_bndFace_prp .and.                                   &
      &     rightPrp.eq.ior(tem_fromFinerFace_prp,tem_remoteFace_prp)) then
      isSendFace = 0

    ! fromCoarser - boundary face (nothing to be done)
    elseif( leftPrp.eq.tem_fromCoarserFace_prp .and.                           &
      &     rightPrp.eq.tem_bndFace_prp ) then
      isSendFace = 0

    ! boundary - fromCoarser face (nothing to be done)
    elseif( leftPrp.eq.tem_bndFace_prp .and.                                   &
      &     rightPrp.eq.tem_fromCoarserFace_prp ) then
      isSendFace = 0

    end if

  end function tem_isSendFace
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Creates dimension by dimension level descriptors.
  !!
  subroutine tem_dimByDim_construction( tree, boundary, commPattern, proc, &
    &                                   levelDescX, levelDescY, levelDescZ )
    ! --------------------------------------------------------------------------
    !> Tree representation of the mesh.
    type(treelmesh_type), intent(inout) :: tree
    !> The boundaries of your simulation domain
    type(tem_bc_prop_type), intent(in) :: boundary
    !> The communication pattern you use for the buffer.
    type(tem_commpattern_type), intent(in) :: commPattern
    !> Level descriptor for each spatial direction and each level of your mesh.
    !! The level descriptor have to be constructed with the dimension by
    !! dimension stencils (+1, 0, -1) for each spatial direction.
    type(tem_levelDesc_type), allocatable, intent(out) :: levelDescX(:)
    !> Level descriptor for each spatial direction and each level of your mesh.
    !! The level descriptor have to be constructed with the dimension by
    !! dimension stencils (+1, 0, -1) for each spatial direction.
    type(tem_levelDesc_type), allocatable, intent(out) :: levelDescY(:)
    !> Level descriptor for each spatial direction and each level of your mesh.
    !! The level descriptor have to be constructed with the dimension by
    !! dimension stencils (+1, 0, -1) for each spatial direction.
    type(tem_levelDesc_type), allocatable, intent(out) :: levelDescZ(:)
    !> Process description to use.
    type(tem_comm_env_type), intent(in) :: proc
    ! --------------------------------------------------------------------------
    ! The dimension-by-dimension stencil list.
    type(tem_stencilHeader_type) :: stencilList(3)
    ! --------------------------------------------------------------------------

    ! Create the dimension by dimension stencils
    write(logUnit(3),*) 'Creating dimension by dimension stencils for face ' &
      &                 //  'descriptor...'
    call tem_define_dimStencil( stencilList(1), tree%nElems, 1 )
    call tem_define_dimStencil( stencilList(2), tree%nElems, 2 )
    call tem_define_dimStencil( stencilList(3), tree%nElems, 3 )

    ! Dump the stencils we used.
    call tem_stencil_dump( stencilList( 1 ))
    call tem_stencil_dump( stencilList( 2 ))
    call tem_stencil_dump( stencilList( 3 ))

    ! Now, build the level descriptors for each of the stencils.
    call tem_create_levelDesc( tree, stencilList(1), boundary, commPattern, &
      &                        levelDescX, proc                             )
    call tem_create_levelDesc( tree, stencilList(2), boundary, commPattern, &
      &                        levelDescY, proc                             )
    call tem_create_levelDesc( tree, stencilList(3), boundary, commPattern, &
      &                        levelDescZ, proc                             )

  end subroutine tem_dimByDim_construction
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Create a face level descriptor
  !!
  subroutine tem_create_levelDesc( tree, stencil, boundary, commPattern, &
    &                              levelDesc, proc                       )
    ! --------------------------------------------------------------------------
    !> Tree representation of the mesh.
    type(treelmesh_type), intent(inout) :: tree
    !> The stencil to create the level descriptor for.
    type(tem_stencilHeader_type), intent(inout) :: stencil(1)
    !> The boundaries of your simulation domain
    type(tem_bc_prop_type), intent(in) :: boundary
    !> The communication pattern you use for the buffer.
    type(tem_commpattern_type), intent(in) :: commPattern
    !> The created level descriptor.
    type(tem_levelDesc_type), allocatable, intent(out) :: levelDesc(:)
    !> Process description to use.
    type(tem_comm_env_type), intent(in) :: proc
    ! --------------------------------------------------------------------------
    integer, allocatable :: levelPointer(:)
    ! --------------------------------------------------------------------------

    call tem_init_elemLevels( me             = levelDesc, &
      &                       boundary       = boundary,  &
      &                       tree           = tree,      &
      &                       stencils       = stencil    )

    call tem_find_allElements( tree           = tree,         &
      &                        levelDesc      = levelDesc,    &
      &                        levelPointer   = levelPointer, &
      &                        computeStencil = stencil,      &
      &                        commpattern    = commPattern,  &
      &                        proc           = proc          )

    ! JQ: allocated in tem_init_elemLevels
    ! do iLevel=tree%global%minlevel,tree%global%maxlevel
    !   allocate(levelDesc(iLevel)%neigh(1))
    ! end do
    call tem_build_HorizontalDependencies( iStencil       = 1,         &
      &                                    levelDesc      = levelDesc, &
      &                                    tree           = tree,      &
      &                                    computeStencil = stencil(1) )

    call tem_build_VerticalDependencies( levelDesc = levelDesc,            &
      &                                  minLevel  = tree%global%minLevel, &
      &                                  maxLevel  = tree%global%maxLevel  )

    call tem_cleanupDependencyArrays( levelDesc = levelDesc )

  end subroutine tem_create_levelDesc
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> summary: Creates the dimension-by-dimension stencils for a given spatial
  !! direction.
  !!
  subroutine tem_define_dimStencil( stencil, nElems, direction)
    ! --------------------------------------------------------------------------
    !> The stencil layout to set.
    type(tem_stencilHeader_type), intent(inout) :: stencil
    !> The number of elements which share this stencil
    integer, intent(in) :: nElems
    !> The spatial direction:
    !! 1 -> x direction
    !! 2 -> y direction
    !! 3 -> z direction
    integer, intent(in) :: direction
    ! --------------------------------------------------------------------------
    integer :: iElem, offsetIndex
    ! --------------------------------------------------------------------------

    stencil%useAll = .true.
    stencil%nElems = nElems
    call init( me = stencil%elem, length = 0 )
    stencil%QQN = 2
    stencil%QQ = 2
    allocate( stencil%cxDir(3, 2) )
    allocate( stencil%cxDirInv(2) )

    do iElem = 1,2
      offsetIndex = (iElem-1)*3 + direction
      stencil%cxDir(1:3, iElem) = qOffset(offsetIndex, 1:3)
      stencil%cxDirInv(iElem) = qInvDir(offsetIndex)
    end do

  end subroutine tem_define_dimStencil
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> summary: Routine to generate the separated lists of faces for
  !! communicated, refined/coarsened and local faces.
  !!
  subroutine tem_build_faceLists( minLevel, maxLevel, nEligibleChildren, faces )
    ! --------------------------------------------------------------------------
    !> Minimum level of your mesh.
    integer, intent(in) :: minLevel
    !> Maximum level of your mesh.
    integer, intent(in) :: maxLevel
    !> The number of eligible children for the vertical face dependency
    integer, intent(in) :: nEligibleChildren
    !> The created face descriptor.
    type(tem_face_type),intent(inout) :: faces(minLevel:maxLevel)
    ! --------------------------------------------------------------------------
    integer :: iLevel
    ! --------------------------------------------------------------------------

    levelLoop: do iLevel = minLevel, maxLevel

      ! Build the list of compute faces for the current level
      call tem_build_computeList( faces(iLevel), nEligibleChildren )

    end do levelLoop

    ! Build the list of from finer faces for the current level
    call tem_build_fromFinerList( minLevel, maxLevel, nEligibleChildren, faces )

  end subroutine tem_build_faceLists
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Extracts all the faces which will be computed by current rank.
  !!
  subroutine tem_build_fromFinerList( minLevel, maxLevel, nEligibleChildren, &
                                    & faces )
    ! --------------------------------------------------------------------------
    !> The min refinement level of the mesh.
    integer, intent(in) :: minLevel
    !> The max refinement level of the mesh.
    integer, intent(in) :: maxLevel
    !> The number of eligible children for the vertical face dependency
    integer, intent(in) :: nEligibleChildren
    !> The created face descriptor.
    type(tem_face_type),intent(inout) :: faces(minLevel:maxLevel)
    ! --------------------------------------------------------------------------
    integer :: nFromFinerFaces(3,2)
    integer :: iDir, iFace, iFromFinerFace(2), iDep, iLevel
    integer :: childFace, childFaceOp, childPos, childPosOp, isFromFiner
    ! --------------------------------------------------------------------------

    ! Loop over all the levels
    levelLoop: do iLevel = minLevel, maxLevel

      ! First, we count all the element which will be covered by this rank.
      nFromFinerFaces = 0
      do iDir = 1,3
        do iFace = 1, faces(iLevel)%faces(iDir)%faceList%faceId%nVals

          ! decide if the element will be computed on this level and this rank.
          isFromFiner = tem_isFromFinerFace(iFace, faces(iLevel)%faces(iDir))
          if( isFromFiner .eq. tem_left) then
            ! If left element is refined, its right face has to be interpolated.
            nFromFinerFaces(iDir,tem_right) =                          &
              &                     nFromFinerFaces(iDir,tem_right) + 1
          elseif( isFromFiner .eq. tem_right) then
            ! If right element is refined, its left face has to be interpolated.
            nFromFinerFaces(iDir,tem_left) =                           &
              &                      nFromFinerFaces(iDir,tem_left) + 1
          end if

        end do
      end do

      ! Allocate memory for the from finer faces
      do iDir = 1,3
        do iFace = 1,2
          allocate( faces(iLevel)%faces(iDir)%fromFinerFace(iFace)%            &
            &                          elemPos( nFromFinerFaces(iDir, iFace) ))
          allocate( faces(iLevel)%faces(iDir)%fromFinerFace(iFace)%            &
            &      childPos( nEligibleChildren, nFromFinerFaces(iDir, iFace) ))
          allocate( faces(iLevel)%faces(iDir)%fromFinerFace(iFace)%            &
            &                        elemPosOp( nFromFinerFaces(iDir, iFace) ))
          allocate( faces(iLevel)%faces(iDir)%fromFinerFace(iFace)%            &
            &    childPosOp( nEligibleChildren, nFromFinerFaces(iDir, iFace) ))
          allocate( faces(iLevel)%faces(iDir)%fromFinerFace(iFace)%            &
            &                          facePos( nFromFinerFaces(iDir, iFace) ))
        end do
      end do

      ! Now, we add all the compute elements to the iterator
      do iDir = 1,3
        iFromFinerFace(:) = 1
        do iFace = 1, faces(iLevel)%faces(iDir)%faceList%faceId%nVals

          ! Decide if the face has to be interpolated.
          isFromFiner = tem_isFromFinerFace(iFace, faces(iLevel)%faces(iDir))
          if( isFromFiner .eq. tem_left) then
            ! The left element is refined -> add the left element to the right
            ! faces that have to be interpolated.
            faces(iLevel)%faces(iDir)%fromFinerFace(tem_right)%        &
              &       elemPos(iFromFinerFace(tem_right)) =             &
              &       faces(iLevel)%faces(iDir)%faceList%leftElemPos%val(iFace)
            faces(iLevel)%faces(iDir)%fromFinerFace(tem_right)%        &
              &      elemPosOp(iFromFinerFace(tem_right)) =            &
              &      faces(iLevel)%faces(iDir)%faceList%rightElemPos%val(iFace)
            faces(iLevel)%faces(iDir)%fromFinerFace(tem_right)%        &
              &              facePos(iFromFinerFace(tem_right)) = iFace
            ! Now, we set the dependencies, i.e. we identify the positions of
            ! the child elements and set them in the from finer description.
            do iDep = 1,nEligibleChildren
              childFace = faces(iLevel)%faces(iDir)%faceDep%                   &
                &                                      childFacePos(iDep,iFace)
              childPos = faces(iLevel+1)%faces(iDir)%faceList%                 &
                &                                    leftElemPos%val(childFace)
              faces(iLevel)%faces(iDir)%fromFinerFace(tem_right)%      &
                &   childPos(iDep,iFromFinerFace(tem_right)) = childPos
              childFaceOp = faces(iLevel)%faces(iDir)%faceDep%                 &
                &                                      childFacePosOp(iDep,iFace)
              childPosOp = faces(iLevel+1)%faces(iDir)%faceList%               &
                &                                   rightElemPos%val(childFaceOp)
              faces(iLevel)%faces(iDir)%fromFinerFace(tem_right)%      &
                &            childPosOp(iDep,iFromFinerFace(tem_right))&
                & = childPosOp
            end do

            ! Count the number of right faces
            iFromFinerFace(tem_right) =                                &
              &                           iFromFinerFace(tem_right) + 1

          elseif( isFromFiner .eq. tem_right) then
            ! The right element is refined -> add the right element to the left
            ! faces that have to be interpolated.
            faces(iLevel)%faces(iDir)%fromFinerFace(tem_left)%         &
              &      elemPos(iFromFinerFace(tem_left)) =               &
              &      faces(iLevel)%faces(iDir)%faceList%rightElemPos%val(iFace)
            faces(iLevel)%faces(iDir)%fromFinerFace(tem_left)%         &
              &       elemPosOp(iFromFinerFace(tem_left)) =            &
              &       faces(iLevel)%faces(iDir)%faceList%leftElemPos%val(iFace)
            faces(iLevel)%faces(iDir)%fromFinerFace(tem_left)%         &
              &               facePos(iFromFinerFace(tem_left)) = iFace
            ! Now, we set the dependencies.
            ! Now, we set the dependencies, i.e. we identify the positions of
            ! the child elements and set them in the from finer description.
            do iDep = 1,nEligibleChildren
              childFace = faces(iLevel)%faces(iDir)%faceDep%                   &
                &                                      childFacePos(iDep,iFace)
              childPos = faces(iLevel+1)%faces(iDir)%faceList%                 &
                &                                   rightElemPos%val(childFace)
              faces(iLevel)%faces(iDir)%fromFinerFace(tem_left)%       &
                &    childPos(iDep,iFromFinerFace(tem_left)) = childPos
              childFaceOp = faces(iLevel)%faces(iDir)%faceDep%                 &
                &                                      childFacePosOp(iDep,iFace)
              childPosOp = faces(iLevel+1)%faces(iDir)%faceList%               &
                &                                    leftElemPos%val(childFaceOp)
              faces(iLevel)%faces(iDir)%fromFinerFace(tem_left)%       &
                & childPosOp(iDep,iFromFinerFace(tem_left)) = childPosOp
            end do

            ! Count the number of left faces
            iFromFinerFace(tem_left) =                                 &
              &                            iFromFinerFace(tem_left) + 1

          end if
        end do
      end do
    end do levelLoop

  end subroutine tem_build_fromFinerList
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> summary: Function to decide if a certain face is constructed from finer
  !! faces on the current rank.
  !!
  function tem_isFromFinerFace(facePos, faces) result( isFromFiner )
    ! --------------------------------------------------------------------------
    integer, intent(in) :: facePos
    type(tem_face_descriptor_type), intent(in) :: faces
    integer :: isFromFiner
    ! --------------------------------------------------------------------------
    integer :: leftPrp, rightPrp
    ! --------------------------------------------------------------------------

    leftPrp = faces%faceList%leftPrp%val(facePos)
    rightPrp = faces%faceList%rightPrp%val(facePos)

    ! Fluid-fluid face (purely local)
    if( leftPrp.eq.tem_fluidFace_prp .and. rightPrp.eq.tem_fluidFace_prp ) then
      isFromFiner = 0

    ! Fluid-ghostFromFiner face (purley local, no remote property set)
    elseif( leftPrp.eq.tem_fluidFace_prp .and.                                 &
      &     rightPrp.eq.tem_fromFinerFace_prp ) then
      isFromFiner = tem_right

    ! GhostFromFiner-fluid face (purley local, no remote property set)
    elseif( leftPrp.eq.tem_fromFinerFace_prp .and.                             &
      &     rightPrp.eq.tem_fluidFace_prp ) then
      isFromFiner = tem_left

    ! GhostFromFiner-GhostFromFiner face (purley local, no remote property set)
    ! --> two ghost elements from finer meet each other somehow. However, we
    !     do not have to interpolate for these face combination.
    elseif( leftPrp.eq.tem_fromFinerFace_prp .and.                             &
      &     rightPrp.eq.tem_fromFinerFace_prp ) then
      isFromFiner = 0

    ! Fluid-ghostFromCoarser face (purely local, no remote property)
    ! --> no from finer face.
    elseif( leftPrp.eq.tem_fluidFace_prp .and.                                 &
      &     rightPrp.eq.tem_fromCoarserFace_prp ) then
      isFromFiner = 0

    ! GhostFromCoarser-fluid face (purely local, no remote property)
    ! --> no from finer face.
    elseif( leftPrp.eq.tem_fromCoarserFace_prp .and.                           &
      &     rightPrp.eq.tem_fluidFace_prp ) then
      isFromFiner = 0

    ! GhostFromCoarser-GhostFromCoarser face (purely local, no remote property)
    ! --> no from finer face. Somehow, two ghost from coarser meet each other,
    !     but it is not from finer face.
    elseif( leftPrp.eq.tem_fromCoarserFace_prp .and.                           &
      &     rightPrp.eq.tem_fromCoarserFace_prp ) then
      isFromFiner = 0

    ! Fluid-remote face (no refinement property set)
    ! --> no refinement, so it is not from finer.
    elseif( leftPrp.eq.tem_fluidFace_prp .and.                                 &
      &     rightPrp.eq.tem_remoteFace_prp ) then
      isFromFiner = 0

    ! Remote-fluid face (no refinement property set)
    ! --> no refinement, so it is not from finer.
    elseif( leftPrp.eq.tem_remoteFace_prp .and.                                &
      &     rightPrp.eq.tem_fluidFace_prp ) then
      isFromFiner = 0

    ! GhostFromCoarser-remote face
    elseif( leftPrp.eq.tem_fromCoarserFace_prp .and.                           &
      &     rightPrp.eq.tem_remoteFace_prp ) then
      isFromFiner = 0

    ! Remote-GhostFromCoarser face
    elseif( leftPrp.eq.tem_remoteFace_prp .and.                                &
      &     rightPrp.eq.tem_fromCoarserFace_prp ) then
      isFromFiner = 0

    ! Remote-remote face (somehow, two halo elements are meeting)
    elseif( leftPrp.eq.tem_remoteFace_prp .and.                                &
      &     rightPrp.eq.tem_remoteFace_prp ) then
      isFromFiner = 0

    ! notExist - remote face (nothing from a finer level)
    elseif( leftPrp.eq.tem_notExist_prp .and.                                  &
      &     rightPrp.eq.tem_remoteFace_prp ) then
      isFromFiner = 0

    ! remote - notExist face (nothing from a finer level)
    elseif( leftPrp.eq.tem_remoteFace_prp .and.                                &
      &     rightPrp.eq.tem_notExist_prp ) then
      isFromFiner = 0

    ! fluid-remote+fromCoaser face (nothing is from a finer level)
    elseif( leftPrp.eq.tem_fluidFace_prp .and.                                 &
      &     rightPrp.eq.ior(tem_remoteFace_prp, tem_fromCoarserFace_prp)) then
      isFromFiner = 0

    ! remote+fromCoarser-fluid face (nothing is from a finer level)
    elseif( leftPrp.eq.ior(tem_remoteFace_prp,tem_fromCoarserFace_prp) .and.   &
      &     rightPrp.eq.tem_fluidFace_prp) then
      isFromFiner = 0

    ! notExist-remote+fromCoaser face (neither my ranks hold an element for this
    ! face, nor elements are from finer. So, this face is not from finer).
    elseif( leftPrp.eq.tem_notExist_prp .and.                                  &
      &     rightPrp.eq.ior(tem_remoteFace_prp, tem_fromCoarserFace_prp)) then
      isFromFiner = 0

    ! remote+fromCoarser-remote face (neither my ranks hold an element for this
    ! face, nor elements are from finer. So, this face is not from finer).
    elseif( leftPrp.eq.ior(tem_remoteFace_prp, tem_fromCoarserFace_prp) .and.  &
      &     rightPrp.eq.tem_notExist_prp ) then
      isFromFiner = 0

    ! remote+fromCoarser-remote+fromCoaser face (neither my ranks hold an
    ! element for this face, nor elements are from finer. So, this face is not
    ! from finer).
    elseif( leftPrp.eq.ior(tem_remoteFace_prp, tem_fromCoarserFace_prp) .and.  &
      &     rightPrp.eq.ior(tem_remoteFace_prp, tem_fromCoarserFace_prp)) then
      isFromFiner = 0

    ! fromFiner - not exist face (this situation happens if a coarser element is
    ! located next to finer elements which are completely located on my rank. So
    ! my rank has to interpolate them upwards before we send them to the coarse
    ! element on the remote rank)
    elseif( leftPrp.eq.tem_fromFinerFace_prp .and.                             &
      &     rightPrp.eq.tem_notExist_prp ) then
      isFromFiner = 0

    ! notExist - fromFiner face (this situation happens if a coarser element is
    ! located next to finer elements which are completely located on my rank. So
    ! my rank has to interpolate them upwards before we send them to the coarse
    ! element on the remote rank)
    elseif( leftPrp.eq.tem_notExist_prp .and.                                  &
      &     rightPrp.eq.tem_fromFinerFace_prp ) then
      isFromFiner = 0

    ! fromCoarser - notExist face (nothing is from finer)
    elseif( leftPrp.eq.tem_fromCoarserFace_prp .and.                           &
      &     rightPrp.eq.tem_notExist_prp ) then
      isFromFiner = 0

    ! notExist - fromCoarser face (nothing is from finer)
    elseif( leftPrp.eq.tem_notExist_prp .and.                                  &
      &     rightPrp.eq.tem_fromCoarserFace_prp ) then
      isFromFiner = 0

    ! fluid-remote+fromFiner face (data is send to the remote rank on this
    ! level. The remote rank will also do the interpolation for me, so my rank
    ! does not have to do anyhting here.)
    elseif( leftPrp.eq.tem_fluidFace_prp .and.                                 &
      &     rightPrp.eq.ior(tem_fromFinerFace_prp, tem_remoteFace_prp)) then
      isFromFiner = tem_right

    ! remote+fromFiner-fluid face (data is send to the remote rank on this
    ! level. The remote rank will also do the interpolation for me, so my rank
    ! does not have to do anyhting here.)
    elseif( leftPrp.eq.ior(tem_fromFinerFace_prp, tem_remoteFace_prp) .and.    &
      &     rightPrp.eq.tem_fluidFace_prp ) then
      isFromFiner = tem_left

    ! remote+fromFiner-remote face (everything is remote)
    elseif( leftPrp.eq.ior(tem_fromFinerFace_prp, tem_remoteFace_prp) .and.    &
      &     rightPrp.eq.tem_remoteFace_prp ) then
      isFromFiner = 0

    ! remote-remote+fromFiner face (everything is remote)
    elseif( leftPrp.eq. tem_remoteFace_prp .and.                               &
      &     rightPrp.eq.ior(tem_remoteFace_prp, tem_fromFinerFace_prp)) then
      isFromFiner = 0

    ! notExist-remote+fromFiner face (nothing is on my rank)
    elseif( leftPrp.eq. tem_notExist_prp .and.                                 &
      &     rightPrp.eq.ior(tem_remoteFace_prp, tem_fromFinerFace_prp)) then
      isFromFiner = 0

    ! remote+fromFiner-notExist face (nothing is on my rank)
    elseif( leftPrp.eq.ior(tem_remoteFace_prp, tem_fromFinerFace_prp) .and.    &
      &     rightPrp.eq.tem_notExist_prp) then
      isFromFiner = 0

    ! fromFiner-fromCoarser face (everything is local, this is an intermediate
    ! face level (i.e. leveljumps of difference larger than 1 occure here) )
    elseif( leftPrp.eq.tem_fromFinerFace_prp .and.                             &
      &     rightPrp.eq.tem_fromCoarserFace_prp ) then
      isFromFiner = tem_left

    ! fromCoarser-fromFiner face (everything is local, this is an intermediate
    ! face level (i.e. leveljumps of difference larger than 1 occure here) )
    elseif( leftPrp.eq.tem_fromCoarserFace_prp .and.                           &
      &     rightPrp.eq.tem_fromFinerFace_prp ) then
      isFromFiner = tem_right

    ! fromFiner+remote - fromCoarser face (intermediate level for domain
    ! boundary and refinement boundary)
    elseif( leftPrp.eq.ior(tem_fromFinerFace_prp,tem_remoteFace_prp) .and.     &
      &     rightPrp.eq.tem_fromCoarserFace_prp) then
      isFromFiner = tem_left

    ! fromCoarser - remote+fromFiner face (intermediate level for domain
    ! boundary and refinement boundary)
    elseif( leftPrp.eq.tem_fromCoarserFace_prp .and.                           &
      &     rightPrp.eq.ior(tem_fromFinerFace_prp,tem_remoteFace_prp)) then
      isFromFiner = tem_right

    ! remote+fromFiner - remote+fromFiner face (everything is remote, nothing
    ! to be done)
    elseif( leftPrp.eq.ior(tem_fromFinerFace_prp,tem_remoteFace_prp)  .and.    &
      &     rightPrp.eq.ior(tem_fromFinerFace_prp,tem_remoteFace_prp)) then
      isFromFiner = 0

    ! remote+fromFiner - fromFiner face (everything is from finer, nothing to
    ! be done)
    elseif( leftPrp.eq.ior(tem_fromFinerFace_prp,tem_remoteFace_prp) .and.     &
      &     rightPrp.eq.tem_fromFinerFace_prp) then
      isFromFiner = 0

    ! fromFiner - remote+fromFiner face (everything is from finer, nothing to
    ! be done)
    elseif( leftPrp.eq.tem_fromFinerFace_prp  .and.                            &
      &     rightPrp.eq.ior(tem_fromFinerFace_prp,tem_remoteFace_prp)) then
      isFromFiner = 0

    ! fromFiner - remote face (nothing to be done)
    elseif(leftPrp.eq.tem_fromFinerFace_prp .and. rightPrp.eq.tem_remoteFace_prp) then
      isFromFiner = 0

    ! remote - fromFiner face (nothing to be done)
    elseif( leftPrp.eq.tem_remoteFace_prp .and.                                &
      &     rightPrp.eq.tem_fromFinerFace_prp) then
      isFromFiner = 0

    ! remote+fromCoarser - fromCoarser face (elements are on coarser
    ! level, here and remote so nothing to be done)
    elseif( leftPrp.eq.ior(tem_fromCoarserFace_prp,tem_remoteFace_prp) .and.   &
      &     rightPrp.eq.tem_fromCoarserFace_prp) then
      isFromFiner = 0

    ! fromCoarser - remote+fromCoarser face (elements are on coarser
    ! level, here and remote so nothing to be done)
    elseif( leftPrp.eq.tem_fromCoarserFace_prp .and.                           &
      &     rightPrp.eq.ior(tem_fromCoarserFace_prp,tem_remoteFace_prp)) then
      isFromFiner = 0

    ! fluid - boundary face (nothing to be done)
    elseif(leftPrp.eq.tem_fluidFace_prp .and. rightPrp.eq.tem_bndFace_prp) then
      isFromFiner = 0

    ! boundary - fluid face (nothing to be done)
    elseif(leftPrp.eq.tem_bndFace_prp .and. rightPrp.eq.tem_fluidFace_prp) then
      isFromFiner = 0

    ! remote - boundary face (nothing to be done)
    elseif(leftPrp.eq.tem_remoteFace_prp .and. rightPrp.eq.tem_bndFace_prp) then
      isFromFiner = 0

    ! boundary - remote face (nothing to be done)
    elseif(leftPrp.eq.tem_bndFace_prp .and. rightPrp.eq.tem_remoteFace_prp) then
      isFromFiner = 0

    ! fromFiner - boundary face (nothing to be done)
    elseif( leftPrp.eq.tem_fromFinerFace_prp .and.                             &
      &     rightPrp.eq.tem_bndFace_prp) then
      isFromFiner = 0

    ! boundary - fromFiner face (nothing to be done)
    elseif( leftPrp.eq.tem_bndFace_prp .and.                                   &
      &     rightPrp.eq.tem_fromFinerFace_prp) then
      isFromFiner = 0

    ! fromFiner+remote - boundary face (nothing to be done)
    elseif( leftPrp.eq.ior(tem_fromFinerFace_prp,tem_remoteFace_prp) .and.     &
      &     rightPrp.eq.tem_bndFace_prp) then
      isFromFiner = 0

    ! boundary - fromFiner+remote face (nothing to be done)
    elseif( leftPrp.eq.tem_bndFace_prp .and.                                   &
      &     rightPrp.eq.ior(tem_fromFinerFace_prp,tem_remoteFace_prp)) then
      isFromFiner = 0

    ! fromCoarser - boundary face (nothing to be done)
    elseif( leftPrp.eq.tem_fromCoarserFace_prp .and.                           &
      &     rightPrp.eq.tem_bndFace_prp) then
      isFromFiner = 0

    ! boundary - fromCoarser face (nothing to be done)
    elseif( leftPrp.eq.tem_bndFace_prp .and.                                   &
      &     rightPrp.eq.tem_fromCoarserFace_prp) then
      isFromFiner = 0

!> \todo JZ: Add more cases for the different types of faces. Especially
!! combinations of finer/coarser elements and remote properties are still
!! missing.
!! For some combinations it might be necessary to determine if vertical
!! dependencies exist or not.

    else
      write(logUnit(1),*) 'ERROR in tem_isFromFinerFace: not able to decide whether '// &
        &        'face is locally from finer or not, stopping'
      write(logUnit(1),*) 'left elem pos: ', faces%faceList%leftElemPos%val(facePos)
      write(logUnit(1),*) 'right elem pos: ', faces%faceList%rightElemPos%val(facePos)
      write(logUnit(1),*) 'left face prp: ', leftPrp
      write(logUnit(1),*) 'right face prp: ', rightPrp
      call tem_abort()
    end if

  end function tem_isFromFinerFace
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Extracts all the faces which will be computed by current rank.
  !!
  subroutine tem_build_computeList( faces, nEligibleChildren )
    ! --------------------------------------------------------------------------
    !> The created face descriptor.
    type(tem_face_type),intent(inout) :: faces
    !> The number of eligible children for the vertical face dependency
    integer, intent(in) :: nEligibleChildren
    ! --------------------------------------------------------------------------
    integer :: nComputeFaces(3)
    integer :: iDir, iFace, iComputeFace
    ! --------------------------------------------------------------------------


    ! First, we count all the element which will be covered by this rank.
    nComputeFaces = 0
    do iDir = 1,3

      do iFace = 1, faces%faces(iDir)%faceList%faceId%nVals

        ! decide if the element will be computed on this level and this rank.
        if ( tem_isComputeFace(iFace,             &
          &                    faces%faces(iDir), &
          &                    nEligibleChildren) ) then
          nComputeFaces(iDir) = nComputeFaces(iDir) + 1
        end if

      end do

    end do

    ! Allocate memory for the compute faces
    do iDir = 1,3
      allocate( faces%faces(iDir)%computeFace%leftPos( nComputeFaces(iDir) ) )
      allocate( faces%faces(iDir)%computeFace%rightPos( nComputeFaces(iDir) ) )
      allocate( faces%faces(iDir)%computeFace%facePos( nComputeFaces(iDir) ) )
    end do

    ! Now, we add all the compute elements to the iterator
    do iDir = 1,3
      iComputeFace = 1
      do iFace = 1, faces%faces(iDir)%faceList%faceId%nVals

        ! Decide if the element will be computed on this level and this rank.
        ! If yes, we add it to the list of compute faces.
        if ( tem_isComputeFace(iFace,             &
          &                    faces%faces(iDir), &
          &                    nEligibleChildren) ) then
          faces%faces(iDir)%computeFace%leftPos(iComputeFace) =                &
            &                 faces%faces(iDir)%faceList%leftElemPos%val(iFace)
          faces%faces(iDir)%computeFace%rightPos(iComputeFace) =               &
            &                faces%faces(iDir)%faceList%rightElemPos%val(iFace)
          faces%faces(iDir)%computeFace%facePos(iComputeFace) = iFace
          iComputeFace = iComputeFace + 1
        end if

      end do
    end do
  end subroutine tem_build_computeList
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> summary: Builds all the vertical dependencies in the face description.
  !!
  subroutine tem_faceDep_vertical(minLevel, maxLevel, faces, nEligibleChildren )
    ! --------------------------------------------------------------------------
    !> Minimum level of your mesh.
    integer, intent(in) :: minLevel
    !> Maximum level of your mesh.
    integer, intent(in) :: maxLevel
    !> Face descriptor where the faces will be appended to.
    type(tem_face_type),intent(inout)  :: faces(minLevel:maxLevel)
    !> The number of eligible children for the vertical face dependency
    integer, intent(in) :: nEligibleChildren
    ! --------------------------------------------------------------------------
    integer :: iLevel, iDir
    ! --------------------------------------------------------------------------

    ! Initialize the container for the vertical face dependencies.
    levelLoop: do iLevel = minLevel, maxLevel
      do iDir = 1, 3
        call tem_init_faceDep( faces(iLevel)%faces(iDir), nEligibleChildren )
      end do
    end do levelLoop

    ! build the dependencies downward (coarse face -> finer faces)
    levelLoopDown: do iLevel = minLevel, maxLevel-1
      do iDir = 1, 3
        call tem_faceDep_verticalDown( faces(iLevel)%faces(iDir),         &
          &                            faces(iLevel+1)%faces(iDir), iDir, &
          &                            nEligibleChildren )
      end do
    end do levelLoopDown

    ! build the dependencies upward (fine face -> coarse face)
    levelLoopUp: do iLevel = maxLevel, minLevel+1, -1
      do iDir = 1, 3
        call tem_faceDep_verticalUp( faces(iLevel)%faces(iDir),                &
          &                          faces(iLevel-1)%faces(iDir) )
      end do
    end do levelLoopUp

  end subroutine tem_faceDep_vertical
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> summary: Subroutine to initialize the container for the vertical
  !! dependencies of a face descriptor.
  !!
  subroutine tem_init_faceDep( faces, nEligibleChildren )
    ! --------------------------------------------------------------------------
    !> The face description in which you want to initialize the vertical
    !! face dependency container.
    type(tem_face_descriptor_type), intent(inout) :: faces
    !> The number of eligible children for the vertical face dependency
    integer, intent(in) :: nEligibleChildren
    ! --------------------------------------------------------------------------
    ! --------------------------------------------------------------------------

    allocate( faces%faceDep%parentFaceId(faces%faceList%faceId%nVals) )
    allocate( faces%faceDep%parentFacePos(faces%faceList%faceId%nVals) )
    allocate( faces%faceDep%childFaceId( nEligibleChildren, &
                                       & faces%faceList%faceId%nVals) )
    allocate( faces%faceDep%childFaceIdOp( nEligibleChildren, &
                                       & faces%faceList%faceId%nVals) )
    allocate( faces%faceDep%childFacePos( nEligibleChildren, &
                                       & faces%faceList%faceId%nVals) )
    allocate( faces%faceDep%childFacePosOp( nEligibleChildren, &
                                       & faces%faceList%faceId%nVals) )

    faces%faceDep%parentFaceId = -1_long_k
    faces%faceDep%parentFacePos = -1
    faces%faceDep%childFaceId = -1_long_k
    faces%faceDep%childFacePos(:,:) = -1

  end subroutine tem_init_faceDep
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> summary: Builds downward face depndencies (coarse face -> finer faces).
  !!
  subroutine tem_faceDep_verticalDown( coarseFaces, fineFaces, dir, &
                                     & nEligibleChildren )
    ! --------------------------------------------------------------------------
    !> Face description on the coarse level. The dependencies to the finer
    !! level will be appended to this face descriptor.
    type(tem_face_descriptor_type), intent(inout) :: coarseFaces
    !> Face description on the finer level.
    type(tem_face_descriptor_type), intent(in) :: fineFaces
    !> The spatial direction of the face to add the downward dependencies for.
    !! 1 --> x direction
    !! 2 --> y direction
    !! 3 --> z direction
    integer, intent(in) :: dir
    !> The number of eligible children for the vertical face dependency
    integer, intent(in) :: nEligibleChildren
    ! --------------------------------------------------------------------------
    integer :: iFace
    integer :: leftPrp, rightPrp
    ! --------------------------------------------------------------------------

    ! Run over all the faces and add the downward depndency (depending on the
    ! properties of the face).
    faceLoop: do iFace = 1, coarseFaces%faceList%faceId%nVals

      ! get its properties from left and right
      call tem_getFace_prp( coarseFaces, iFace, leftPrp, rightPrp )

      ! check if this face requires downwards dependencies. If yes, we just add
      ! this dependency to the coarse face descriptor.
      if( tem_reqDep_down(leftPrp, rightPrp) ) then
        call tem_addDep_down( iFace, coarseFaces, fineFaces, &
                            & dir, nEligibleChildren )
      end if

    end do faceLoop

  end subroutine tem_faceDep_verticalDown
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> summary: Builds upward face depndencies (fine face -> coarse faces).
  !!
  subroutine tem_faceDep_verticalUp( fineFaces, coarseFaces )
    ! --------------------------------------------------------------------------
    !> Face description on the fine level. The dependencies to the coarser
    !! level will be appended to this face descriptor.
    type(tem_face_descriptor_type), intent(inout) :: fineFaces
    !> Face description on the coarse level.
    type(tem_face_descriptor_type), intent(in) :: coarseFaces
    ! --------------------------------------------------------------------------
    integer :: iFace
    integer :: leftPrp, rightPrp
    ! --------------------------------------------------------------------------

    ! Run over all the faces and add the downward depndency (depending on the
    ! properties of the face).
    faceLoop: do iFace = 1, fineFaces%faceList%faceId%nVals

      ! get its properties from left and right
      call tem_getFace_prp( fineFaces, iFace, leftPrp, rightPrp )

      ! check if this face requires upwards dependencies. If yes, we just add
      ! this dependency to the coarse face descriptor.
      if( tem_reqDep_up(leftPrp, rightPrp) ) then
        call tem_addDep_up( iFace, fineFaces, coarseFaces )
      end if

    end do faceLoop

  end subroutine tem_faceDep_verticalUp
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Subroutine returns left and right propery of a face.
  !!
  subroutine tem_getFace_prp( faces, facePos, leftPrp, rightPrp)
    ! --------------------------------------------------------------------------
    !> The face descriptor the face is located in.
    type(tem_face_descriptor_type), intent(in) :: faces
    !> The position of the face in the face descriptor.
    integer, intent(in) :: facePos
    !> The left property of the face.
    integer, intent(out) :: leftPrp
    !> The right property of the face.
    integer, intent(out) :: rightPrp
    ! --------------------------------------------------------------------------
    ! --------------------------------------------------------------------------

    ! read the left and right property
    leftPrp = faces%faceList%leftPrp%val(facePos)
    rightPrp = faces%faceList%rightPrp%val(facePos)

  end subroutine tem_getFace_prp
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Function to check if a face with given left and right property requires
  !! a downward (coarse->fine) dependency.
  function tem_reqDep_down(leftPrp, rightPrp) result( reqDep )
    ! --------------------------------------------------------------------------
    !> Left property of the face.
    integer, intent(in) :: leftPrp
    !> Right property of the face.
    integer, intent(in) :: rightPrp
    !> Does the face requires a vertical dependency.
    logical :: reqDep
    ! --------------------------------------------------------------------------

    reqDep = .false.

    ! Check if one of the face properties (left or right) has the from finer
    ! property.
    if ( iand(tem_fromFinerFace_prp,leftPrp).ne.0 &
      & .or.                                      &
      & iand(tem_fromFinerFace_prp,rightPrp).ne.0 ) then
      reqDep = .true.
    end if

  end function tem_reqDep_down
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Function to check if a face with given left and right property requires
  !! a upward (fine->coarse) dependency.
  !!
  function tem_reqDep_up(leftPrp, rightPrp) result( reqDep )
    ! --------------------------------------------------------------------------
    !> Left property of the face.
    integer, intent(in) :: leftPrp
    !> Right property of the face.
    integer, intent(in) :: rightPrp
    !> Does the face requires a vertical dependency.
    logical :: reqDep
    ! --------------------------------------------------------------------------

    reqDep = .false.

    ! Check if one of the face properties (left or right) has the from coarser
    ! property.
    if ( iand(tem_fromCoarserFace_prp,leftPrp)/=0  &
      &  .or.                                      &
      &  iand(tem_fromCoarserFace_prp,rightPrp)/=0 ) then
      reqDep = .true.
    end if

  end function tem_reqDep_up
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> summary: Adds an downward face dependency (coarse->fine) to the coarse face
  !! descriptor.
  !!
  subroutine tem_addDep_down( coarseFacePos, coarseFaces, fineFaces, dir, &
                            & nEligibleChildren )
    ! --------------------------------------------------------------------------
    !> The position of the coarse face in coarseFaces you want to add
    !! the child dependencies for.
    integer, intent(in) :: coarseFacePos
    !> The description of the faces on the coarse level. The dependency will be
    !! added to this face descriptor.
    type(tem_face_descriptor_type), intent(inout) :: coarseFaces
    !> The descriptor of the faces on the fine level.
    type(tem_face_descriptor_type), intent(in) :: fineFaces
    !> The spatial direction of the face to add the downward dependencies for.
    !! 1 --> x direction
    !! 2 --> y direction
    !! 3 --> z direction
    integer, intent(in) :: dir
    !> The number of eligible children for the vertical face dependency
    integer, intent(in) :: nEligibleChildren
    ! --------------------------------------------------------------------------
    integer(kind=long_k) :: faceId, rightElemId, childFaceId, childFaceIdOp, &
                          & leftOf_childFaceIdOp
    integer :: leftOf_coord(4)
    integer(kind=long_k) :: children(8), childrenOp(8)
    integer, allocatable :: eligibleChildren(:), eligibleChildrenOp(:)
    integer :: eligChildPar, eligChildParOp, iChild, childFacePos, childFacePosOp
    ! --------------------------------------------------------------------------

    ! Left element id of the face.
    faceId = coarseFaces%faceList%faceId%val(coarseFacePos)
    rightElemId = coarseFaces%faceList%rightElemId%val(coarseFacePos)

    ! Identify the all the eligible child elements.
    select case(nEligibleChildren)
    case(4)
      select case(dir)
      case(1) ! x - direction
        eligChildPar = q__E
        eligChildParOp = q__W
      case(2) ! y - direction
        eligChildPar = q__N
        eligChildParOp = q__S
      case(3) ! z - direction
        eligChildPar = q__T
        eligChildParOp = q__B
      case default
        write(*,*) 'ERROR in tem_addDep_down: unknown spatial direction, '//     &
          &        'stopping ...'
        call tem_abort()
      end select
      call tem_eligibleChildren(eligibleChildren, eligChildPar)
      call tem_eligibleChildren(eligibleChildrenOp, eligChildParOp)
      children = tem_directChildren(faceId)
      childrenOp = tem_directChildren(rightElemId)
    case(8)
      allocate(eligibleChildren(nEligibleChildren))
      allocate(eligibleChildrenOp(nEligibleChildren))
      eligibleChildren =  (/ 1, 2, 3, 4, 5, 6, 7, 8 /)
      eligibleChildrenOp =  (/ 1, 2, 3, 4, 5, 6, 7, 8 /)
      children = tem_directChildren(faceId)
      childrenOp = tem_directChildren(rightElemId)
    case default
      write(logUnit(1),*) 'ERROR in tem_addDep_down: Unknown number of'
      write(logUnit(1),*) ' eligible children, stopping ...'
      call tem_abort()
    end select

    ! Loop over the relevant, eligible children.
    childLoop: do iChild = 1,nEligibleChildren

      ! The tree id of the current eligible child
      childFaceId = children(eligibleChildren(iChild))
      childFaceIdOp = childrenOp(eligibleChildrenOp(iChild))

      ! Get the TreeID which is left of childFaceIdOp for face identification
      leftOf_coord = tem_CoordOfId( childFaceIdOp  )
      leftOf_coord(dir) = leftOf_coord(dir)-1
      leftOf_childFaceIdOp = tem_IdOfCoord(leftOf_coord)

      ! Get its position in the fine face descriptor.
      childFacePos = PositionOfVal( me  = fineFaces%faceList%faceId,           &
        &                           val = childFaceId )

      ! Get the op position in the fine face descriptor.
      childFacePosOp = PositionOfVal( me  = fineFaces%faceList%faceId,           &
        &                           val = leftOf_childFaceIdOp )

      ! If the child element exists in the fine face descriptor, add a
      ! dependency from coarse element to its child element.
      if(childFacePos<=0) then
        coarseFaces%faceDep%childFaceId(iChild, coarseFacePos) = -1_long_k
        coarseFaces%faceDep%childFacePos(iChild, coarseFacePos) = -1
      else
        coarseFaces%faceDep%childFaceId(iChild, coarseFacePos) = childFaceId
        coarseFaces%faceDep%childFacePos(iChild, coarseFacePos) = childFacePos
      end if

      if(childFacePosOp<=0) then
        coarseFaces%faceDep%childFaceIdOp(iChild, coarseFacePos) = -1_long_k
        coarseFaces%faceDep%childFacePosOp(iChild, coarseFacePos) = -1
      else
        coarseFaces%faceDep%childFaceIdOp(iChild, coarseFacePos) = leftOf_childFaceIdOp
        coarseFaces%faceDep%childFacePosOp(iChild, coarseFacePos) = childFacePosOp
      end if

    end do childLoop

  end subroutine tem_addDep_down
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> summary: Adds an upward face dependency (fine->coarse) to the fine face
  !! descriptor.
  !!
  subroutine tem_addDep_up( fineFacePos, fineFaces, coarseFaces )
    ! --------------------------------------------------------------------------
    !> The position of the face in the fineFaces descriptor you want to add.
    integer, intent(in) :: fineFacePos
    !> The description of the faces on the fine level. The dependency will be
    !! added to this face descriptor.
    type(tem_face_descriptor_type), intent(inout) :: fineFaces
    !> The descriptor of the faces on the coarse level.
    type(tem_face_descriptor_type), intent(in) :: coarseFaces
    ! --------------------------------------------------------------------------
    integer(kind=long_k) :: faceId, parentId
    integer :: parentPos
    ! --------------------------------------------------------------------------

    ! Left element id of the face.
    faceId = fineFaces%faceList%faceId%val(fineFacePos)

    ! Identify the parent element.
    parentId = tem_parentOf(faceId)

    ! Get its position in the coarse face descriptor.
    parentPos = PositionOfVal( me  = coarseFaces%faceList%faceId, &
      &                        val = parentId                     )

    ! If the parent element exists in the coarse face descriptor, add a
    ! dependency from fine element to its parent element.
    if( parentPos<=0 ) then
      fineFaces%faceDep%parentFaceId(fineFacePos) = -1_long_k
      fineFaces%faceDep%parentFacePos(fineFacePos) = -1
    else
      fineFaces%faceDep%parentFaceId(fineFacePos) = parentId
      fineFaces%faceDep%parentFacePos(fineFacePos) = parentPos
    end if

  end subroutine tem_addDep_up
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> summary: Collects all the faces which have at least one neighboring
  !! element that exists in the level descriptor.
  !!
  subroutine tem_collect_faces(levelDesc, minLevel, maxLevel, faces )
    ! --------------------------------------------------------------------------
    !> Minimum level of your mesh.
    integer, intent(in) :: minLevel
    !> Maximum level of your mesh.
    integer, intent(in) :: maxLevel
    !> Level descriptor for each spatial direction and each level of your mesh.
    type(tem_levelDesc_type), intent(in) :: levelDesc(1:3,minLevel:maxLevel)
    !> Face descriptor where the faces will be appended to.
    type(tem_face_type),intent(inout)  :: faces(minLevel:maxLevel)
    ! --------------------------------------------------------------------------
    integer :: iLevel, iDir
    ! --------------------------------------------------------------------------

    levelLoop: do iLevel = minLevel, maxLevel
      directionLoop: do iDir = 1, 3

        ! For this level and this direction we collect all the faces.
        call tem_get_faces( levelDesc = levelDesc(iDir, iLevel),  &
          &                 direction = iDir,                     &
          &                 faces     = faces(iLevel)%faces(iDir) )

      end do directionLoop
    end do levelLoop

  end subroutine tem_collect_faces
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Subroutine to initialize the dyanmic content of a face list.
  !!
  subroutine tem_init_faceList( faceList )
    ! --------------------------------------------------------------------------
    !> The face list to initialize.
    type(tem_faceList_type), intent(inout) :: faceList
    ! --------------------------------------------------------------------------

    call init( me = faceList%faceId, length = 32 )
    call init( me = faceList%rightElemId, length = 32 )
    call init( me = faceList%leftElemPos, length = 32 )
    call init( me = faceList%rightElemPos, length = 32 )
    call init( me = faceList%leftPrp, length = 32 )
    call init( me = faceList%rightPrp, length = 32 )

  end subroutine tem_init_faceList
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> summary: Collects all the faces on a given level of the mesh in a given
  !! direction. It attaches also two-sided properties to the face.
  !!
  subroutine tem_get_faces(levelDesc, direction, faces)
    ! --------------------------------------------------------------------------
    !> Level descriptor of the level of the mesh you want to collect the
    !! faces for.
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> The spatial direction to collect the faces for:
    !! 1 -> x direction
    !! 2 -> y direction
    !! 3 -> z direction
    integer, intent(in) :: direction
    !> Description of the faces on the current level.
    type(tem_face_descriptor_type), intent(inout) :: faces
    ! --------------------------------------------------------------------------
    integer :: iElem
    integer :: neighIndex
    integer :: neighPos
    integer :: elemPrp, neighPrp
    integer(kind=long_k) :: elemId, neighId
    ! --------------------------------------------------------------------------

    ! Initialize the facelist (with a dynamic array for the faceIDs)
    call tem_init_faceList( faces%faceList )

    ! Loop over all the elements (includes fluid, ghost and halos) have a look
    ! to the left and to the right face.
    ! Try to add both faces of the element to the list of all faces.
    elemLoop: do iElem = 1, levelDesc%nElems

      ! Get the current element property
      elemPrp = tem_get_elemPrp(levelDesc, iElem)
      elemId = levelDesc%total(iElem)

      ! Look for the left face of the current element.
      ! We use the compute stencil to determine the position of the neighboring
      ! left element.
      neighIndex = tem_left
      call tem_get_faceNeigh( levelDesc, iElem, direction, neighIndex, &
        &                     neighId, neighPos                        )
      neighPrp = tem_get_elemPrp(levelDesc, neighPos)

      ! For the left face, the current element will be on its right side.
      call tem_addFace( leftElemId   = neighId,  &
        &               leftElemPos  = neighPos, &
        &               leftElemPrp  = neighPrp, &
        &               rightElemId  = elemId,   &
        &               rightElemPos = iElem,    &
        &               rightElemPrp = elemPrp,  &
        &               faces        = faces     )

      ! Look for the right face of the current element.
      ! Get the right neighbor of the element.
      neighIndex = tem_right
      call tem_get_faceNeigh( levelDesc, iElem, direction, neighIndex, &
        &                     neighId, neighPos                        )
      neighPrp = tem_get_elemPrp(levelDesc, neighPos)

      ! For the right face, the current element will be on its left side.
      call tem_addFace( leftElemId   = elemId,   &
        &               leftElemPos  = iElem,    &
        &               leftElemPrp  = elemPrp,  &
        &               rightElemId  = neighId,  &
        &               rightElemPos = neighPos, &
        &               rightElemPrp = neighPrp, &
        &               faces        = faces     )

    end do elemLoop

  end subroutine tem_get_faces
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Function to get the face neighbor of a certain element in the
  !! level descriptor. Even non-existing face neighbors can be handled
  !! by this routine.
  !!
  subroutine tem_get_faceNeigh( levelDesc, elemPos, dir, leftOrRight, neighId, &
    &                           neighPos )
    ! --------------------------------------------------------------------------
    !> The level descriptor the element is located in.
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> The element position in the level descriptor.
    integer, intent(in) :: elemPos
    !> The spatial direction of the face.
    !!
    !! 1. -> x direction
    !! 2. -> y direction
    !! 3. -> z direction
    integer, intent(in) :: dir
    !> Find left or right face neighbor of the element.
    !! Use
    !! \ref tem_face_module::tem_left
    !! or
    !! \ref tem_face_module::tem_right.
    integer, intent(in) :: leftOrRight
    !> The treeid of the face neighbor element.
    integer(kind=long_k), intent(out) :: neighId
    !> Position of the neighbor in the total list, or 0 if element does not
    !! exist.
    integer, intent(out) :: neighPos
    ! --------------------------------------------------------------------------
    integer(kind=long_k) :: elemId
    integer :: elemCoord(4)
    integer :: dirOffset
    ! --------------------------------------------------------------------------

    ! We check the neighbor position of the element.
    neighPos = levelDesc%neigh(1)%nghElems( leftOrRight, elemPos )
    elemId = levelDesc%total(elemPos)

    ! Can we find the neighbor in the level descriptor. If not, we have
    ! to build the tree id our own.
    if (neighPos>0 .and. neighPos<=levelDesc%nElems) then
      neighId = levelDesc%total(neighPos)
    else
      ! So, the neighbor element does not exist in level descriptor.
      ! Let's create the tree id (of the neighbor on this level) by our own.
      elemCoord = tem_coordOfId(elemId)
      select case(leftOrRight)
      case(tem_left)
        dirOffset = -1
      case(tem_right)
        dirOffset = 1
      case default
        write(*,*) 'ERROR in tem_get_faceNeigh: unknown face neighbor '//      &
          &        'direction stopping ... '
        write(*,*) 'Given direction was: ', dir
        call tem_abort()
      end select
      elemCoord(dir) = elemCoord(dir) + dirOffset
      neighId = tem_idOfCoord(elemCoord)

      ! Check if it is boundary element or not
      if (neighPos<0) then
        ! In case of a boundary element, we keep the negative boundary
        ! unmodified.
      else
        ! If no, we try to locate the neighbor in the total list of the level
        ! descriptor
        neighPos = tem_treeIDinTotal( neighId, levelDesc)
      end if
    endif

  end subroutine tem_get_faceNeigh
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Adds a new face to the face description.
  !!
  !! If the face already exists in the face description,
  !! only the properties of the already existing face will be overwritten.
  subroutine tem_addFace( leftElemId, leftElemPos, leftElemPrp,          &
    &                     rightElemId, rightElemPos, rightElemPrp, faces )
    ! --------------------------------------------------------------------------
    !> Element id of the left element
    integer(kind=long_k), intent(in) :: leftElemId
    !> Element id of the right element
    integer(kind=long_k), intent(in) :: rightElemId
    !> Position of the left element in the level descriptor's total list.
    integer, intent(in) :: leftElemPos
    !> Properties of the left element.
    integer, intent(in) :: leftElemPrp
    !> Position of the right element in the level desriptor's total list.
    integer, intent(in) :: rightElemPos
    !> Properties of the right element
    integer, intent(in) :: rightElemPrp
    !> The face description the new face will be added to. If the face already
    !! exists in this face description. The existing face's property will be
    !! overwritten by the new ones.
    type(tem_face_descriptor_type), intent(inout) :: faces
    ! --------------------------------------------------------------------------
    logical :: wasAdded
    integer :: pos
    ! --------------------------------------------------------------------------


    ! Try to add the left tree id to the face id list of the faces
    call append( me       = faces%faceList%faceId, &
      &          val      = leftElemId,            &
      &          wasAdded = wasAdded,              &
      &          pos      = pos                    )

    ! Check if the face already existed in the face descriptor.
    ! If yes, we overwrite the properties. If no, we overwrite
    ! the existing properties.
    if(wasAdded) then ! New face
      ! Append right element, element positions and new properties
      call append( me = faces%faceList%rightElemId, val = rightElemId )
      call append( me = faces%faceList%leftElemPos, val = leftElemPos )
      call append( me = faces%faceList%rightElemPos, val = rightElemPos )
      call append( me = faces%faceList%leftPrp, val = leftElemPrp )
      call append( me = faces%faceList%rightPrp, val = rightElemPrp )
    else ! Face already existed
      ! Overwrite existing properies
      faces%faceList%leftPrp%val(pos) = leftElemPrp
      faces%faceList%rightPrp%val(pos) = rightElemPrp
    end if

  end subroutine tem_addFace
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> summary: Extends the properties of the faces which are marked as
  !! communication faces (at least by one of the neighboring elements). It
  !! checks if the communicated elements are from finer/coarser level and
  !! attaches the fromFiner/fromCoarser property to this face if necessary.
  !!
  subroutine tem_extend_remotePrp(levelDesc, minLevel, maxLevel, faces )
    ! --------------------------------------------------------------------------
    !> Minimum level of your mesh.
    integer, intent(in) :: minLevel
    !> Maximum level of your mesh.
    integer, intent(in) :: maxLevel
    !> Level descriptor for each level of your mesh.
    type(tem_levelDesc_type), intent(in) :: levelDesc(1:3,minLevel:maxLevel)
    !> The face descriptor to be corrected.
    type(tem_face_type),intent(inout)  :: faces(minLevel:maxLevel)
    ! --------------------------------------------------------------------------
    integer :: iLevel, iDir
    ! --------------------------------------------------------------------------

    ! We correct the properties for each level and each direction independently
    levelLoop: do iLevel = minLevel, maxLevel
      directionLoop: do iDir = 1, 3

        ! Check for from finer property of the communicated faces.
        call tem_extend_commFromFinerPrp( levelDesc(iDir, iLevel), iDir, &
          &                               faces(iLevel)%faces(iDir)      )

        ! Check for from coarser property of the communicated faces.
        call tem_extend_commFromCoarserPrp( levelDesc(iDir, iLevel), iDir, &
          &                                 faces(iLevel)%faces(iDir)      )

      end do directionLoop
    end do levelLoop

  end subroutine tem_extend_remotePrp
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> summary: Extend communication property for faces by the finer property if
  !! the neighboring halo element is from finer on the remote partition.
  !!
  subroutine tem_extend_commFromFinerPrp( levelDesc, direction, faces )
    ! --------------------------------------------------------------------------
    !> Level descriptor of the level of the mesh you want to collect the
    !! faces for.
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> The face direction.
    integer, intent(in) :: direction
    !> Description of the faces on the current level.
    type(tem_face_descriptor_type), intent(inout) :: faces
    ! --------------------------------------------------------------------------
    integer :: iProc, iElem
    integer :: neighIndex
    integer(kind=long_k) :: elemId, neighId
    integer :: elemPos, neighPos
    integer :: elemPrp
    ! --------------------------------------------------------------------------

    ! We run over all the elements in the from finer buffer and attach
    ! the from finer property to the face the element is aligned to.
    processLoop: do iProc = 1, levelDesc%recvbufferFromFiner%nProcs

      elementLoop: do iElem = 1, levelDesc%recvbufferFromFiner%nElemsProc(iProc)

        ! The element position in the total list of the level descriptor.
        elemPos = levelDesc%recvbufferFromFiner%elemPos(iProc)%val(iElem)
        elemId = levelDesc%total(elemPos)
        elemPrp = tem_get_elemPrp(levelDesc, elemPos)

        ! Attach elem property for its right face (checks also if face exists).
        neighIndex = tem_right
        call tem_get_faceNeigh( levelDesc, elemPos, direction, neighIndex, &
          &                     neighId, neighPos                          )
        call tem_appendFace_prp( elemId, tem_fromFinerFace_prp, faces, &
          &                      tem_left                      )

        ! Attach elem property for its left face (checks also if face exists).
        neighIndex = tem_left
        call tem_get_faceNeigh( levelDesc, elemPos, direction, neighIndex, &
          &                     neighId, neighPos                          )
        call tem_appendFace_prp( neighId, tem_fromFinerFace_prp, faces, &
          &                      tem_right                      )

    end do elementLoop
  end do processLoop

  end subroutine tem_extend_commFromFinerPrp
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> summary: Extend communication property for faces by the coarser property if
  !! the neighboring halo element is from coarser on the remote partition.
  !!
  subroutine tem_extend_commFromCoarserPrp( levelDesc, direction, faces )
    ! --------------------------------------------------------------------------
    !> Level descriptor of the level of the mesh you want to collect the
    !! faces for.
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> The face direction.
    integer, intent(in) :: direction
    !> Description of the faces on the current level.
    type(tem_face_descriptor_type), intent(inout) :: faces
    ! --------------------------------------------------------------------------
    integer :: iElem, iProc
    integer :: neighIndex
    integer(kind=long_k) :: elemId, neighId
    integer :: elemPos, neighPos
    integer :: elemPrp
    ! --------------------------------------------------------------------------

    ! We run over all the elements in the from coarser buffer and attach
    ! the from coarser property to the face the element is aligned to.
    processLoop: do iProc = 1, levelDesc%recvbufferFromCoarser%nProcs

      elementLoop: do iElem = 1, levelDesc%recvbufferFromCoarser%              &
        &                                                     nElemsProc(iProc)

        ! The element position in the total list of the level descriptor.
        elemPos = levelDesc%recvbufferFromCoarser%elemPos(iProc)%val(iElem)
        elemId = levelDesc%total(elemPos)
        elemPrp = tem_get_elemPrp(levelDesc, elemPos)

        ! Attach elem property for its right face (checks also if face exists).
        neighIndex = tem_right
        call tem_get_faceNeigh( levelDesc, elemPos, direction, neighIndex, &
          &                     neighId, neighPos                          )
        call tem_appendFace_prp( elemId, tem_fromCoarserFace_prp, faces, &
          &                      tem_left                        )

        ! Attach elem property for its left face (checks also if face exists).
        neighIndex = tem_left
        call tem_get_faceNeigh( levelDesc, elemPos, direction, neighIndex, &
          &                     neighId, neighPos                          )
        call tem_appendFace_prp( neighId, tem_fromCoarserFace_prp, faces, &
          &                      tem_right                        )

      end do elementLoop
    end do processLoop

  end subroutine tem_extend_commFromCoarserPrp
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Attaches another property to a given face (from left or right). If the face
  !! does not exist this routine will not do anything.
  !!
  subroutine tem_appendFace_prp( faceId, prp, faces, prp_dir )
    ! --------------------------------------------------------------------------
    !> The face identifier to append the face property to.
    integer(kind=long_k), intent(in) :: faceId
    !> The property to attach
    integer, intent(in) :: prp
    !> The face desriptor to update.
    type(tem_face_descriptor_type), intent(inout) :: faces
    !> Attach the property to the left or right property of the face.
    integer, intent(in) :: prp_dir
    ! --------------------------------------------------------------------------
    integer :: facePos, currentPrp, newPrp
    ! --------------------------------------------------------------------------

    ! Check if the face element exists and attach the property.
    facePos = PositionOfVal( faces%faceList%faceId, faceId )

    ! If the element is found: Check if the element has already the
    ! property we want to attach to it. If not, add it.
    if(facePos.gt.0) then

      ! Get the current face propery.
      if( prp_dir .eq. tem_left ) then
        currentPrp = faces%faceList%leftPrp%val(facePos)
      else
        currentPrp = faces%faceList%rightPrp%val(facePos)
      end if

      ! Melt the property we want to attach with the current one.
      newPrp = tem_melt_facePrp(currentPrp, prp)

      ! Set the melted face propery.
      if( prp_dir .eq. tem_left ) then
        faces%faceList%leftPrp%val(facePos) = newPrp
      else
        faces%faceList%rightPrp%val(facePos) = newPrp
      end if

    end if

  end subroutine tem_appendFace_prp
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Function to melt two properties together. The resulting property
  !! holds the union of firstPrp and secondPrp.
  !!
  function tem_melt_facePrp(firstPrp, secondPrp) result(meltedPrp)
    ! --------------------------------------------------------------------------
    !> The first property to be melted.
    integer, intent(in) :: firstPrp
    !> The second property to be melted.
    integer, intent(in) :: secondPrp
    !> The resulting property (the union of the first and second property)
    integer :: meltedPrp
    ! --------------------------------------------------------------------------

    ! Apply a bitwise OR operation to melt the two properties.
    meltedPrp = ior(firstPrp, secondPrp)

  end function tem_melt_facePrp
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> summary: Returns the property of a given element in the level descriptor
  !! which can be assigned to one of its face properties.
  !!
  function tem_get_elemPrp(levelDesc, elemPos) result(elemPrp)
    ! --------------------------------------------------------------------------
    !> The level descriptor that contains the investigated element.
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> The position of the element in the total list of the level descriptor.
    integer, intent(in) :: elemPos
    !> The element property of the element.
    integer :: elemPrp
    ! --------------------------------------------------------------------------

    elemPrp = 0

    ! If the element does not exist, we can only assign the non-existing
    ! property.
    if(elemPos == 0) then
      elemPrp = tem_notExist_prp
      return
    elseif(elemPos < 0) then
      elemPrp = tem_bndFace_prp
      return
    end if

    ! Check for fluid
    if( levelDesc%offset(1,1)+1 <= elemPos .and.                               &
      & elemPos<=levelDesc%offset(1,1)+levelDesc%elem%nElems( eT_fluid )  ) then
      elemPrp = tem_fluidFace_prp
    ! Check for ghost from coarser
    else if (levelDesc%offset(1,2)+1<=elemPos .and.                          &
      & elemPos<=levelDesc%offset(1,2)+levelDesc%elem%nElems( eT_ghostFromCoarser) ) then
      elemPrp = tem_fromCoarserFace_prp
    ! Check for ghost from finer
    else if (levelDesc%offset(1,3)+1<=elemPos .and.                          &
      & elemPos<=levelDesc%offset(1,3)+levelDesc%elem%nElems( eT_ghostFromFiner )) then
      elemPrp = tem_fromFinerFace_prp
    ! Check for halo
    else if (levelDesc%offset(1,4)+1<=elemPos .and.                          &
      & elemPos<=levelDesc%offset(1,4)+levelDesc%elem%nElems( eT_halo ) ) then
      elemPrp = tem_remoteFace_prp
    ! Everything else should not occur
    else
      write(*,*) 'ERROR in tem_get_elemPrp: Unknown element property, '//      &
        &        'stopping ...'
      call tem_abort()
    end if

  end function tem_get_elemPrp
  ! ************************************************************************** !


end module tem_face_module
! **************************************************************************** !
