tem_shape_findElemByBCLabels Subroutine

  subroutine tem_shape_findElemByBCLabels( bcLabels, cutOffQVal, bc_prop,     &
    &                                      foundAny, map2global, inTree,      &
    &                                      countPoints, grwPnts, storePnts,   &
    &                                      bcIDs, stencil )
    ! --------------------------------------------------------------------------
    !> bcLabels
    character(len=labelLen), intent(in) :: bcLabels(:)
    !> Cut off qValue
    real(kind=rk), intent(in) :: cutOffQVal
    !> bc property
    type(tem_bc_prop_type), intent(in) :: bc_prop
    !> if any element be identified
    logical,                 intent(out) :: foundAny
    !> dynamic array. Elements positions in bc_prop%property%elemID
    type(dyn_intArray_type), intent(inout) :: map2global
    !> Global tree
    type(treelmesh_type), intent(in) :: inTree
    !> How many points there will be
    integer, intent( inout ) :: countPoints
    !> growing array to store tracking points
    type(tem_grwPoints_type), intent(inout) :: grwPnts
    !> to Store points in grwPnts
    logical, intent(in) :: storePnts
    !> id of boundary condition to be tracked
    integer, allocatable, intent(out) :: bcIDs(:)
    !> stencil required to get useful links
    type( tem_stencilHeader_type ), optional, intent(in) :: stencil
    ! --------------------------------------------------------------------------
    integer :: dPos, iElem, iBC, nBCs, iBCtype, posInTree, QQN
    logical :: wasAdded, found
    integer, allocatable :: map(:)
    real(kind=rk), allocatable :: qValWeight(:), weightedQVal(:)
    real(kind=rk) :: dx, bary(3), point(3), minQVal
    integer :: iDir, length, iElem_qVal, minQValiDir
    ! --------------------------------------------------------------------------

    if( present( stencil ) ) then
      allocate( map( size(stencil%map) ) )
      map = stencil%map
      qqn = stencil%QQN
    else
      allocate( map(qQQQ) )
      map = qDir
      qqn = qQQQ
    end if

    ! calculate weights on qValue for each direction. 
    ! It depends on the link length.
    allocate(weightedQVal(qqn))
    allocate(qValWeight(qqn))
    do iDir = 1, qqn
      length   = qOffset( map(iDir), 1 )**2 &
        &      + qOffset( map(iDir), 2 )**2 &
        &      + qOffset( map(iDir), 3 )**2

      qValWeight(iDir) = sqrt(real(length, kind=rk))
    end do

    foundAny = .false.
    ! convert bcLabels to bcIDs
    nBCs = size(bcLabels)
    allocate( bcIDs( nBCs ) )

    ! loop over all required bcLabels
    do iBC = 1, nBCs

      found = .false.

      ! loop over all BCs in the mesh
      do iBCtype = 1, bc_prop%nBCtypes
        if ( trim(bcLabels(iBC)) == bc_prop%bc_label( iBCtype ) ) then
          bcIDs( iBC ) = iBCtype
          found = .true.
          exit
        end if
      end do

      if ( .not. found ) then
        write(logUnit(1),*) 'Required BC label: '//trim(bcLabels(iBC))
        write(logUnit(1),*) 'can not be found in given mesh!'
        stop
      end if
    end do ! iBC

    iElem_qVal = 0
    ! Loop over all element with boundary property
    do iElem = 1, bc_prop%property%nElems
      posInTree = bc_prop%property%elemID(iElem)
      ! Count the elem with qVal property to access the qVal from
      ! bc_prop%qVal
      if ( btest( inTree%elemPropertyBits( posInTree ), prp_hasQVal) ) then
        iElem_qVal = iElem_qVal + 1
      end if

      do iBC = 1, nBCs
        if ( any(int(bc_prop%boundary_ID( map(:qqn), iElem )) == bcIDs(iBC) ) ) then

          ! If element has qValues then consider element which satisfy
          ! cutOffQVal
          if (btest(inTree%ElemPropertyBits(posInTree), prp_hasQVal)) then
            do iDir = 1, qqn
              weightedQVal(iDir) = qValWeight(iDir)                  &
                &                * bc_prop%qVal(map(iDir), iElem_qVal)
            end do
            minQValiDir = minloc(weightedQVal, DIM=1, MASK=(weightedQVal>0))
            minQVal = bc_prop%qVal(map(minQValiDir), iElem_qVal)
          else
            minQVal = cutOffQVal
          end if

          if (minQVal <= cutOffQVal) then

            ! Append to treeID list (note that already existing ones are
            ! omitted)
            call append( me       = map2global,   &
              &          pos      = dPos,         &
              &          val      = posInTree,    &
              &          wasAdded = wasAdded )

            ! Count up if it was added
            if( wasAdded ) then
              ! tLevel = tem_levelOf( inTree%treeID(posInTree) )
              ! countElems( tLevel ) = countElems( tLevel ) + 1
              foundAny = .true.

              if (storePnts) then
                ! Create growing array of points for boundary elements with 
                ! qValues.
                ! Point is created along minimum qValue direction.
                bary = tem_BaryOfId( inTree, inTree%treeID(posInTree) )
                if (btest(inTree%ElemPropertyBits(posInTree), prp_hasQVal)) then
                  dx = tem_ElemSize( inTree, inTree%treeID(posInTree) ) 
                  point = bary + dx * qOffset(map(minQValiDir), :)
                  ! append the physical points to the growing array of points
                  call append( me  = grwPnts, &
                    &          val = point    )
                else
                  call append( me  = grwPnts, &
                    &          val = bary     )
                end if
                countPoints = countPoints + 1
              end if
            end if !storePnts

            exit ! continue with next element
          end if ! wasAdded

        end if ! boundary_ID == bcIDs
      end do

    end do ! iElem

  end subroutine tem_shape_findElemByBCLabels