tem_cano_storePntsInSubTree Subroutine

  subroutine tem_cano_storePntsInSubTree( me, inTree, map2global, countPoints, &
    &                                     grwPnts )
    ! --------------------------------------------------------------------------
    !> canonicalND objects on which to work
    type(tem_canonicalND_type ),intent(in) :: me(:)
    !> Global tree
    type(treelmesh_type), intent(in) :: inTree
    !> How many points there will be
    integer, intent(inout) :: countPoints
    !> growing array for the map2global
    type(dyn_intArray_type), intent(in) :: map2global
    !> growing array to store tracking points
    type(tem_grwPoints_type), intent(inout) :: grwPnts
    ! --------------------------------------------------------------------------
    integer :: nElems, nPoints, maxLevel, elemPos, neighPos, coordOfId(4)
    integer :: iCano, iPnt, iQQN
    real(kind=rk) :: coord(3), offset_a, offset_b, offset_c
    real(kind=rk) :: unit_vec_a(3),unit_vec_b(3), unit_vec_c(3)
    integer(kind=long_k) :: treeID, tOffset, neighID
    integer(kind=long_k), allocatable :: subTreeID(:)
    ! --------------------------------------------------------------------------
    maxLevel = inTree%global%maxLevel
    ! Append the physical points to the growing array of points
    nElems = map2global%nVals
    allocate(subTreeID(nElems))
    subTreeID = inTree%treeID(map2global%val(map2global%sorted(1:nElems)))
    do iCano = 1, size(me)
      ! total number of elements
      nPoints = me(iCano)%segments(1) &
        &     * me(iCano)%segments(2) &
        &     * me(iCano)%segments(3)

      unit_vec_a =  me(iCano)%vec(:,1)                      &
        &        / real( max(me(iCano)%segments(1)-1,1), rk )
      unit_vec_b =  me(iCano)%vec(:,2)                      &
        &        / real( max(me(iCano)%segments(2)-1,1), rk )
      unit_vec_c =  me(iCano)%vec(:,3)                      &
        &        / real( max(me(iCano)%segments(3)-1,1), rk )

      ! Generate points and append only the points available in tree
      do iPnt = 1, nPoints
        offset_a = real(mod((iPnt-1),me(iCano)%segments(1)), kind=rk)
        offset_b = real(mod((iPnt-1)/me(iCano)%segments(1), &
          &                 me(iCano)%segments(2)), kind=rk)
        offset_c = real(mod((iPnt-1)/(me(iCano)%segments(1)  &
          &                          *me(iCano)%segments(2)), &
          &                 me(iCano)%segments(3)), kind=rk)

        coord = me(iCano)%origin + offset_a * unit_vec_a &
          &                      + offset_b * unit_vec_b &
          &                      + offset_c * unit_vec_c

        ! Get the treeID on the highest level
        treeID = tem_IdOfCoord( tem_CoordOfReal(inTree, coord(1:3), maxLevel) )
        ! get position of the treeID in subTree
        elemPos = tem_PosOfId( treeID, subTreeID )

        neighPos = 0
        if (elemPos <= 0) then
          ! Point must be outside fluid domain, its neighbor must be in subTreeID
          ! since it was added in tem_cano_initSubTree
          coordOfId = tem_CoordOfId( treeID )
          tOffset = tem_FirstIdAtLevel( coordOfId(4) )
          directionLoop: do iQQN = 1, qQQQ
            neighID = tem_IdOfCoord(                          &
              &         [ coordOfId(1) + qOffset( iQQN, 1 ),  &
              &           coordOfId(2) + qOffset( iQQN, 2 ),  &
              &           coordOfId(3) + qOffset( iQQN, 3 ),  &
              &           coordOfId(4) ], tOffset)
            neighPos = tem_PosOfId( neighID, subTreeID)
            if (neighPos > 0) exit directionLoop
          end do directionLoop
        end if

        if( elempos > 0 .or. neighPos > 0 ) then
          ! append the physical points to the growing array of points
          call append( me  = grwpnts, &
            &          val = coord    )

          countPoints = countPoints + 1
        end if
      end do !iPoint
    end do !iCano
    deallocate(subTreeID)

  end subroutine tem_cano_storePntsInSubTree