FaceNeighbor_QuadTreeMesh2D Subroutine

  recursive subroutine FaceNeighbor_QuadTreeMesh2D(this,node,s,nbr,ns,nf)
    !! Find the equal-or-larger face neighbour of `node` across its local side s using the
    !! classic quadtree ascend/descend search. Returns:
    !!   nbr - neighbour node id (0 = physical domain boundary). It is either a LEAF at any level
    !!         <= level(node), or an INTERNAL node at exactly level(node) (meaning the shared face
    !!         is subdivided on the neighbour side, i.e. finer neighbours exist).
    !!   ns  - the neighbour's local side that faces `node`.
    !!   nf  - the flip between the two shared edges (0 same direction, 1 reversed), inherited
    !!         from the base-mesh face where the search crosses a root boundary.
    !! With this, a 2:1 hanging face is exactly "nbr is a leaf with level(nbr) = level(node)-1",
    !! and finer neighbours are exactly "nbr is internal".
    implicit none
    class(QuadTreeMesh2D),intent(in) :: this
    integer,intent(in) :: node,s
    integer,intent(out) :: nbr,ns,nf
    ! Local
    integer :: p,c,pnbr,ps,pf,t,tq

    if(this%level(node) == 0) then
      ! Cross a base-mesh face: neighbour root, its side, and the base flip.
      nbr = this%rootNbr(s,this%rootElem(node))
      ns = this%rootNbrSide(s,this%rootElem(node))
      nf = this%rootFlip(s,this%rootElem(node))
      return
    endif

    p = this%parent(node)
    c = this%quadrant(node)

    if(qt_internal(s,c)) then
      ! Neighbour is the sibling on the other side of an interior face of the parent.
      nbr = this%child(qt_reflect(s,c),p)
      ns = qt_opposite(s)
      nf = 0
      return
    endif

    ! Otherwise ascend: find the parent's neighbour across the same side, then descend.
    call FaceNeighbor_QuadTreeMesh2D(this,p,s,pnbr,ps,pf)
    if(pnbr == 0) then
      nbr = 0; ns = 0; nf = 0
      return
    endif

    if(this%child(1,pnbr) == 0) then
      ! Parent's neighbour is a leaf (equal or larger than the parent) -> our larger neighbour.
      nbr = pnbr; ns = ps; nf = pf
      return
    endif

    ! Parent's neighbour is internal (at level(node)-1): descend one level to the child that
    ! borders `node`, matching sub-positions across the face through the flip.
    t = qt_subpos(c,s)
    if(pf == 0) then
      tq = t
    else
      tq = 3-t
    endif
    nbr = this%child(childOfSide(tq,ps),pnbr)
    ns = ps
    nf = pf

  endsubroutine FaceNeighbor_QuadTreeMesh2D