adtSearch.F90

Go to the documentation of this file.

module adtsearch
    !
    !      Module which contains the subroutines for the global search.
    !
    use constants
    use adtlocalsearch, only: containmenttreesearch, mindistancetreesearch
    use adtutils, only: adtterminate
    use adtdata
    implicit none
 
    !=================================================================
 
contains
 
    !===============================================================
 
    subroutine containmentsearch(nCoor, coor, &
                                 adtID, procID, &
                                 elementType, elementID, &
                                 uvw, nInterpol, &
                                 arrDonor, arrInterpol)
        !
        !        This routine attempts for every coordinate to find the
        !        element in the given ADT, which contains that coordinate.
        !        If no element is found the corresponding entry in procID is
        !        set to -1 to indicate failure.
        !        Subroutine intent(in) arguments.
        !        --------------------------------
        !        nCoor:     Number of coordinates for which the element must
        !                   be determined.
        !        coor:      The coordinates of these points.
        !        adtID:     The ADT to be searched.
        !        nInterpol: Number of variables to be interpolated.
        !        arrDonor:  Array with the donor data; needed to obtain the
        !                   interpolated data.
        !        Subroutine intent(out) arguments.
        !        ---------------------------------
        !        procID:      The ID of the processor in the group of the ADT
        !                     where the element containing the point is
        !                     stored. If no element is found for a given
        !                     point the corresponding entry in procID is set
        !                     to -1 to indicate failure. Remember that the
        !                     processor ID's start at 0 and not at 1.
        !        elementType: The type of element which contains the point.
        !        elementID:   The entry in the connectivity of this element
        !                     which contains the point.
        !        uvw:         The parametric coordinates of the point in the
        !                     transformed element; this transformation is
        !                     such that every element is transformed into a
        !                     standard element in parametric space. The u, v
        !                     and w coordinates can be used to determine the
        !                     actual interpolation weights.
        !        arrInterpol: Array with the interpolated data.
        !
        implicit none
        !
        !       Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: nCoor, nInterpol
        character(len=*), intent(in) :: adtID
 
        real(kind=realtype), dimension(:, :), intent(in) :: coor
        real(kind=realtype), dimension(:, :), intent(in) :: arrdonor
 
        integer, dimension(:), intent(out) :: procID
        integer(kind=intType), dimension(:), intent(out) :: elementID
 
        integer(kind=adtElementType), dimension(:), intent(out) :: &
            elementType
        real(kind=realtype), dimension(:, :), intent(out) :: uvw
        real(kind=realtype), dimension(:, :), intent(out) :: arrinterpol
        !
        !       Local variables.
        !
        integer :: ierr
 
        type(adttype), pointer :: ADT
 
        integer(kind=intType) :: jj, nAlloc
 
        real(kind=realtype), dimension(1) :: dummy
 
        ! Determine the index in the array ADTs, which stores the given
        ! ID. As the number of trees stored is limited, a linear search
        ! algorithm is okay.
 
        nalloc = ubound(adts, 1)
        do jj = 1, nalloc
            if (adtid == adts(jj)%adtID) exit
        end do
 
        ! Check if the ADT to be searched exists. If not stop.
        ! Note that adtTerminate is not called. The reason is that the
        ! processor ID is not known.
 
        if (jj > nalloc) stop "ADT to be searched does not exist."
 
        ! Set pointer to the ADT we will be working with
 
        adt => adts(jj)
 
        ! Check if the ADT corresponds to a volume grid. If not terminate.
 
        if (adt%adtType /= adtvolumeadt) then
            if (adt%myID == 0) &
                call adtterminate(adt, "containmentSearch", &
                                  "ADT does not contain a volume mesh.")
            call mpi_barrier(adt%comm, ierr)
        end if
 
        ! Initialize the search, i.e. determine the number of
        ! coordinates to be searched in each of the local ADT's.
 
        call initsearch(ncoor, coor, dummy, adt, .true.)
 
        ! Perform the actual search.
 
        call search(ncoor, coor, procid, elementtype, &
                    elementid, uvw, dummy, adt, &
                    .true., ninterpol, arrdonor, arrinterpol)
 
    end subroutine containmentsearch
 
    subroutine failsafesearch(nCoor, coor, &
                              adtID, procID, &
                              elementType, elementID, &
                              uvw, dist2, &
                              nInterpol, arrDonor, &
                              arrInterpol)
        !
        !        This routine attempts for every coordinate to find the
        !        element in the given ADT, which contains that coordinate.
        !        If no element is found a minimum distance search is
        !        performed, such that always an interpolation can be
        !        performed. To indicate that the element does not contain the
        !        point the element ID is negated.
        !        Subroutine intent(in) arguments.
        !        --------------------------------
        !        nCoor: Number of coordinates for which the element must be
        !               determined.
        !        coor:  The coordinates of these points.
        !        adtID: The ADT to be searched.
        !        nInterpol: Number of variables to be interpolated.
        !        arrDonor:  Array with the donor data; needed to obtain the
        !                   interpolated data.
        !        Subroutine intent(out) arguments.
        !        ---------------------------------
        !        procID:      The ID of the processor in the group of the ADT
        !                     where the element containing the point is
        !                     stored. If no element is found for a given
        !                     point the corresponding entry in procID is set
        !                     to -1 to indicate failure. Remember that the
        !                     processor ID's start at 0 and not at 1.
        !        elementType: The type of element which contains the point.
        !        elementID:   The entry in the connectivity of this element
        !                     which contains the point. The ID is negative if
        !                     the coordinate is outside the element, i.e. if
        !                     a minimum distance search had to be used.
        !        uvw:         The parametric coordinates of the point in the
        !                     transformed element; this transformation is
        !                     such that every element is transformed into a
        !                     standard element in parametric space. The u, v
        !                     and w coordinates can be used to determine the
        !                     actual interpolation weights.
        !        arrInterpol: Array with the interpolated data.
        !        Subroutine intent(inout) arguments.
        !        -----------------------------------
        !        dist2: Minimum distance squared of the coordinates to the
        !               elements of the ADT. On input it should be
        !               initialized by the calling program, possibly to a
        !               large value. In this way it is possible to handle
        !               periodic problems as efficiently as possible.
        !
        implicit none
        !
        !       Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: nCoor, nInterpol
        character(len=*), intent(in) :: adtID
 
        real(kind=realtype), dimension(:, :), intent(in) :: coor
        real(kind=realtype), dimension(:, :), intent(in) :: arrdonor
 
        integer, dimension(:), intent(out) :: procID
        integer(kind=intType), dimension(:), intent(out) :: elementID
 
        integer(kind=adtElementType), dimension(:), intent(out) :: &
            elementType
 
        real(kind=realtype), dimension(:, :), intent(out) :: uvw
        real(kind=realtype), dimension(:, :), intent(out) :: arrinterpol
 
        real(kind=realtype), dimension(:), intent(inout) :: dist2
        !
        !       Local variables.
        !
        integer :: ierr
        integer, dimension(:), allocatable :: tmpProcID
 
        type(adttype), pointer :: ADT
 
        integer(kind=intType) :: i, j, ii, jj, nAlloc, nFail
        integer(kind=intType), dimension(:), allocatable :: tmpElementID
        integer(kind=intType), dimension(:), allocatable :: coorIDs
 
        integer(kind=adtElementType), dimension(:), allocatable :: &
            tmpElementType
 
        real(kind=realtype), dimension(1) :: dummy
 
        real(kind=realtype), dimension(:), allocatable :: tmpdist2
        real(kind=realtype), dimension(:, :), allocatable :: tmpcoor
        real(kind=realtype), dimension(:, :), allocatable :: tmpuvw
        real(kind=realtype), dimension(:, :), allocatable :: tmparrint
 
        ! Determine the index in the array ADTs, which stores the given
        ! ID. As the number of trees stored is limited, a linear search
        ! algorithm is okay.
 
        nalloc = ubound(adts, 1)
        do jj = 1, nalloc
            if (adtid == adts(jj)%adtID) exit
        end do
 
        ! Check if the ADT to be searched exists. If not stop.
        ! Note that adtTerminate is not called. The reason is that the
        ! processor ID is not known.
 
        if (jj > nalloc) stop "ADT to be searched does not exist."
 
        ! Set pointer to the ADT we will be working with
        adt => adts(jj)
 
        ! Check if the ADT corresponds to a volume grid.
        ! If not terminate.
 
        if (adt%adtType /= adtvolumeadt) then
            if (adt%myID == 0) &
                call adtterminate(adt, "failSafeSearch", &
                                  "ADT does not contain a volume mesh.")
            call mpi_barrier(adt%comm, ierr)
        end if
 
        ! Perform the containment search.
 
        call initsearch(ncoor, coor, dummy, adt, .true.)
 
        call search(ncoor, coor, procid, elementtype, &
                    elementid, uvw, dummy, adt, &
                    .true., ninterpol, arrdonor, arrinterpol)
 
        ! Determine the number of coordinates for which the containment
        ! search failed. Set for the other coordinates the distance to
        ! zero.
 
        nfail = 0
        do i = 1, ncoor
            if (procid(i) == -1) then
                nfail = nfail + 1
            else
                dist2(i) = zero
            end if
        end do
 
        ! Determine the global sum of nFail, which is stored in ii.
 
        call mpi_allreduce(nfail, ii, 1, adflow_integer, mpi_max, &
                           adt%comm, ierr)
 
        ! Return if ii == 0, because the minimum distance search
        ! is not needed.
 
        if (ii == 0) return
 
        ! Allocate the memory for the arrays needed for the minimum
        ! distance search.
 
        allocate (tmpcoor(3, nfail), tmpprocid(nfail), &
                  tmpelementtype(nfail), tmpelementid(nfail), &
                  tmpuvw(3, nfail), tmpdist2(nfail), &
                  coorids(nfail), tmparrint(ninterpol, nfail), &
                  stat=ierr)
        if (ierr /= 0) &
             call adtterminate(adt, "failSafeSearch", &
             "Memory allocation failure for the arrays &
             &for the minimum distance search.")
 
        ! Store the information needed for minimum distance search in
        ! the arrays tmpCoor, tmpDist2 and coorIDs.
 
        ii = 0
        do i = 1, ncoor
            if (procid(i) == -1) then
                ii = ii + 1
 
                tmpcoor(1, ii) = coor(1, i)
                tmpcoor(2, ii) = coor(2, i)
                tmpcoor(3, ii) = coor(3, i)
                tmpdist2(ii) = dist2(i)
                coorids(ii) = i
            end if
        end do
 
        ! Perform the minimum distance search for the coordinates for
        ! which the containment search failed.
 
        call initsearch(nfail, tmpcoor, tmpdist2, adt, .false.)
 
        call search(nfail, tmpcoor, tmpprocid, tmpelementtype, &
                    tmpelementid, tmpuvw, tmpdist2, adt, &
                    .false., ninterpol, arrdonor, tmparrint)
 
        ! Copy for the successful searches the data into the arrays to
        ! be returned by this subroutine. Note that the element IDs are
        ! negated to indicate that the coordinate is outside the element.
 
        do i = 1, nfail
            if (tmpprocid(i) /= -1) then
                ii = coorids(i)
 
                procid(ii) = tmpprocid(i)
                elementtype(ii) = tmpelementtype(i)
                elementid(ii) = -tmpelementid(i)
                uvw(1, ii) = tmpuvw(1, i)
                uvw(2, ii) = tmpuvw(2, i)
                uvw(3, ii) = tmpuvw(3, i)
                dist2(ii) = tmpdist2(i)
 
                do j = 1, ninterpol
                    arrinterpol(j, ii) = tmparrint(j, i)
                end do
            end if
        end do
 
        ! Release the memory used in the minimum distance search.
 
        deallocate (tmpcoor, tmpprocid, tmpelementtype, tmpelementid, &
                    tmpuvw, tmpdist2, coorids, tmparrint, &
                    stat=ierr)
        if (ierr /= 0) &
             call adtterminate(adt, "failSafeSearch", &
             "Deallocation failure for the arrays &
             &for the minimum distance search.")
 
    end subroutine failsafesearch
 
    subroutine initsearch(nCoor, coor, dist2, ADT, containmentSearch)
        !
        !        This routine performs the initialization tasks before the
        !        actual search takes place. It determines the number and the
        !        ID's of the coordinates every local tree may have to
        !        interpolate. From this info this routine also determines the
        !        number of rounds needed in the actual algorithm to avoid a
        !        memory bottleneck.
        !        Subroutine intent(in) arguments.
        !        --------------------------------
        !        nCoor:             Number of local coordinates for which the
        !                           element must be determined.
        !        coor:              The coordinates of these points.
        !        ADT:               ADT type whose ADT must be searched
        !        containmentSearch: Whether or not a containment search must
        !                           be performed. If not a minimum distance
        !                           search algorithm is used, which is more
        !                           expensive.
        !        Subroutine intent(inout) arguments.
        !        -----------------------------------
        !        dist2: Guaranteed minimum distance squared of the
        !               coordinates to the elements of the ADT. This array
        !               should be initialized by the calling routine,
        !               possibly to a large value. It is only used for a
        !               minimum distance search.
        !
        implicit none
        !
        !       Subroutine arguments.
        !
        type(adttype), intent(inout) :: adt
        integer(kind=intType), intent(in) :: ncoor
 
        real(kind=realtype), dimension(:, :), intent(in) :: coor
        real(kind=realtype), dimension(:), intent(inout) :: dist2
 
        logical, intent(in) :: containmentsearch
        !
        !       Local variables.
        !
        integer :: ierr, nrootleaves, comm, nprocs, myid
        integer :: myentryinrootprocs
 
        integer, dimension(:), pointer :: rootleavesprocs
 
        integer(kind=intType) :: i, j, k, mm, nn
 
        integer(kind=intType), dimension(:), allocatable :: ncoorperproc
        integer(kind=intType), dimension(:), allocatable :: ncoorfromproc
 
        real(kind=realtype) :: d1, d2, dx, dy, dz
 
        real(kind=realtype), dimension(:, :, :), pointer :: rootbboxes
 
        ! Set some pointers to make the code more readable.
 
        nrootleaves = adt%nRootLeaves
        myentryinrootprocs = adt%myEntryInRootProcs
        rootleavesprocs => adt%rootLeavesProcs
        rootbboxes => adt%rootBBoxes
 
        comm = adt%comm
        nprocs = adt%nProcs
        myid = adt%myID
 
        ! Determine the global maximum of nCoor. This number will serve
        ! as an upper bound for the number of points to be searched
        ! during a round.
 
        call mpi_allreduce(ncoor, ncoormax, 1, adflow_integer, mpi_max, &
                           comm, ierr)
        ncoormax = max(ncoormax, ncoormaxlowerlimit)
        !
        !        Determine for every root leaf the local coordinates, which
        !        should be searched in the corresponding ADT. The criterion
        !        for a containment search is of course different from a
        !        minimum distance search and therefore a distinction must be
        !        made between the two methods.
        !
        ! Allocate the memory for nCoorPerRootLeaf and mCoorPerRootLeaf.
        ! Initially the latter is a copy of the former, but its data
        ! will change in the outer loop in the iterative algorithm, see
        ! adtSearch. The numbering starts at 0, because these arrays
        ! will be cumulative storage format arrays.
 
        allocate (ncoorperrootleaf(0:nprocs), &
                  mcoorperrootleaf(0:nprocs), stat=ierr)
        if (ierr /= 0) &
             call adtterminate(adt, "initSearch", &
             "Memory allocation failure for &
             &nCoorPerRootLeaf and mCoorPerRootLeaf.")
 
        ! Determine for a minimum distance search the guaranteed minimum
        ! distance squared to each of the root leaves and store the
        ! absolute minimum.
 
        testmindistance: if (.not. containmentsearch) then
 
            do j = 1, nrootleaves
                do i = 1, ncoor
                    d1 = abs(coor(1, i) - rootbboxes(1, 1, j))
                    d2 = abs(coor(1, i) - rootbboxes(1, 2, j))
                    dx = max(d1, d2)
 
                    d1 = abs(coor(2, i) - rootbboxes(2, 1, j))
                    d2 = abs(coor(2, i) - rootbboxes(2, 2, j))
                    dy = max(d1, d2)
 
                    d1 = abs(coor(3, i) - rootbboxes(3, 1, j))
                    d2 = abs(coor(3, i) - rootbboxes(3, 2, j))
                    dz = max(d1, d2)
 
                    d2 = dx * dx + dy * dy + dz * dz
                    dist2(i) = min(dist2(i), d2)
                end do
            end do
 
        end if testmindistance
 
        ! Determine the local number of coordinates to be searched in
        ! each of the local trees, i.e. the array nCoorPerRootLeaf. Note
        ! that the processor storing the ADT stores this number rather
        ! than the root leaf. Furthermore nCoorPerRootLeaf is stored in
        ! cumulative storage format and therefore an outer loop over the
        ! number of root leaves is the most logical thing to do.
 
        ncoorperrootleaf(0) = 0
        mm = 0
 
        loop1rootleaves: do j = 1, nrootleaves
 
            ! Determine the next processor which stores a tree; initialize
            ! in the same loop nCoorPerRootLeaf of the processors in
            ! between. The counter mm contains the entry in
            ! nCoorPerRootLeaf where the number is stored. Remember that
            ! the processor ID's start at 0.
 
            do k = (mm + 1), (rootleavesprocs(j) + 1)
                ncoorperrootleaf(k) = ncoorperrootleaf(mm)
            end do
 
            mm = rootleavesprocs(j) + 1
 
            ! Make a distinction between a containment and a
            ! minimum distance search.
 
            test1containment: if (containmentsearch) then
 
                ! Containment search. Loop over the local number of
                ! coordinates and check if the coordinates are within the
                ! bounding box of the current root leaf. If so, update
                ! the counter.
 
                do i = 1, ncoor
                    if (coor(1, i) >= rootbboxes(1, 1, j) .and. &
                        coor(1, i) <= rootbboxes(1, 2, j) .and. &
                        coor(2, i) >= rootbboxes(2, 1, j) .and. &
                        coor(2, i) <= rootbboxes(2, 2, j) .and. &
                        coor(3, i) >= rootbboxes(3, 1, j) .and. &
                        coor(3, i) <= rootbboxes(3, 2, j)) &
                        ncoorperrootleaf(mm) = ncoorperrootleaf(mm) + 1
                end do
 
                else test1containment
 
                ! Minimum distance search. Loop over the local number of
                ! coordinates and determine the possible minimum distance
                ! squared to the bounding box of the current root leaf.
                ! If less than the currently stored guaranteed minimum
                ! distance squared, update the counter.
 
                do i = 1, ncoor
                    if (coor(1, i) < rootbboxes(1, 1, j)) then
                        dx = coor(1, i) - rootbboxes(1, 1, j)
                    else if (coor(1, i) > rootbboxes(1, 2, j)) then
                        dx = coor(1, i) - rootbboxes(1, 2, j)
                    else
                        dx = zero
                    end if
 
                    if (coor(2, i) < rootbboxes(2, 1, j)) then
                        dy = coor(2, i) - rootbboxes(2, 1, j)
                    else if (coor(2, i) > rootbboxes(2, 2, j)) then
                        dy = coor(2, i) - rootbboxes(2, 2, j)
                    else
                        dy = zero
                    end if
 
                    if (coor(3, i) < rootbboxes(3, 1, j)) then
                        dz = coor(3, i) - rootbboxes(3, 1, j)
                    else if (coor(3, i) > rootbboxes(3, 2, j)) then
                        dz = coor(3, i) - rootbboxes(3, 2, j)
                    else
                        dz = zero
                    end if
 
                    d2 = dx * dx + dy * dy + dz * dz
                    if (d2 < dist2(i)) &
                        ncoorperrootleaf(mm) = ncoorperrootleaf(mm) + 1
 
                end do
 
            end if test1containment
 
        end do loop1rootleaves
 
        ! Fill the rest of the array nCoorPerRootLeaf.
 
        do k = (mm + 1), nprocs
            ncoorperrootleaf(k) = ncoorperrootleaf(mm)
        end do
 
        ! Copy nCoorPerRootLeaf to mCoorPerRootLeaf.
 
        do j = 0, nprocs
            mcoorperrootleaf(j) = ncoorperrootleaf(j)
        end do
 
        ! Allocate the memory for coorPerRootLeaf.
 
        nn = ncoorperrootleaf(nprocs)
        allocate (coorperrootleaf(nn), stat=ierr)
        if (ierr /= 0) &
             call adtterminate(adt, "initSearch", &
             "Memory allocation failure for &
             &coorPerRootLeaf.")
 
        ! Repeat the loop over the number of root leaves, but now store
        ! the ID's of the local coordinates.
 
        nn = 0
        loop2rootleaves: do j = 1, nrootleaves
 
            ! Make a distinction between a containment and a
            ! minimum distance search.
 
            test2containment: if (containmentsearch) then
 
                ! Containment search. Store the nodes which are within the
                ! bounding box of the root leaf.
 
                do i = 1, ncoor
                    if (coor(1, i) >= rootbboxes(1, 1, j) .and. &
                        coor(1, i) <= rootbboxes(1, 2, j) .and. &
                        coor(2, i) >= rootbboxes(2, 1, j) .and. &
                        coor(2, i) <= rootbboxes(2, 2, j) .and. &
                        coor(3, i) >= rootbboxes(3, 1, j) .and. &
                        coor(3, i) <= rootbboxes(3, 2, j)) then
                        nn = nn + 1
                        coorperrootleaf(nn) = i
                    end if
                end do
 
                else test2containment
 
                ! Minimum distance search. Store the nodes which have a
                ! smaller possible minimum distance squared than the
                ! currently  stored value.
 
                do i = 1, ncoor
                    if (coor(1, i) < rootbboxes(1, 1, j)) then
                        dx = coor(1, i) - rootbboxes(1, 1, j)
                    else if (coor(1, i) > rootbboxes(1, 2, j)) then
                        dx = coor(1, i) - rootbboxes(1, 2, j)
                    else
                        dx = zero
                    end if
 
                    if (coor(2, i) < rootbboxes(2, 1, j)) then
                        dy = coor(2, i) - rootbboxes(2, 1, j)
                    else if (coor(2, i) > rootbboxes(2, 2, j)) then
                        dy = coor(2, i) - rootbboxes(2, 2, j)
                    else
                        dy = zero
                    end if
 
                    if (coor(3, i) < rootbboxes(3, 1, j)) then
                        dz = coor(3, i) - rootbboxes(3, 1, j)
                    else if (coor(3, i) > rootbboxes(3, 2, j)) then
                        dz = coor(3, i) - rootbboxes(3, 2, j)
                    else
                        dz = zero
                    end if
 
                    d2 = dx * dx + dy * dy + dz * dz
                    if (d2 < dist2(i)) then
                        nn = nn + 1
                        coorperrootleaf(nn) = i
                    end if
                end do
 
            end if test2containment
 
        end do loop2rootleaves
        !
        !        Determine for every tree the number of coordinates from
        !        other processors it should search and from that information
        !        the number of rounds in the outer loop of the search
        !        algorithm in search.
        !
        ! Allocate the memory for some help arrays.
 
        allocate (procrecv(nprocs - 1), ncoorprocrecv(nprocs - 1), &
                  ncoorperproc(0:nprocs - 1), ncoorfromproc(0:nprocs - 1), &
                  stat=ierr)
        if (ierr /= 0) &
            call adtterminate(adt, "initSearch", &
                              "Memory allocation failure for help arrays.")
 
        ! Determine the number of coordinates I want to be searched in
        ! the trees of other processors. Store the number I have to
        ! search in my own tree and set it to 0 afterwards; the local
        ! interpolations are handled differently.
 
        do i = 0, (nprocs - 1)
            ncoorperproc(i) = ncoorperrootleaf(i + 1) - ncoorperrootleaf(i)
        end do
 
        nlocalinterpol = ncoorperproc(myid)
        ncoorperproc(myid) = 0
 
        ! Communicate these numbers to the other processors.
 
        call mpi_alltoall(ncoorperproc, 1, adflow_integer, &
                          ncoorfromproc, 1, adflow_integer, comm, ierr)
 
        ! Determine the total number of coordinates I have to interpolate
        ! and the processor ID's I will receive coordinates from.
 
        nn = 0
        nprocrecv = 0
 
        do i = 0, (nprocs - 1)
            if (ncoorfromproc(i) > 0) then
                nprocrecv = nprocrecv + 1
                procrecv(nprocrecv) = i
                nn = nn + ncoorfromproc(i)
            end if
        end do
 
        ! Determine the number of rounds in the outer loop of the
        ! interpolation algorithm.
 
        i = real(nn, realtype) / real(ncoormax, realtype)
        if (i * ncoormax < nn) i = i + 1
        i = max(i, 1_inttype)
 
        call mpi_allreduce(i, nrounds, 1, adflow_integer, mpi_max, &
                           comm, ierr)
 
        ! Modify the sequence of procRecv, such that number of
        ! coordinates sent from processors during a certain round
        ! is more balanced.
 
        nn = max(nrootleaves, 1_inttype)
        j = (real(myentryinrootprocs, realtype) / real(nn, realtype)) &
            * nprocrecv
 
        do i = 1, nprocrecv
            j = j + 1
            if (j > nprocrecv) j = 1
 
            nn = procrecv(i)
            procrecv(i) = procrecv(j)
            procrecv(j) = nn
        end do
 
        ! Store the number of coordinates to be received in
        ! nCoorProcRecv.
 
        do i = 1, nprocrecv
            ncoorprocrecv(i) = ncoorfromproc(procrecv(i))
        end do
 
        ! Release the memory of the local help arrays.
 
        deallocate (ncoorperproc, ncoorfromproc, stat=ierr)
        if (ierr /= 0) &
             call adtterminate(adt, "initSearch", &
             "Deallocation failure for nCoorPerProc and &
             &nCoorFromProc")
 
    end subroutine initsearch
 
    subroutine mindistancesearch(nCoor, coor, &
                                 adtID, procID, &
                                 elementType, elementID, &
                                 uvw, dist2, &
                                 nInterpol, arrDonor, &
                                 arrInterpol)
        !
        !        This routine attempts for every coordinate to find the
        !        element in the given ADT which minimizes the distance to
        !        this point.
        !        Subroutine intent(in) arguments.
        !        --------------------------------
        !        nCoor: Number of coordinates for which the element must be
        !               determined.
        !        coor:  The coordinates of these points.
        !        adtID: The ADT to be searched.
        !        nInterpol: Number of variables to be interpolated.
        !        arrDonor:  Array with the donor data; needed to obtain the
        !                   interpolated data.
        !        Subroutine intent(out) arguments.
        !        ---------------------------------
        !        procID:      The ID of the processor in the group of the ADT
        !                     where the element containing the point is
        !                     stored. If no element is found for a given
        !                     point the corresponding entry in procID is set
        !                     to -1 to indicate failure. Remember that the
        !                     processor ID's start at 0 and not at 1.
        !        elementType: The type of element which contains the point.
        !        elementID:   The entry in the connectivity of this element
        !                     which contains the point. The ID is negative if
        !                     the coordinate is outside the element.
        !        uvw:         The parametric coordinates of the point in the
        !                     transformed element; this transformation is
        !                     such that every element is transformed into a
        !                     standard element in parametric space. The u, v
        !                     and w coordinates can be used to determine the
        !                     actual interpolation weights. If the tree
        !                     corresponds to a surface mesh the third entry
        !                     of this array will not be filled.
        !        arrInterpol: Array with the interpolated data.
        !        Subroutine intent(inout) arguments.
        !        -----------------------------------
        !        dist2: Minimum distance squared of the coordinates to the
        !               elements of the ADT. On input it should be
        !               initialized by the calling program, possibly to a
        !               large value. In this way it is possible to handle
        !               periodic problems as efficiently as possible.
        !
        implicit none
        !
        !       Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: ncoor, ninterpol
        character(len=*), intent(in) :: adtid
 
        real(kind=realtype), dimension(:, :), intent(in) :: coor
        real(kind=realtype), dimension(:, :), intent(in) :: arrdonor
 
        integer, dimension(:), intent(out) :: procid
        integer(kind=intType), dimension(:), intent(out) :: elementid
 
        integer(kind=adtElementType), dimension(:), intent(out) :: &
            elementtype
 
        real(kind=realtype), dimension(:, :), intent(out) :: uvw
        real(kind=realtype), dimension(:, :), intent(out) :: arrinterpol
 
        real(kind=realtype), dimension(:), intent(inout) :: dist2
        !
        !       Local variables.
        !
        integer(kind=intType) :: jj, nalloc
 
        type(adttype), pointer :: ADT
 
        ! Determine the index in the array ADTs, which stores the given
        ! ID. As the number of trees stored is limited, a linear search
        ! algorithm is okay.
 
        nalloc = ubound(adts, 1)
        do jj = 1, nalloc
            if (adtid == adts(jj)%adtID) exit
        end do
 
        ! Check if the ADT to be searched exists. If not stop.
        ! Note that adtTerminate is not called. The reason is that the
        ! processor ID is not known.
 
        if (jj > nalloc) stop "ADT to be searched does not exist."
 
        ! Set pointer to the ADT we will be working with
 
        adt => adts(jj)
 
        ! Initialize the search, i.e. determine the number of
        ! coordinates to be searched in each of the local ADT's.
 
        call initsearch(ncoor, coor, dist2, adt, .false.)
 
        ! Perform the actual search.
 
        call search(ncoor, coor, procid, elementtype, &
                    elementid, uvw, dist2, adt, &
                    .false., ninterpol, arrdonor, arrinterpol)
 
        ! Negate the elementID if the coordinate is outside the element,
        ! i.e. if the distance is larger than zero.
 
        do jj = 1, ncoor
            if (dist2(jj) > zero) elementid(jj) = -elementid(jj)
        end do
 
    end subroutine mindistancesearch
 
    subroutine search(nCoor, coor, procID, &
                      elementType, elementID, uvw, &
                      dist2, ADT, containmentSearch, &
                      nInterpol, arrDonor, arrInterpol)
        !
        !        This routine implements the parallel part of the search
        !        algorithm and calls the appropriate local tree searches.
        !        Subroutine intent(in) arguments.
        !        --------------------------------
        !        nCoor:             Number of local coordinates for which the
        !                           element must be determined.
        !        coor:              The coordinates of these points.
        !        ADT:                ADT type whose ADT must be searched
        !        containmentSearch: Whether or not a containment search must
        !                           be performed. If not a minimum distance
        !                           search algorithm is used, which is more
        !                           expensive.
        !        nInterpol:         Number of variables to be interpolated.
        !        arrDonor:          Array with the donor data; needed to
        !                           obtain the interpolated data.
        !        Subroutine intent(inout) arguments.
        !        -----------------------------------
        !        dist2: Minimum distance squared of the coordinates to the
        !               elements of the ADT. On input it contains the
        !               guarenteed distance squared to one of the root
        !               leaves. On output it contains the distance squared to
        !               the nearest element of the global tree. It is only
        !               used for a minimum distance search.
        !        Subroutine intent(out) arguments.
        !        ---------------------------------
        !        procID:      The ID of the processor in the group of the ADT
        !                     where the element containing the point is
        !                     stored.
        !        elementType: The type of element which contains the point or
        !                     minimizes the distance to the point.
        !        elementID:   The entry in the connectivity of this element.
        !        uvw:         The parametric coordinates of (the projection
        !                     of) the point in the transformed element; this
        !                     transformation is such that every element is
        !                     transformed into a standard element in
        !                     parametric space. The u, v and w coordinates
        !                     can be used to determine the actual
        !                     interpolation weights.
        !        arrInterpol: Array with the interpolated data.
        !
        implicit none
        !
        !       Subroutine arguments.
        !
        type(adttype), intent(inout) :: ADT
        integer(kind=intType), intent(in) :: nCoor
        integer(kind=intType), intent(in) :: nInterpol
 
        real(kind=realtype), dimension(:, :), intent(in) :: coor
        real(kind=realtype), dimension(:), intent(inout) :: dist2
 
        real(kind=realtype), dimension(:, :), intent(in) :: arrdonor
 
        integer, dimension(:), intent(out) :: procID
        integer(kind=intType), dimension(:), intent(out) :: elementID
 
        integer(kind=adtElementType), dimension(:), intent(out) :: &
            elementType
 
        real(kind=realtype), dimension(:, :), intent(out) :: uvw
        real(kind=realtype), dimension(:, :), intent(out) :: arrinterpol
 
        logical, intent(in) :: containmentSearch
        !
        !       Local variables.
        !
        integer :: ierr
        integer :: comm, nProcs, myID
        integer :: nProcRecvCur, nProcSendCur, nVarCoor, nVarUVW
        integer :: startProcRecv, procCur, sizeMessage
 
        integer, dimension(mpi_status_size) :: mpiStatus
 
        integer, dimension(:), allocatable :: procSendCur
        integer, dimension(:), allocatable :: sendRequest
        integer, dimension(:, :), allocatable :: sendRecvRequest
 
        integer(kind=intType) :: i, j, k, k1, l, m, ii, mm, nn
        integer(kind=intType) :: nLocalInterpolRound
        integer(kind=intType) :: iStartLocal, iEndLocal, nCoorRecv
 
        integer(kind=intType), dimension(:), allocatable :: nCoorPerProc
        integer(kind=intType), dimension(:), allocatable :: nCoorFromProc
 
        integer(kind=intType), dimension(:, :), allocatable :: intRecv
        integer(kind=intType), dimension(:, :), allocatable :: intBuf
 
        real(kind=realtype), dimension(:, :), allocatable :: coorbuf
        real(kind=realtype), dimension(:, :), allocatable :: coorrecv
        real(kind=realtype), dimension(:, :), allocatable :: uvwrecv
        real(kind=realtype), dimension(:, :), allocatable :: uvwbuf
 
        logical, dimension(:), allocatable :: coorRequested
 
        ! Some abbreviations to make the code more readable.
 
        comm = adt%comm
        nprocs = adt%nProcs
        myid = adt%myID
 
        ! Determine the number of variables stored in the coordinate
        ! buffers. For a minimum distance search also the distance
        ! squared is stored, i.e. 4 variables instead of 3.
 
        if (containmentsearch) then
            nvarcoor = 3
        else
            nvarcoor = 4
        end if
 
        ! And the size of the uvw buffers. These contain nVarCoor plus
        ! the number of variables to be interpolated.
 
        nvaruvw = nvarcoor + max(ninterpol, 0_inttype)
 
        ! Initialize procID to -1, which indicates failure.
 
        do i = 1, ncoor
            procid(i) = -1
        end do
 
        ! Allocate the memory for some help arrays used in the search
        ! algorithm.
 
        nn = ncoorperrootleaf(nprocs)
        allocate (procsendcur(nprocs - 1), sendrequest(nprocs - 1), &
                  ncoorperproc(0:nprocs - 1), ncoorfromproc(0:nprocs - 1), &
                  sendrecvrequest(2, nprocs - 1), coorrequested(nn), &
                  stat=ierr)
        if (ierr /= 0) &
            call adtterminate(adt, "search", &
                              "Memory allocation failure for help arrays.")
 
        ! Initialize coorRequested to .false. This indicates that
        ! the corresponding entry in coorPerRootLeaf has not been
        ! requested for interpolation.
 
        do j = 1, nn
            coorrequested(j) = .false.
        end do
 
        ! Initialize the starting position in the array procRecv to 1.
        ! This variable indicates the starting position in procRecv for
        ! the current round. Also initializes nCoorFromProc to 0.
 
        startprocrecv = 1
 
        do i = 0, (nprocs - 1)
            ncoorfromproc(i) = 0
        end do
 
        ! Determine the number of local interpolations per round and
        ! initialize the iStartLocal and iEndLocal, the start and end
        ! indices for the local interpolation of the current round.
 
        nn = nlocalinterpol / nrounds
        if (nn * nrounds < nlocalinterpol) nn = nn + 1
        nlocalinterpolround = nn
 
        istartlocal = 0
        iendlocal = nlocalinterpolround
        !
        !        Iterative algorithm to determine the elements containing the
        !        coordinates or the elements which minimize the distance. The
        !        algorithm consists of a synchronous outer loop over the
        !        number of times the inner loop should be executed. This
        !        inner loop is asynchronous and performs the actual ADT
        !        search. The outer loop is present to avoid that too much
        !        data is communicated to a single processor at once such that
        !        a memory bottleneck occurs.
        !
        outerloop: do mm = 1, nrounds
 
            ! Determine the processors I want data from in this round
            ! as well as the number of coordinates from these nodes.
 
            nprocrecvcur = 0
            ncoorrecv = 0
 
            do i = startprocrecv, nprocrecv
 
                ! Exit the loop if the maximum number of nodes has been
                ! reached.
 
                if (ncoorrecv == ncoormax) exit
 
                ! Update the number of processors from which I will receive
                ! data during this round.
 
                nprocrecvcur = nprocrecvcur + 1
 
                ! Determine the amount I can receive from this processor.
                ! If the upper limit is exceeded cut it off.
 
                j = ncoorprocrecv(i)
                if ((ncoorrecv + j) > ncoormax) j = ncoormax - ncoorrecv
                ncoorrecv = ncoorrecv + j
 
                ! Store the amount j in the appropriate place of
                ! nCoorFromProc and determine the number of nodes to be
                ! sent to this processor in the next round (possibly 0).
 
                ncoorfromproc(procrecv(i)) = j
                ncoorprocrecv(i) = ncoorprocrecv(i) - j
 
                ! Exit the loop if still some data should be received from
                ! this processor in a next round. The difference between this
                ! exit and the one in the beginning of the do loop is that
                ! here the counter i is not update yet and therefore
                ! startProcRecv will be set correctly for the next round.
 
                if (ncoorprocrecv(i) > 0) exit
            end do
 
            ! Do an all to all communication such that every processor
            ! knows the amount of data it should send to other processors.
 
            call mpi_alltoall(ncoorfromproc, 1, adflow_integer, &
                              ncoorperproc, 1, adflow_integer, comm, ierr)
 
            ! Set the non-zero entries of nCoorFromProc to zero again
            ! for the next round.
 
            nn = min(i, nprocrecv)
            do j = startprocrecv, nn
                ncoorfromproc(procrecv(j)) = 0
            end do
 
            ! Set the starting index for the next round.
 
            startprocrecv = i
 
            ! Determine the number of messages I have to send, the
            ! corresponding processors and the total number of points to
            ! be sent.
 
            nprocsendcur = 0
            nn = 0
 
            do i = 0, (nprocs - 1)
                if (ncoorperproc(i) > 0) then
                    nprocsendcur = nprocsendcur + 1
                    procsendcur(nprocsendcur) = i
                    nn = nn + ncoorperproc(i)
                end if
            end do
 
            ! Allocate the memory for the send buffer of the coordinates
            ! and possibly the minimum distance squared.
 
            allocate (coorbuf(nvarcoor, nn), stat=ierr)
            if (ierr /= 0) &
                call adtterminate(adt, "search", &
                                  "Memory allocation failure for coorBuf.")
 
            ! Send the coordinates to the appropriate processors.
            ! Initialize the counter k in the coordinate buffer to 0.
 
            k = 0
            sendcoorloop: do i = 1, nprocsendcur
 
                ! Store the processor ID and the starting entry in
                ! coorPerRootLeaf a bit easier and initialize k1 to k.
 
                proccur = procsendcur(i)
                nn = mcoorperrootleaf(proccur)
                k1 = k
 
                ! Loop to fill to buffer to this processor.
                ! A distinction is made between a containment search and a
                ! minimum distance search. In the former case a coordinate is
                ! only sent if it has not been interpolated yet; in the
                ! latter case it is sent if the distance is larger than zero.
                ! For a minimum distance search also the current distance
                ! squared is stored.
 
                test1containment: if (containmentsearch) then
 
                    ! Containment search. Store the coordinates of the points
                    ! to be searched in coorBuf. Note that the counter j is
                    ! updated inside the loop. This is done to be able to send
                    ! the maximum number of coordinates possible.
 
                    j = 0
                    do
                        ! Exit the loop if the maximum number or if the last
                        ! coordinate for this processor has been reached.
 
                        if (j == ncoorperproc(proccur) .or. &
                            nn == ncoorperrootleaf(proccur + 1)) exit
 
                        ! Update the counter nn and check if this coordinate
                        ! still needs to be interpolated.
 
                        nn = nn + 1
                        l = coorperrootleaf(nn)
 
                        if (procid(l) == -1) then
 
                            ! Coordinate needs to be interpolated. Update the
                            ! counters j and k and copy the coordinate in the
                            ! send buffer.
 
                            j = j + 1
                            k = k + 1
                            coorbuf(1, k) = coor(1, l)
                            coorbuf(2, k) = coor(2, l)
                            coorbuf(3, k) = coor(3, l)
 
                            ! Set the entry in coorRequested to .true. to indicate
                            ! that a request was sent.
 
                            coorrequested(nn) = .true.
 
                        end if
                    end do
 
                    else test1containment
 
                    ! Minimum distance search. Even if an earlier minimum
                    ! distance was found, it is possible that an even smaller
                    ! distance can be found. Therefore only points with a
                    ! minimum distance squared of zero will not be sent. Both
                    ! the coordinates and the current minimum distance squared
                    ! is sent. Again the counter j is updated inside the loop.
 
                    j = 0
                    do
                        ! Exit the loop if the maximum number or if the last
                        ! coordinate for this processor has been reached.
 
                        if (j == ncoorperproc(proccur) .or. &
                            nn == ncoorperrootleaf(proccur + 1)) exit
 
                        ! Update the counter nn and check if this coordinate still
                        ! needs to be interpolated.
 
                        nn = nn + 1
                        l = coorperrootleaf(nn)
 
                        if (dist2(l) > zero) then
 
                            ! Coordinate needs to be interpolated. Update the
                            ! counters j and k and store the coordinates and the
                            ! distance squared in the send buffer.
 
                            j = j + 1
                            k = k + 1
                            coorbuf(1, k) = coor(1, l)
                            coorbuf(2, k) = coor(2, l)
                            coorbuf(3, k) = coor(3, l)
                            coorbuf(4, k) = dist2(l)
 
                            ! Set the entry in coorRequested to .true. to indicate
                            ! that a request was sent.
 
                            coorrequested(nn) = .true.
 
                        end if
                    end do
 
                end if test1containment
 
                ! Determine the size of the message and send it.
                ! Use nonblocking sends to avoid deadlock.
 
                sizemessage = nvarcoor * (k - k1)
                call mpi_isend(coorbuf(1, k1 + 1), sizemessage, adflow_real, &
                               proccur, proccur, comm, &
                               sendrequest(i), ierr)
 
            end do sendcoorloop
            !
            !          Perform the local interpolations. This is done here to
            !          have an overlap between communication and computation.
            !
            ! Determine the local number to be interpolated in this round.
 
            nn = iendlocal - istartlocal
 
            ! Allocate the memory for the coorRecv, the integers for
            ! storing the processor ID, element type and element ID, and
            ! uvwRecv. This memory is also used for the local
            ! interpolation, which explains the max test for intRecv and
            ! uvwRecv. The coordinates are different, because they will be
            ! released after each received message.
 
            i = max(nn, ncoorrecv)
            allocate (coorrecv(nvarcoor, nn), intrecv(3, i), &
                      uvwrecv(nvaruvw, i), stat=ierr)
            if (ierr /= 0) &
                 call adtterminate(adt, "search", &
                 "Memory allocation failure for &
                 &recv arrays")
 
            ! Initialize the counters i and j, which are used to fill the
            ! buffer coorRecv.
 
            j = 0
            i = mcoorperrootleaf(myid)
 
            ! Make a distinction between a containment and a minimum
            ! distance search.
 
            test2containment: if (containmentsearch) then
 
                ! Containment search. Copy the local coordinates to be
                ! searched in coorRecv. As it is possible that in an earlier
                ! round the coordinate has already been found, only take
                ! coordinates whose element has not been found yet.
 
                do
                    ! Exit the loop if the maximum number or if the last
                    ! coordinate of the local interpolation has been reached.
 
                    if (j == nn .or. i == ncoorperrootleaf(myid + 1)) exit
 
                    ! Update the counter i and check if the corresponding
                    ! coordinate still needs to be interpolated. If so, store it
                    ! in coorRecv and set the corresponding entry in
                    ! coorRequested to .true.
 
                    i = i + 1
                    l = coorperrootleaf(i)
 
                    if (procid(l) == -1) then
                        j = j + 1
 
                        coorrecv(1, j) = coor(1, l)
                        coorrecv(2, j) = coor(2, l)
                        coorrecv(3, j) = coor(3, l)
 
                        coorrequested(i) = .true.
                    end if
                end do
 
                ! Perform the local interpolations. Note that j contains the
                ! actual number of coordinates to be searched.
 
                nn = j
                call containmenttreesearch(adt, coorrecv, intrecv, &
                                           uvwrecv, arrdonor, nn, &
                                           ninterpol)
 
                ! Store the interpolation data at the correct location in
                ! the corresponding arrays.
 
                i = mcoorperrootleaf(myid)
                do j = 1, nn
 
                    ! Determine the coordinate entry, which corresponds to the
                    ! counter j. Remember that in the calls to search routines
                    ! only nodes are given which were not interpolated.
 
                    do
                        i = i + 1
                        if (coorrequested(i)) exit
                    end do
                    l = coorperrootleaf(i)
 
                    ! Copy the data if an actual element was found.
 
                    if (intrecv(1, j) >= 0) then
                        procid(l) = intrecv(1, j)
                        elementtype(l) = intrecv(2, j)
                        elementid(l) = intrecv(3, j)
                        uvw(1, l) = uvwrecv(1, j)
                        uvw(2, l) = uvwrecv(2, j)
                        uvw(3, l) = uvwrecv(3, j)
 
                        do m = 1, ninterpol
                            arrinterpol(m, l) = uvwrecv(m + nvarcoor, j)
                        end do
                    end if
 
                end do
 
                else test2containment
 
                ! Minimum distance search. Store the local coordinates and
                ! distance in coorRecv for the coordinates with non-zero
                ! distance.
 
                do
                    ! Exit the loop if the maximum number or if the last
                    ! coordinate of the local interpolation has been reached.
 
                    if (j == nn .or. i == ncoorperrootleaf(myid + 1)) exit
 
                    ! Update the counter i and check if the corresponding
                    ! coordinate still needs to be interpolated. If so, store
                    ! its coordinates and distance squared in coorRecv and set
                    ! the corresponding entry in coorRequested to .true.
 
                    i = i + 1
                    l = coorperrootleaf(i)
 
                    if (dist2(l) > zero) then
                        j = j + 1
 
                        coorrecv(1, j) = coor(1, l)
                        coorrecv(2, j) = coor(2, l)
                        coorrecv(3, j) = coor(3, l)
                        coorrecv(4, j) = dist2(l)
 
                        coorrequested(i) = .true.
                    end if
                end do
 
                ! Perform the local interpolations. Note that j contains the
                ! actual number of coordinates to be searched.
 
                nn = j
                call mindistancetreesearch(adt, coorrecv, intrecv, &
                                           uvwrecv, arrdonor, nn, &
                                           ninterpol)
 
                ! Store the interpolation data at the correct location in
                ! the corresponding arrays.
 
                i = mcoorperrootleaf(myid)
 
                do j = 1, nn
 
                    ! Determine the next coordinate entry, whose request was
                    ! sent to be interpolated. Remember that it is possible
                    ! that some nodes were not sent to be interpolated (if
                    ! their distance squared is zero already).
 
                    do
                        i = i + 1
                        if (coorrequested(i)) exit
                    end do
                    l = coorperrootleaf(i)
 
                    ! Copy the data if an actual element was found and if
                    ! the corresponding distance squared is less than the
                    ! currently stored value.
 
                    if (intrecv(1, j) >= 0 .and. uvwrecv(4, j) < dist2(l)) then
 
                        procid(l) = intrecv(1, j)
                        elementtype(l) = intrecv(2, j)
                        elementid(l) = intrecv(3, j)
                        uvw(1, l) = uvwrecv(1, j)
                        uvw(2, l) = uvwrecv(2, j)
                        uvw(3, l) = uvwrecv(3, j)
                        dist2(l) = uvwrecv(4, j)
 
                        do m = 1, ninterpol
                            arrinterpol(m, l) = uvwrecv(m + nvarcoor, j)
                        end do
                    end if
 
                end do
 
            end if test2containment
 
            ! The buffer coorRecv is not needed anymore.
            ! Release the memory.
 
            deallocate (coorrecv, stat=ierr)
            if (ierr /= 0) &
                call adtterminate(adt, "search", &
                                  "Deallocation failure for coorRecv.")
 
            ! Set iStartLocal and iEndLocal for the next round.
            ! Also update mCoorPerRootLeaf(myID).
 
            istartlocal = iendlocal
            iendlocal = iendlocal + nlocalinterpolround
            iendlocal = max(iendlocal, nlocalinterpol)
 
            mcoorperrootleaf(myid) = i
            !
            !          Perform the interpolations from the other processors.
            !          Their coordinates are received in an arbitrary sequence
            !          using blocking receives and the interpolated data is sent
            !          back using nonblocking sends.
            !
            ! Loop over the number of messages I will receive. The counter
            ! ii contains the current starting position for the buffers
            ! to be sent back to the requesting processors.
 
            ii = 1
            recvsendloop: do i = 1, nprocrecvcur
 
                ! Block until a message arrives and find the source and size
                ! of the message.
 
                call mpi_probe(mpi_any_source, myid, comm, mpistatus, ierr)
 
                proccur = mpistatus(mpi_source)
                call mpi_get_count(mpistatus, adflow_real, sizemessage, ierr)
 
                ! Check in debug mode that the message is of correct size.
 
                if (debug) then
                    if (sizemessage == mpi_undefined .or. &
                        mod(sizemessage, nvarcoor) /= 0) &
                        call adtterminate(adt, "search", &
                                          "Unexpected size of message")
                end if
 
                ! Allocate the memory for the coordinates to be received and
                ! receive them using a blocking receive; the message has
                ! already arrived.
 
                nn = sizemessage / nvarcoor
                allocate (coorrecv(nvarcoor, nn), stat=ierr)
                if (ierr /= 0) &
                     call adtterminate(adt, "search", &
                     "Memory allocation failure for &
                     &coorRecv.")
 
                call mpi_recv(coorrecv, sizemessage, adflow_real, proccur, &
                              myid, comm, mpistatus, ierr)
 
                ! Search the corresponding elements in the local tree and
                ! release coorRecv afterwards.
 
                if (containmentsearch) then
                    call containmenttreesearch(adt, coorrecv, &
                                               intrecv(:, ii:), uvwrecv(:, ii:), &
                                               arrdonor, nn, &
                                               ninterpol)
                else
                    call mindistancetreesearch(adt, coorrecv, &
                                               intrecv(:, ii:), uvwrecv(:, ii:), &
                                               arrdonor, nn, &
                                               ninterpol)
                end if
 
                deallocate (coorrecv, stat=ierr)
                if (ierr /= 0) &
                    call adtterminate(adt, "search", &
                                      "Deallocation failure for coorRecv.")
 
                ! Send the integer and the floating point information back to
                ! the requesting processor.
 
                sizemessage = 3 * nn
                call mpi_isend(intrecv(1, ii), sizemessage, adflow_integer, &
                               proccur, proccur + 1, comm, &
                               sendrecvrequest(1, i), ierr)
 
                sizemessage = nvaruvw * nn
                call mpi_isend(uvwrecv(1, ii), sizemessage, adflow_real, &
                               proccur, proccur + 2, comm, &
                               sendrecvrequest(2, i), ierr)
 
                ! Update the counter ii for the next message.
 
                ii = ii + nn
 
            end do recvsendloop
 
            ! Complete the nonblocking coordinate sends and release the
            ! memory of coorBuf.
 
            do i = 1, nprocsendcur
                call mpi_waitany(nprocsendcur, sendrequest, sizemessage, &
                                 mpistatus, ierr)
            end do
 
            deallocate (coorbuf, stat=ierr)
            if (ierr /= 0) &
                call adtterminate(adt, "search", &
                                  "Deallocation failure for coorBuf.")
 
            ! Loop over the number of processors to which I sent requests.
            ! Now it is time to receive the information they interpolated.
            ! The sequence of receiving messages is arbitrary.
 
            recvloop: do ii = 1, nprocsendcur
 
                ! Block until an integer message arrives. These messages have
                ! tags of myID+1. Also determine the sending processor and
                ! the size of the message.
 
                call mpi_probe(mpi_any_source, myid + 1, comm, mpistatus, ierr)
 
                proccur = mpistatus(mpi_source)
                call mpi_get_count(mpistatus, adflow_integer, sizemessage, ierr)
 
                ! Check in debug mode that the message is of correct size.
 
                if (debug) then
                    if (sizemessage == mpi_undefined .or. &
                        mod(sizemessage, 3) /= 0) &
                        call adtterminate(adt, "search", &
                                          "Unexpected size of message")
                end if
 
                ! Allocate the memory for the integer and uvw buffers, such
                ! that the messages can be received.
 
                nn = sizemessage / 3
                allocate (intbuf(3, nn), uvwbuf(nvaruvw, nn), stat=ierr)
                if (ierr /= 0) &
                     call adtterminate(adt, "search", &
                     "Memory allocation failure for intBuf &
                     &and uvwBuf.")
 
                call mpi_recv(intbuf, sizemessage, adflow_integer, proccur, &
                              myid + 1, comm, mpistatus, ierr)
 
                sizemessage = nvaruvw * nn
                call mpi_recv(uvwbuf, sizemessage, adflow_real, proccur, &
                              myid + 2, comm, mpistatus, ierr)
 
                ! Store the interpolation data at the correct location in the
                ! corresponding arrays. A distinction must be made between
                ! containment search and minimum distance search, because for
                ! the latter it is possible that a better candidate is
                ! already stored.
 
                i = mcoorperrootleaf(proccur)
 
                test3containment: if (containmentsearch) then
 
                    ! Containment search. Loop over the number of points
                    ! requested on the other processor.
 
                    do j = 1, nn
 
                        ! Determine the next coordinate entry, whose request was
                        ! sent to be interpolated. Remember that only nodes are
                        ! sent which were not interpolated.
 
                        do
                            i = i + 1
                            if (coorrequested(i)) exit
                        end do
 
                        ! Copy the data if an actual element was found.
 
                        if (intbuf(1, j) >= 0) then
                            l = coorperrootleaf(i)
 
                            procid(l) = intbuf(1, j)
                            elementtype(l) = intbuf(2, j)
                            elementid(l) = intbuf(3, j)
                            uvw(1, l) = uvwbuf(1, j)
                            uvw(2, l) = uvwbuf(2, j)
                            uvw(3, l) = uvwbuf(3, j)
 
                            do m = 1, ninterpol
                                arrinterpol(m, l) = uvwbuf(m + nvarcoor, j)
                            end do
                        end if
 
                    end do
 
                    else test3containment
 
                    ! Minimum distance search. Loop over the number of points
                    ! requested on the other processor.
 
                    do j = 1, nn
 
                        ! Determine the next coordinate entry, whose request was
                        ! sent to be interpolated. Remember that it is possible
                        ! that some nodes were not sent to be interpolated (if
                        ! their distance squared is zero already).
 
                        do
                            i = i + 1
                            if (coorrequested(i)) exit
                        end do
 
                        ! Copy the data if an actual element was found and if
                        ! the corresponding distance squared is less than the
                        ! currently stored value.
 
                        l = coorperrootleaf(i)
                        if (intbuf(1, j) >= 0 .and. uvwbuf(4, j) < dist2(l)) then
 
                            procid(l) = intbuf(1, j)
                            elementtype(l) = intbuf(2, j)
                            elementid(l) = intbuf(3, j)
                            uvw(1, l) = uvwbuf(1, j)
                            uvw(2, l) = uvwbuf(2, j)
                            uvw(3, l) = uvwbuf(3, j)
                            dist2(l) = uvwbuf(4, j)
 
                            do m = 1, ninterpol
                                arrinterpol(m, l) = uvwbuf(m + nvarcoor, j)
                            end do
                        end if
 
                    end do
 
                end if test3containment
 
                ! Update mCoorPerRootLeaf(procCur) for the next round and
                ! release the memory of intBuf and uvwBuf.
 
                mcoorperrootleaf(proccur) = i
 
                deallocate (intbuf, uvwbuf, stat=ierr)
                if (ierr /= 0) &
                     call adtterminate(adt, "search", &
                     "Deallocation failure for intBuf &
                     &and uvwBuf.")
            end do recvloop
 
            ! Complete the nonblocking sends of the interpolated data.
 
            !nProcRecvCur = 2*nProcRecvCur
            do i = 1, nprocrecvcur
                call mpi_waitany(nprocrecvcur, sendrecvrequest(1, :), sizemessage, &
                                 mpistatus, ierr)
                call mpi_waitany(nprocrecvcur, sendrecvrequest(2, :), sizemessage, &
                                 mpistatus, ierr)
            end do
 
            ! Release the memory of the buffers used in the nonblocking
            ! sends of the interpolated data.
 
            deallocate (intrecv, uvwrecv, stat=ierr)
            if (ierr /= 0) &
                 call adtterminate(adt, "search", &
                 "Deallocation failure for intRecv and &
                 &uvwRecv")
 
            ! Synchronize the processors, because wild cards have been
            ! used in the communication.
 
            call mpi_barrier(comm, ierr)
 
        end do outerloop
 
        ! Release the memory of the help arrays allocated in this
        ! routine.
 
        deallocate (procsendcur, sendrequest, ncoorperproc, &
                    ncoorfromproc, sendrecvrequest, coorrequested, &
                    stat=ierr)
        if (ierr /= 0) &
            call adtterminate(adt, "search", &
                              "Deallocation failure for help arrays.")
 
        ! Release the memory of the help arrays stored in the module
        ! adtData.
 
        deallocate (procrecv, ncoorprocrecv, ncoorperrootleaf, &
                    mcoorperrootleaf, coorperrootleaf, stat=ierr)
        if (ierr /= 0) &
             call adtterminate(adt, "search", &
             "Deallocation failure for help arrays &
             &stored in the module adtData.")
 
    end subroutine search
 
end module adtsearch