tem_CoordOfReal Function

  pure function tem_CoordOfReal(mesh, point, level) result(coord)
    ! -------------------------------------------------------------------- !
    type(treelmesh_type), intent(in) :: mesh !< Mesh to locate the point in
    real(kind=rk), intent(in) :: point(3) !< Point to look up
    integer, intent(in), optional :: level !< optional level to return the
                                           !! coordinate on
    integer :: coord(4) !< x,y,z,level
    ! -------------------------------------------------------------------- !
    real(kind=rk) :: locInCube(3)
    real(kind=rk) :: meshDensity
    integer :: dimLen
    integer :: coordlevel
    ! -------------------------------------------------------------------- !

    if (present(level)) then
      coordlevel = min(level, globalMaxLevels)
    else
      coordlevel = globalMaxLevels
    end if

    locInCube = point - mesh%global%Origin
    dimLen = 2**coordlevel
    meshDensity = real(dimLen, kind=rk) / mesh%global%BoundingCubeLength

    ! Look up the real coordinate on the finest possible resolution.
    coord(4) = coordlevel

    ! Coordinates range from 0 to dimLen-1 in each direction
    ! Do not use periodic domain here, instead use the first
    ! and last element to capture all points outside the
    ! domain to deal with numerical inaccuracies.
    ! That is all elements include their lower boundaries, except
    ! the last one in a given direction, which includes the lower
    ! as well as the upper.
    coord(1:3) = max( min( int(locInCube*meshDensity), dimLen-1 ), 0 )

  end function tem_CoordOfReal