oversetCommUtilites.F90

Go to the documentation of this file.

module oversetcommutilities
 
contains
 
    subroutine getcommpattern(oMat, sendList, nSend, recvList, nRecv)
 
        use constants
        use oversetdata, only: cumdomproc, ndomproc, ndomtotal, csrmatrix
        use blockpointers, only: ndom
        use communication, only: nproc, myid
        use sorting, only: unique
        implicit none
 
        ! Input/output
        type(csrmatrix), intent(in) :: oMat
        integer(kind=intType), intent(out) :: sendList(:, :), recvList(:, :)
        integer(kind=intType), intent(out) :: nSend, nRecv
 
        ! Working:
        integer(kind=intType) :: nn, iDom, nnRow, i, jj, ii, nUniqueProc, iProc
        integer(kind=intType), dimension(:), allocatable :: procsForThisRow, inverse, blkProc
        logical :: added
 
        ! Generic routine to determine what I need to send/recv based on the
        ! data provided in the overlap matrix. The '2' in the send and
        ! receive lists will record the processor and the global 'idom'
        ! index, which is suffient to use for the subsequent communication
        ! structure.
 
        nsend = 0
        nrecv = 0
 
        ! These variables are used to compact the sending of
        ! blocks/fringes. The logic is as follows: A a different
        ! processor may need a block/fringe for more than 1 search. This
        ! is manifested by having two or more entries is the rows (or
        ! columns) I own. It would be inefficient to send the same data
        ! to the same processor more than once, so we "uniquify" the
        ! processors before we send. There is also another salient
        ! reason: If we were to send the same data twice, and the other
        ! processor started using the data, we could get a race condition
        ! as it was modified the received fringes (during a search) while
        ! the same fringes were being overwritten by the receive operation.
 
        allocate (procsforthisrow(ndomtotal), inverse(ndomtotal), blkproc(ndomtotal))
 
        ii = 0
        do iproc = 0, nproc - 1
            do i = 1, ndomproc(iproc)
                ii = ii + 1
                blkproc(ii) = iproc
            end do
        end do
 
        ! Loop over the owned rows of the normal matrix
        do nn = 1, ndom
            idom = cumdomproc(myid) + nn
            nnrow = omat%rowPtr(idom + 1) - omat%rowPtr(idom)
            procsforthisrow(1:nnrow) = omat%assignedProc(omat%rowPtr(idom):omat%rowPtr(idom + 1) - 1)
            call unique(procsforthisrow, nnrow, nuniqueproc, inverse)
 
            do jj = 1, nuniqueproc
                if (procsforthisrow(jj) /= myid) then
                    ! This intersection requires a row quantity from me
                    nsend = nsend + 1
                    sendlist(1, nsend) = procsforthisrow(jj)
                    sendlist(2, nsend) = idom
                end if
            end do
        end do
 
        ! Now we loop back through the whole matrix looking at what I have
        ! to do. If there is a row I don't own, I will have to receive it:
 
        do idom = 1, omat%nRow
            added = .false.
            rowloop: do jj = omat%rowPtr(idom), omat%rowPtr(idom + 1) - 1
 
                ! I have to do this intersection
                if (omat%assignedProc(jj) == myid) then
 
                    ! But I don't know the row entry
                    if (.not. (idom > cumdomproc(myid) .and. idom <= cumdomproc(myid + 1))) then
 
                        ! Need to back out what proc the iDom correponds to:
                        nrecv = nrecv + 1
                        recvlist(1, nrecv) = blkproc(idom)
                        recvlist(2, nrecv) = idom
                        added = .true.
                    end if
                end if
 
                ! Just move on to the next row since we only need to receive it once.
                if (added) then
                    exit rowloop
                end if
            end do rowloop
        end do
 
        deallocate (procsforthisrow, inverse, blkproc)
    end subroutine getcommpattern
 
    subroutine getosurfcommpattern(oMat, oMatT, sendList, nSend, &
                                   recvList, nRecv, rBufSize)
 
        ! This subroutine get the the comm pattern to send the oWall types.
        use constants
        use blockpointers, only: ndom
        use oversetdata, only: ndomtotal, csrmatrix, cumdomproc, ndomproc
        use communication, only: myid, nproc
        use sorting, only: unique
        implicit none
 
        ! Input/output
        type(csrmatrix), intent(in) :: oMat, oMatT
        integer(kind=intType), intent(out) :: sendList(:, :), recvList(:, :)
        integer(kind=intType), intent(out) :: nSend, nRecv
        integer(kind=intType), intent(in) :: rBufSize(nDomTotal)
 
        ! Working:
        integer(kind=intType) :: nn, iDom, jDom, nnRow, nnRowT, i, jj, ii, nUniqueProc, iProc
        integer(kind=intType), dimension(:), allocatable :: blkProc, toRecv
        integer(kind=intType), dimension(:), allocatable :: procsForThisRow, inverse
 
        nsend = 0
        nrecv = 0
 
        allocate (procsforthisrow(2 * ndomtotal), inverse(2 * ndomtotal), blkproc(ndomtotal))
 
        ii = 0
        do iproc = 0, nproc - 1
            do i = 1, ndomproc(iproc)
                ii = ii + 1
                blkproc(ii) = iproc
            end do
        end do
 
        ! Loop over the owned rows of the regular matrix and the rows
        ! transposed matrix (ie columns of the regular matrix)
        do nn = 1, ndom
 
            idom = cumdomproc(myid) + nn
 
            ! Only deal with this block if the rbuffer size for the oWall is
            ! greater than zero. If it is zero, it is empty we don't need to
            ! deal with it.
            if (rbufsize(idom) > 0) then
 
                nnrow = omat%rowPtr(idom + 1) - omat%rowPtr(idom)
                procsforthisrow(1:nnrow) = omat%assignedProc(omat%rowPtr(idom):omat%rowPtr(idom + 1) - 1)
 
                nnrowt = omatt%rowPtr(idom + 1) - omatt%rowPtr(idom)
                procsforthisrow(nnrow + 1:nnrow + nnrowt) = &
                    omatt%assignedProc(omatt%rowPtr(idom):omatt%rowPtr(idom + 1) - 1)
 
                call unique(procsforthisrow, nnrow + nnrowt, nuniqueproc, inverse)
 
                do jj = 1, nuniqueproc
                    if (procsforthisrow(jj) /= myid) then
                        ! This intersection requires a row quantity from me
                        nsend = nsend + 1
                        sendlist(1, nsend) = procsforthisrow(jj)
                        sendlist(2, nsend) = idom
                    end if
                end do
            end if
        end do
 
        ! Now we loop back through the whole matrix looking at what I have
        ! to do. If there is a row or column I don't own, I will have to receive it:
        allocate (torecv(ndomtotal))
        torecv = 0
        do idom = 1, ndomtotal
            do jj = omat%rowPtr(idom), omat%rowPtr(idom + 1) - 1
                jdom = omat%colInd(jj)
                if (omat%assignedProc(jj) == myid) then
                    ! I don't have the row entry:
                    if (.not. (idom > cumdomproc(myid) .and. idom <= cumdomproc(myid + 1))) then
                        torecv(idom) = 1
                    end if
                    ! Don't have the column entry:
                    if (.not. (jdom > cumdomproc(myid) .and. jdom <= cumdomproc(myid + 1))) then
                        torecv(jdom) = 1
                    end if
                end if
            end do
        end do
 
        ! Now loop back through and set my recvList. Only add if the
        ! rBufferSize is larger than zero.
        do idom = 1, ndomtotal
            if (torecv(idom) == 1 .and. rbufsize(idom) > 0) then
                nrecv = nrecv + 1
                recvlist(1, nrecv) = blkproc(idom)
                recvlist(2, nrecv) = idom
            end if
        end do
 
        deallocate (procsforthisrow, inverse, blkproc, torecv)
    end subroutine getosurfcommpattern
 
    subroutine sendoblock(oBlock, iDom, iProc, tagOffset, sendCount)
 
        use constants
        use communication, only: adflow_comm_world, sendrequests
        use oversetdata, only: oversetblock
        use utils, only: echk
        implicit none
 
        ! Input/Output
        type(oversetblock), intent(inout) :: oBlock
        integer(kind=intType), intent(in) :: iProc, iDom, tagOffset
        integer(kind=intType), intent(inout) :: sendCount
 
        ! Working
        integer(kind=intType) :: tag, ierr
 
        tag = tagoffset + idom
        sendcount = sendcount + 1
        call mpi_isend(oblock%rBuffer, size(oblock%rbuffer), adflow_real, &
                       iproc, tag, adflow_comm_world, sendrequests(sendcount), ierr)
        call echk(ierr, __file__, __line__)
 
        sendcount = sendcount + 1
        call mpi_isend(oblock%iBuffer, size(oblock%iBuffer), adflow_integer, &
                       iproc, tag, adflow_comm_world, sendrequests(sendcount), ierr)
        call echk(ierr, __file__, __line__)
 
    end subroutine sendoblock
 
    subroutine sendofringe(oFringe, iDom, iProc, tagOffset, sendCount)
 
        use constants
        use communication, only: adflow_comm_world, sendrequests
        use oversetdata, only: oversetfringe
        use utils, only: echk
        implicit none
 
        ! Input/Output
        type(oversetfringe), intent(inout) :: oFringe
        integer(kind=intType), intent(in) :: iProc, iDom, tagOffset
        integer(kind=intType), intent(inout) :: sendCount
 
        ! Working
        integer(kind=intType) :: tag, ierr
 
        tag = idom + tagoffset
        sendcount = sendcount + 1
        call mpi_isend(ofringe%rBuffer, size(ofringe%rbuffer), adflow_real, &
                       iproc, tag, adflow_comm_world, sendrequests(sendcount), ierr)
        call echk(ierr, __file__, __line__)
 
        sendcount = sendcount + 1
        call mpi_isend(ofringe%iBuffer, size(ofringe%iBuffer), adflow_integer, &
                       iproc, tag, adflow_comm_world, sendrequests(sendcount), ierr)
        call echk(ierr, __file__, __line__)
 
    end subroutine sendofringe
 
    subroutine sendosurf(oWall, iDom, iProc, tagOffset, sendCount)
 
        use constants
        use communication, only: sendrequests, adflow_comm_world
        use oversetdata, only: oversetwall
        use utils, only: echk
        implicit none
 
        ! Input/Output
        type(oversetwall), intent(inout) :: oWall
        integer(kind=intType), intent(in) :: iProc, iDom, tagOffset
        integer(kind=intType), intent(inout) :: sendCount
 
        ! Working
        integer(kind=intType) :: tag, ierr
 
        tag = idom + tagoffset
        sendcount = sendcount + 1
        call mpi_isend(owall%rBuffer, size(owall%rbuffer), adflow_real, &
                       iproc, tag, adflow_comm_world, sendrequests(sendcount), ierr)
        call echk(ierr, __file__, __line__)
 
        sendcount = sendcount + 1
        call mpi_isend(owall%iBuffer, size(owall%iBuffer), adflow_integer, &
                       iproc, tag, adflow_comm_world, sendrequests(sendcount), ierr)
        call echk(ierr, __file__, __line__)
 
    end subroutine sendosurf
 
    subroutine recvoblock(oBlock, iDom, iProc, tagOffset, iSize, rSize, &
                          recvCount, recvInfo)
 
        use constants
        use communication, only: adflow_comm_world, recvrequests
        use oversetdata, only: oversetblock
        use utils, only: echk
        implicit none
 
        ! Input/Output
        type(oversetblock), intent(inout) :: oBlock
        integer(kind=intType), intent(in) :: iDom, iProc, tagOffset, rSize, iSize
        integer(kind=intType), intent(inout) :: recvCount
        integer(kind=intType), intent(inout) :: recvInfo(2, recvCount + 2)
 
        ! Working
        integer(kind=intType) :: tag, ierr
 
        tag = tagoffset + idom
        allocate (oblock%rBuffer(rsize), oblock%iBuffer(isize))
 
        recvcount = recvcount + 1
        call mpi_irecv(oblock%rBuffer, rsize, adflow_real, &
                       iproc, tag, adflow_comm_world, recvrequests(recvcount), ierr)
        call echk(ierr, __file__, __line__)
        recvinfo(:, recvcount) = (/idom, 1/)
 
        recvcount = recvcount + 1
        call mpi_irecv(oblock%iBuffer, isize, adflow_integer, &
                       iproc, tag, adflow_comm_world, recvrequests(recvcount), ierr)
        call echk(ierr, __file__, __line__)
        recvinfo(:, recvcount) = (/idom, 2/)
 
    end subroutine recvoblock
 
    subroutine recvofringe(oFringe, iDom, iProc, tagOffset, iSize, rSize, &
                           recvCount, recvInfo)
 
        use constants
        use communication, only: adflow_comm_world, recvrequests
        use oversetdata, only: oversetfringe
        use utils, only: echk
        implicit none
 
        ! Input/Output
        type(oversetfringe), intent(inout) :: oFringe
        integer(kind=intType), intent(in) :: iDom, iProc, tagOffset, rSize, iSize
        integer(kind=intType), intent(inout) :: recvCount
        integer(kind=intType), intent(inout) :: recvInfo(2, recvCount + 2)
 
        ! Working
        integer(kind=intType) :: tag, ierr
 
        tag = tagoffset + idom
        allocate (ofringe%rBuffer(rsize), ofringe%iBuffer(isize))
 
        recvcount = recvcount + 1
        call mpi_irecv(ofringe%rBuffer, rsize, adflow_real, &
                       iproc, tag, adflow_comm_world, recvrequests(recvcount), ierr)
        call echk(ierr, __file__, __line__)
        recvinfo(:, recvcount) = (/idom, 3/)
 
        recvcount = recvcount + 1
        call mpi_irecv(ofringe%iBuffer, isize, adflow_integer, &
                       iproc, tag, adflow_comm_world, recvrequests(recvcount), ierr)
        call echk(ierr, __file__, __line__)
        recvinfo(:, recvcount) = (/idom, 4/)
 
    end subroutine recvofringe
 
    subroutine recvosurf(oWall, iDom, iProc, tagOffset, iSize, rSize, &
                         recvCount, recvInfo)
 
        use constants
        use communication, only: adflow_comm_world, recvrequests
        use oversetdata, only: oversetwall
        use utils, only: echk
        implicit none
 
        ! Input/Output
        type(oversetwall), intent(inout) :: oWall
        integer(kind=intType), intent(in) :: iDom, iProc, tagOffset, rSize, iSize
        integer(kind=intType), intent(inout) :: recvCount
        integer(kind=intType), intent(inout) :: recvInfo(2, recvCount + 2)
 
        ! Working
        integer(kind=intType) :: tag, ierr
 
        tag = tagoffset + idom
        allocate (owall%rBuffer(rsize), owall%iBuffer(isize))
 
        recvcount = recvcount + 1
        call mpi_irecv(owall%rBuffer, rsize, adflow_real, &
                       iproc, tag, adflow_comm_world, recvrequests(recvcount), ierr)
        call echk(ierr, __file__, __line__)
        recvinfo(:, recvcount) = (/idom, 5/)
 
        recvcount = recvcount + 1
 
        call mpi_irecv(owall%iBuffer, isize, adflow_integer, &
                       iproc, tag, adflow_comm_world, recvrequests(recvcount), ierr)
        call echk(ierr, __file__, __line__)
        recvinfo(:, recvcount) = (/idom, 6/)
 
    end subroutine recvosurf
 
    subroutine getfringereturnsizes(oFringeSendList, oFringeRecvList, &
                                    nOFringeSend, nOfringeRecv, oFringes, &
                                    fringeRecvSizes, cumFringeRecv)
 
        ! For this data exchange we use the exact *reverse* of fringe
        ! communication pattern. This communiation simply determines the
        ! number of fringes that must be returned to the owning process.
 
        use constants
        use communication, only: sendrequests, recvrequests, adflow_comm_world
        use utils, only: echk
        use oversetdata, onlY: oversetfringe
        implicit none
 
        ! Input/output
        type(oversetfringe), dimension(:) :: oFringes
        integer(kind=intType), dimension(:, :) :: oFringeSendList, oFringeRecvList
        integer(kind=intType), dimension(:), allocatable :: cumFringeRecv, fringeRecvSizes
        integer(kind=intType) :: nOFringeSend, nOfringeRecv
        ! Working
        integer(kind=intType) :: sendCount, recvCount
        integer(kind=intType) :: iDom, iProc, jj, ierr, index, i
        integer mpiStatus(MPI_STATUS_SIZE)
 
        ! Post all the fringe iSends
        sendcount = 0
        do jj = 1, nofringerecv
 
            iproc = ofringerecvlist(1, jj)
            idom = ofringerecvlist(2, jj)
            sendcount = sendcount + 1
            call mpi_isend(ofringes(idom)%fringeReturnSize, 1, adflow_integer, &
                           iproc, idom, adflow_comm_world, sendrequests(sendcount), ierr)
            call echk(ierr, __file__, __line__)
        end do
 
        allocate (fringerecvsizes(nofringesend))
 
        ! Non-blocking receives
        recvcount = 0
        do jj = 1, nofringesend
 
            iproc = ofringesendlist(1, jj)
            idom = ofringesendlist(2, jj)
            recvcount = recvcount + 1
 
            call mpi_irecv(fringerecvsizes(jj), 1, adflow_integer, &
                           iproc, idom, adflow_comm_world, recvrequests(recvcount), ierr)
            call echk(ierr, __file__, __line__)
        end do
 
        ! Last thing to do wait for all the sends and receives to finish
        do i = 1, sendcount
            call mpi_waitany(sendcount, sendrequests, index, mpistatus, ierr)
            call echk(ierr, __file__, __line__)
        end do
 
        do i = 1, recvcount
            call mpi_waitany(recvcount, recvrequests, index, mpistatus, ierr)
            call echk(ierr, __file__, __line__)
        end do
 
        ! Compute the cumulative form of the fringeRecvSizes
 
        allocate (cumfringerecv(1:nofringesend + 1))
        cumfringerecv(1) = 1
        do jj = 1, nofringesend ! These are the fringes we *sent*
            ! originally, now are going to receive them
            ! back
            cumfringerecv(jj + 1) = cumfringerecv(jj) + fringerecvsizes(jj)
        end do
 
    end subroutine getfringereturnsizes
 
    !
    !       oversetLoadBalance determine the deistributation of donor and
    !       receiver blocks that will result in approximate even load
    !       balancing. The sparse matrix structrue of the overla is
    !       provided. This computation runs on all processors.
 
    subroutine oversetloadbalance(overlap)
 
        use constants
        use communication, only: nproc
        use oversetdata, only: csrmatrix
        implicit none
 
        ! Input/Output
        type(csrmatrix), intent(inout) :: overlap
 
        ! Working paramters
        integer(kind=intType) :: curRow, jj, jj1, iProc, iRow
        real(kind=realtype) :: evencost, potentialsum, targetcost
        real(Kind=realtype) :: totalsearch, totalbuild
 
        real(kind=realtype), dimension(0:nProc - 1) :: proccosts
        real(kind=realtype), dimension(0:nProc) :: cumproccosts
        real(kind=realtype), dimension(overlap%nRow) :: buildcost
        real(kind=realtype), parameter :: tol = 0.1_realtype
        !  real(kind=realType), parameter :: K=10_realType
        logical, dimension(overlap%nnz) :: blockTaken
        logical :: increment
 
        ! Pointers to make code a litte easier to read
        integer(kind=intType), pointer, dimension(:) :: rowPtr, assignedProc
        real(kind=realtype), pointer, dimension(:) :: data
 
        ! Set the couple  of pointers
        rowptr => overlap%rowPtr
        assignedproc => overlap%assignedProc
        data => overlap%data
 
        ! Determine the total search cost:
        totalsearch = sum(overlap%data)
 
        ! Target amount of work for each processor
        evencost = totalsearch / nproc
 
        ! Initialize the taken processor to False
        blocktaken = .false.
 
        ! Initialzie assignedProc to -1 since there could be entries we can
        ! ignore.
        assignedproc(:) = -1
        proccosts = zero
        cumproccosts(0) = zero
 
        ! Initialize the starting point
        jj = 1
        iproc = 0
 
        ! Find the first row with non-zeros
        currow = 1
        do while (rowptr(currow + 1) - rowptr(currow) == 0)
            currow = currow + 1
        end do
 
        masterloop: do while (currow <= overlap%nRow .and. iproc <= nproc)
 
            ! Normally we increment
            increment = .true.
 
            ! This is our current target cost.
            targetcost = evencost * (iproc + 1)
 
            ! It is still possible that data(jj) is zero. That's ok...we'll
            ! explictly ignore them.
            if (data(jj) /= zero .and. .not. (blocktaken(jj))) then
 
                if (proccosts(iproc) == 0 .or. iproc == nproc - 1) then
                    ! Must be added
                    proccosts(iproc) = proccosts(iproc) + data(jj)
                    blocktaken(jj) = .true.
                    assignedproc(jj) = iproc
 
                else
 
                    ! There is already something in there. See what the
                    !  potential sum will be:
                    potentialsum = cumproccosts(iproc) + proccosts(iproc) + data(jj)
 
                    if (potentialsum < targetcost - tol * evencost) then
                        !  We are not close to our limit yet so just add it normally
                        proccosts(iproc) = proccosts(iproc) + data(jj)
                        blocktaken(jj) = .true.
                        assignedproc(jj) = iproc
 
                    else if (potentialsum >= targetcost - tol * evencost .and. &
                             potentialsum <= targetcost + tol * evencost) then
 
                        ! This one looks perfect. Call it a day...add it and
                        !  move on to the next proc
 
                        proccosts(iproc) = proccosts(iproc) + data(jj)
                        blocktaken(jj) = .true.
                        assignedproc(jj) = iproc
 
                        ! Processor can be incremented
                        cumproccosts(iproc + 1) = cumproccosts(iproc) + proccosts(iproc)
                        iproc = iproc + 1
                    else
                        ! This means potentialSum > targetCost + tol*evenCost
 
                        ! This is somewhat bad news...this may be *horrendly*
                        ! load balanced. The algorithm dictates we *MUST*
                        ! finish this proc no matter what before we go back to
                        ! the outer loop. Essentially we know jj is bad,
                        ! instead scan over the rest of the row and see if we
                        ! can add something else that is decent.
                        increment = .false.
 
                        restofrow: do jj1 = jj + 1, rowptr(currow + 1) - 1
 
                            potentialsum = cumproccosts(iproc) + proccosts(iproc) + data(jj1)
 
                            if (data(jj1) /= zero .and. .not. (blocktaken(jj1))) then
 
                                if (potentialsum < targetcost - tol * evencost) then
                                    !Huh...that one fit in without going
                                    ! over....add it and kep going in the loop
 
                                    proccosts(iproc) = proccosts(iproc) + data(jj1)
                                    blocktaken(jj1) = .true.
                                    assignedproc(jj1) = iproc
 
                                else if (potentialsum >= targetcost - tol * evencost .and. &
                                         potentialsum <= targetcost + tol * evencost) then
 
                                    ! This one fit in perfectly.
                                    proccosts(iproc) = proccosts(iproc) + data(jj1)
                                    blocktaken(jj1) = .true.
                                    assignedproc(jj1) = iproc
 
                                    ! No need to keep going
                                    exit restofrow
 
                                end if
                            end if
                        end do restofrow
 
                        ! Well, the loop finished, we may or may not have
                        ! added something. If so great...if not, oh well. We
                        ! just keep going to the next proc. That's the greedy
                        ! algorithm for you.
 
                        ! Processor can be incremented
                        cumproccosts(iproc + 1) = cumproccosts(iproc) + proccosts(iproc)
                        iproc = iproc + 1
                    end if
                end if
            end if
 
            ! Move 1 in jj, until we reach the end and wrap around.
            if (increment) then
                jj = jj + 1
 
                ! Switch to the next row:
                if (jj == rowptr(currow + 1)) then
 
                    ! This is really tricky...we know we're at the end of the
                    ! row, but we have to SKIP OVER THE EMPTY rows, or else the
                    ! algorithm will crap out. Keep incrementing the curRow
                    ! until we get a row with something in it. Make sure we
                    ! don't go out the end, so check again nRow
 
                    findnextnonzerorow: do while (jj == rowptr(currow + 1))
                        currow = currow + 1
                        if (currow > overlap%nRow) then
                            exit findnextnonzerorow
                        end if
                    end do findnextnonzerorow
                end if
            end if
        end do masterloop
 
    end subroutine oversetloadbalance
 
    subroutine exchangefringes(level, sps, commPattern, internal)
        !
        !       ExchangeFringes exchanges the donorInformation of the fringes:
        !       donorProc, donorBlock, dIndex and donorFrac. It does this
        !       the 1:1 halos for the given level and spectral instance. Since
        !       we have real values and integer values we will do all the ints
        !       first and then the reals.
        !
        use constants
        use block, only: fringetype
        use blockpointers, only: flowdoms
        use communication, only: commtype, internalcommtype, recvbuffer, sendbuffer, myid, &
                                 adflow_comm_world, sendrequests, recvrequests
        use oversetutilities, only: addtofringelist, windindex
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: level, sps
 
        type(commtype), dimension(*), intent(in) :: commPattern
        type(internalcommtype), dimension(*), intent(in) :: internal
        !
        !      Local variables.
        !
        integer :: size, procId, ierr, index
        integer, dimension(mpi_status_size) :: mpiStatus
 
        integer(kind=intType) :: i, j, ii, jj, nVar, iFringe, jFringe
        integer(kind=intType) :: d1, i1, j1, k1, d2, i2, j2, k2
        integer(kind=intType) :: il, jl, kl, myIndex
        type(fringetype) :: fringe
        integer(kind=intType), dimension(:), allocatable :: sendBufInt
        integer(kind=intType), dimension(:), allocatable :: recvBufInt
 
        ! Allocate the memory for the sending and receiving buffers.
        nvar = 3
        ii = commpattern(level)%nProcSend
        ii = commpattern(level)%nsendCum(ii)
        jj = commpattern(level)%nProcRecv
        jj = commpattern(level)%nrecvCum(jj)
 
        allocate (sendbufint(ii * nvar), recvbufint(jj * nvar), stat=ierr)
 
        ! Send the variables. The data is first copied into
        ! the send buffer after which the buffer is sent asap.
 
        ii = 1
        intsends: do i = 1, commpattern(level)%nProcSend
 
            ! Store the processor id and the size of the message
            ! a bit easier.
 
            procid = commpattern(level)%sendProc(i)
            size = nvar * commpattern(level)%nsend(i)
 
            ! Copy the data in the correct part of the send buffer.
 
            jj = ii
            do j = 1, commpattern(level)%nsend(i)
 
                ! Store the block id and the indices of the donor
                ! a bit easier.
 
                d1 = commpattern(level)%sendList(i)%block(j)
                i1 = commpattern(level)%sendList(i)%indices(j, 1)
                j1 = commpattern(level)%sendList(i)%indices(j, 2)
                k1 = commpattern(level)%sendList(i)%indices(j, 3)
 
                ! Copy integer values to buffer
                ifringe = flowdoms(d1, level, sps)%fringePtr(1, i1, j1, k1)
                if (ifringe > 0) then
                    sendbufint(jj) = flowdoms(d1, level, sps)%fringes(ifringe)%donorProc
                    sendbufint(jj + 1) = flowdoms(d1, level, sps)%fringes(ifringe)%donorBlock
                    sendbufint(jj + 2) = flowdoms(d1, level, sps)%fringes(ifringe)%dIndex
                else
                    sendbufint(jj) = -1
                    sendbufint(jj + 1) = 0
                    sendbufint(jj + 2) = 0
                end if
 
                jj = jj + nvar
 
            end do
 
            ! Send the data.
 
            call mpi_isend(sendbufint(ii), size, adflow_integer, procid, &
                           procid, adflow_comm_world, sendrequests(i), &
                           ierr)
 
            ! Set ii to jj for the next processor.
 
            ii = jj
 
        end do intsends
 
        ! Post the nonblocking receives.
 
        ii = 1
        intreceives: do i = 1, commpattern(level)%nProcRecv
 
            ! Store the processor id and the size of the message
            ! a bit easier.
 
            procid = commpattern(level)%recvProc(i)
            size = nvar * commpattern(level)%nrecv(i)
 
            ! Post the receive.
 
            call mpi_irecv(recvbufint(ii), size, adflow_integer, procid, &
                           myid, adflow_comm_world, recvrequests(i), ierr)
 
            ! And update ii.
 
            ii = ii + size
 
        end do intreceives
 
        ! Copy the local data.
 
        intlocalcopy: do i = 1, internal(level)%ncopy
 
            ! Store the block and the indices of the donor a bit easier.
 
            d1 = internal(level)%donorBlock(i)
            i1 = internal(level)%donorIndices(i, 1)
            j1 = internal(level)%donorIndices(i, 2)
            k1 = internal(level)%donorIndices(i, 3)
 
            ! Idem for the halo's.
 
            d2 = internal(level)%haloBlock(i)
            i2 = internal(level)%haloIndices(i, 1)
            j2 = internal(level)%haloIndices(i, 2)
            k2 = internal(level)%haloIndices(i, 3)
 
            ifringe = flowdoms(d1, level, sps)%fringePtr(1, i1, j1, k1)
            if (ifringe > 0) then
                ! The sender has an actual fringe. Nowe check if the
                ! receiver has somewhere already to put it:
 
                jfringe = flowdoms(d2, level, sps)%fringePtr(1, i2, j2, k2)
 
                ! Setup the new fringe:
                fringe%myBlock = d2
 
                il = flowdoms(d2, level, sps)%il
                jl = flowdoms(d2, level, sps)%jl
                kl = flowdoms(d2, level, sps)%kl
                fringe%myIndex = windindex(i2, j2, k2, il, jl, kl)
 
                fringe%donorProc = flowdoms(d1, level, sps)%fringes(ifringe)%donorProc
                fringe%donorBlock = flowdoms(d1, level, sps)%fringes(ifringe)%donorBlock
                fringe%dIndex = flowdoms(d1, level, sps)%fringes(ifringe)%dIndex
                fringe%donorFrac = flowdoms(d1, level, sps)%fringes(ifringe)%donorFrac
 
                if (jfringe > 0) then
                    ! Just copy the fringe into the slot. No need to update
                    ! the pointer since it is already correct.
                    flowdoms(d2, level, sps)%fringes(jfringe) = fringe
                else
 
                    ! There is no slot available yet. Tack the fringe onto
                    ! the end of the d2 fringe list and set the pointers
                    ! accordingly.
 
                    call addtofringelist(flowdoms(d2, level, sps)%fringes, &
                                         flowdoms(d2, level, sps)%nDonors, fringe)
 
                    ! Note that all three values (pointer, start and end) are
                    ! all the same.
                    flowdoms(d2, level, sps)%fringePtr(:, i2, j2, k2) = &
                        flowdoms(d2, level, sps)%nDonors
                end if
            else
                ! The donor isn't a receiver so make sure the halo isn't
                ! either. Just set the fringePtr to 0
 
                flowdoms(d2, level, sps)%fringePtr(1, i2, j2, k2) = 0
            end if
 
        end do intlocalcopy
 
        ! Complete the nonblocking receives in an arbitrary sequence and
        ! copy the variables from the buffer into the halo's.
 
        size = commpattern(level)%nProcRecv
        intcompleterecvs: do i = 1, commpattern(level)%nProcRecv
 
            ! Complete any of the requests.
 
            call mpi_waitany(size, recvrequests, index, mpistatus, ierr)
 
            ! Copy the data just arrived in the halo's.
 
            ii = index
            jj = nvar * commpattern(level)%nrecvCum(ii - 1)
            do j = 1, commpattern(level)%nrecv(ii)
 
                ! Store the block and the indices of the halo a bit easier.
 
                d2 = commpattern(level)%recvList(ii)%block(j)
                i2 = commpattern(level)%recvList(ii)%indices(j, 1)
                j2 = commpattern(level)%recvList(ii)%indices(j, 2)
                k2 = commpattern(level)%recvList(ii)%indices(j, 3)
 
                fringe%myBlock = d2
                ! Recompute my Index:
                il = flowdoms(d2, level, sps)%il
                jl = flowdoms(d2, level, sps)%jl
                kl = flowdoms(d2, level, sps)%kl
                fringe%myIndex = windindex(i2, j2, k2, il, jl, kl)
 
                fringe%donorProc = recvbufint(jj + 1)
                fringe%donorBlock = recvbufint(jj + 2)
                fringe%dIndex = recvbufint(jj + 3)
 
                ifringe = flowdoms(d2, level, sps)%fringePtr(1, i2, j2, k2)
                if (ifringe > 0) then
                    ! We have somehwere to to put the data already:
                    flowdoms(d2, level, sps)%fringes(ifringe) = fringe
                else
                    ! We don't somehwhere to put the fringe to add to the list:
                    call addtofringelist(flowdoms(d2, level, sps)%fringes, &
                                         flowdoms(d2, level, sps)%nDonors, fringe)
 
                    ! Note that all three values (pointer, start and end) are
                    ! all the same.
                    flowdoms(d2, level, sps)%fringePtr(:, i2, j2, k2) = &
                        flowdoms(d2, level, sps)%nDonors
                end if
                jj = jj + nvar
            end do
 
        end do intcompleterecvs
 
        ! Complete the nonblocking sends.
 
        size = commpattern(level)%nProcSend
        do i = 1, commpattern(level)%nProcSend
            call mpi_waitany(size, sendrequests, index, mpistatus, ierr)
        end do
 
        ! Done with the integer memory.
 
        deallocate (sendbufint, recvbufint)
 
        ! Now do the real exchange. We can use the regular real buffers here
        ! since they are large enough
 
        ! ================================================================================
 
        ! Allocate the memory for the sending and receiving buffers.
        nvar = 3
 
        ! Send the variables. The data is first copied into
        ! the send buffer after which the buffer is sent asap.
 
        ii = 1
        sends: do i = 1, commpattern(level)%nProcSend
 
            ! Store the processor id and the size of the message
            ! a bit easier.
 
            procid = commpattern(level)%sendProc(i)
            size = nvar * commpattern(level)%nsend(i)
 
            ! Copy the data in the correct part of the send buffer.
 
            jj = ii
            do j = 1, commpattern(level)%nsend(i)
 
                ! Store the block id and the indices of the donor
                ! a bit easier.
 
                d1 = commpattern(level)%sendList(i)%block(j)
                i1 = commpattern(level)%sendList(i)%indices(j, 1)
                j1 = commpattern(level)%sendList(i)%indices(j, 2)
                k1 = commpattern(level)%sendList(i)%indices(j, 3)
 
                ! Copy real values to buffer
                ifringe = flowdoms(d1, level, sps)%fringePtr(1, i1, j1, k1)
                if (ifringe > 0) then
                    sendbuffer(jj:jj + 2) = flowdoms(d1, level, sps)%fringes(ifringe)%donorFrac
                else
                    sendbuffer(jj:jj + 2) = zero
                end if
                jj = jj + nvar
 
            end do
 
            ! Send the data.
 
            call mpi_isend(sendbuffer(ii), size, adflow_real, procid, &
                           procid, adflow_comm_world, sendrequests(i), &
                           ierr)
 
            ! Set ii to jj for the next processor.
 
            ii = jj
 
        end do sends
 
        ! Post the nonblocking receives.
 
        ii = 1
        receives: do i = 1, commpattern(level)%nProcRecv
 
            ! Store the processor id and the size of the message
            ! a bit easier.
 
            procid = commpattern(level)%recvProc(i)
            size = nvar * commpattern(level)%nrecv(i)
 
            ! Post the receive.
 
            call mpi_irecv(recvbuffer(ii), size, adflow_real, procid, &
                           myid, adflow_comm_world, recvrequests(i), ierr)
 
            ! And update ii.
 
            ii = ii + size
 
        end do receives
 
        ! ***********************************************************
        ! No local copy since we copied the fringes directly and the
        ! donorFrac info is already there.
        ! ***********************************************************
 
        ! Complete the nonblocking receives in an arbitrary sequence and
        ! copy the variables from the buffer into the halo's.
 
        size = commpattern(level)%nProcRecv
        completerecvs: do i = 1, commpattern(level)%nProcRecv
 
            ! Complete any of the requests.
 
            call mpi_waitany(size, recvrequests, index, mpistatus, ierr)
 
            ! Copy the data just arrived in the halo's.
 
            ii = index
            jj = nvar * commpattern(level)%nrecvCum(ii - 1)
            do j = 1, commpattern(level)%nrecv(ii)
 
                ! Store the block and the indices of the halo a bit easier.
 
                d2 = commpattern(level)%recvList(ii)%block(j)
                i2 = commpattern(level)%recvList(ii)%indices(j, 1)
                j2 = commpattern(level)%recvList(ii)%indices(j, 2)
                k2 = commpattern(level)%recvList(ii)%indices(j, 3)
 
                ! Now, there should already be a spot available for the
                ! donorFrac since it was created if necessary in the integer exchange.
                ifringe = flowdoms(d2, level, sps)%fringePtr(1, i2, j2, k2)
                flowdoms(d2, level, sps)%fringes(ifringe)%donorFrac = recvbuffer(jj + 1:jj + 3)
 
                jj = jj + nvar
            end do
 
        end do completerecvs
 
        ! Complete the nonblocking sends.
 
        size = commpattern(level)%nProcSend
        do i = 1, commpattern(level)%nProcSend
            call mpi_waitany(size, sendrequests, index, mpistatus, ierr)
        end do
 
    end subroutine exchangefringes
 
    subroutine exchangestatus(level, sps, commPattern, internal)
        !
        !       ExchangeIsCompute exchanges the isCompute flag for the 1 to 1
        !       connectivity for the given level and sps instance.
        !
        use constants
        use blockpointers, only: ndom, ib, jb, kb, flowdoms
        use communication, only: commtype, internalcommtype
        use utils, only: setpointers
        use haloexchange, only: whalo1to1intgeneric
 
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: level, sps
 
        type(commtype), dimension(*), intent(in) :: commPattern
        type(internalcommtype), dimension(*), intent(in) :: internal
        integer(kind=intType) :: nn
 
        domainloop: do nn = 1, ndom
            flowdoms(nn, level, sps)%intCommVars(1)%var => flowdoms(nn, level, sps)%status(:, :, :)
        end do domainloop
 
        ! Run the generic integer exchange
        call whalo1to1intgeneric(1, level, sps, commpattern, internal)
 
    end subroutine exchangestatus
 
    subroutine exchangestatustranspose(level, sps, commPattern, internal)
 
        ! exchangeStatusTranspose performs the *TRANSPOSE* of the normal
        ! halo exchange. That means it takes information *in the halo cells*
        ! and accumulate it into the *owned cells*. In this particular case,
        ! we are transmitting the isDonor and isWallDonor information from
        ! the halos to the owned cells. The "accumulate" operation will be
        ! an MPI_LOR.  Note that this actually hast he same comm structure
        ! as 'whalo1to1_b'.
 
        use constants
        use blockpointers, only: flowdoms
        use communication, only: commtype, internalcommtype, recvbuffer, sendbuffer, myid, &
                                 adflow_comm_world, sendrequests, recvrequests
        use oversetutilities, only: getstatus, setisdonor, setiswalldonor
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: level, sps
        type(commtype), dimension(*), intent(in) :: commPattern
        type(internalcommtype), dimension(*), intent(in) :: internal
        !
        !      Local variables.
        !
        integer :: size, procID, ierr, index
        integer, dimension(mpi_status_size) :: mpiStatus
 
        integer(kind=intType) :: mm
        integer(kind=intType) :: i, j, k, ii, jj
        integer(kind=intType) :: d1, i1, j1, k1, d2, i2, j2, k2
        integer(kind=intType), dimension(:), allocatable :: sendBuf, recvBuf
        logical :: CisDonor, CisHole, CisCompute, CisFloodSeed, CisFlooded, CisWallDonor, CisReceiver
        logical :: DisDonor, DisHole, DisCompute, DisFloodSeed, DisFlooded, DisWallDonor, DisReceiver
        integer(kind=intType) :: cellStatus, donorStatus
 
        ii = commpattern(level)%nProcSend
        ii = commpattern(level)%nsendCum(ii)
        jj = commpattern(level)%nProcRecv
        jj = commpattern(level)%nrecvCum(jj)
 
        ! We are exchanging 1 piece of information
        allocate (sendbuf(ii), recvbuf(jj), stat=ierr)
 
        ! Gather up the seeds into the *recv* buffer. Note we loop
        ! over nProcRECV here! After the buffer is assembled it is
        ! sent off.
 
        jj = 1
        ii = 1
        recvs: do i = 1, commpattern(level)%nProcRecv
 
            ! Store the processor id and the size of the message
            ! a bit easier.
 
            procid = commpattern(level)%recvProc(i)
            size = commpattern(level)%nrecv(i)
 
            ! Copy the data into the buffer
 
            do j = 1, commpattern(level)%nrecv(i)
 
                ! Store the block and the indices of the halo a bit easier.
 
                d2 = commpattern(level)%recvList(i)%block(j)
                i2 = commpattern(level)%recvList(i)%indices(j, 1)
                j2 = commpattern(level)%recvList(i)%indices(j, 2)
                k2 = commpattern(level)%recvList(i)%indices(j, 3)
 
                recvbuf(jj) = flowdoms(d2, level, sps)%status(i2, j2, k2)
                jj = jj + 1
 
            end do
 
            ! Send the data.
            call mpi_isend(recvbuf(ii), size, adflow_integer, procid, &
                           procid, adflow_comm_world, sendrequests(i), &
                           ierr)
 
            ! Set ii to jj for the next processor.
 
            ii = jj
 
        end do recvs
 
        ! Post the nonblocking receives.
 
        ii = 1
        sends: do i = 1, commpattern(level)%nProcSend
 
            ! Store the processor id and the size of the message
            ! a bit easier.
 
            procid = commpattern(level)%sendProc(i)
            size = commpattern(level)%nsend(i)
 
            ! Post the receive.
 
            call mpi_irecv(sendbuf(ii), size, adflow_integer, procid, &
                           myid, adflow_comm_world, recvrequests(i), ierr)
 
            ! And update ii.
 
            ii = ii + size
 
        end do sends
 
        ! Copy the local data.
 
        localcopy: do i = 1, internal(level)%ncopy
 
            ! Store the block and the indices of the donor a bit easier.
 
            d1 = internal(level)%donorBlock(i)
            i1 = internal(level)%donorIndices(i, 1)
            j1 = internal(level)%donorIndices(i, 2)
            k1 = internal(level)%donorIndices(i, 3)
 
            ! Idem for the halo's.
 
            d2 = internal(level)%haloBlock(i)
            i2 = internal(level)%haloIndices(i, 1)
            j2 = internal(level)%haloIndices(i, 2)
            k2 = internal(level)%haloIndices(i, 3)
 
            ! OR operation. Note we modify the '1' values ie. the 'donors'
            ! which are now receivers because of the transpose operation.
            cellstatus = flowdoms(d1, level, sps)%status(i1, j1, k1)
            call getstatus(cellstatus, cisdonor, cishole, ciscompute, &
                           cisfloodseed, cisflooded, ciswalldonor, cisreceiver)
 
            donorstatus = flowdoms(d2, level, sps)%status(i2, j2, k2)
            call getstatus(donorstatus, disdonor, dishole, discompute, &
                           disfloodseed, disflooded, diswalldonor, disreceiver)
 
            call setisdonor(flowdoms(d1, level, sps)%status(i1, j1, k1), &
                            cisdonor .or. disdonor)
 
            call setiswalldonor(flowdoms(d1, level, sps)%status(i1, j1, k1), &
                                ciswalldonor .or. diswalldonor)
 
        end do localcopy
 
        ! Complete the nonblocking receives in an arbitrary sequence and
        ! copy the variables from the buffer into the halo's.
 
        size = commpattern(level)%nProcSend
        completesends: do i = 1, commpattern(level)%nProcSend
 
            ! Complete any of the requests.
 
            call mpi_waitany(size, recvrequests, index, mpistatus, ierr)
 
            ! Copy the data just arrived in the halo's.
 
            ii = index
 
            jj = commpattern(level)%nsendCum(ii - 1)
 
            do j = 1, commpattern(level)%nsend(ii)
 
                ! Store the block and the indices of the halo a bit easier.
 
                d1 = commpattern(level)%sendList(ii)%block(j)
                i1 = commpattern(level)%sendList(ii)%indices(j, 1)
                j1 = commpattern(level)%sendList(ii)%indices(j, 2)
                k1 = commpattern(level)%sendList(ii)%indices(j, 3)
 
                cellstatus = flowdoms(d1, level, sps)%status(i1, j1, k1)
                call getstatus(cellstatus, cisdonor, cishole, ciscompute, &
                               cisfloodseed, cisflooded, ciswalldonor, cisreceiver)
                jj = jj + 1
                donorstatus = sendbuf(jj)
                call getstatus(donorstatus, disdonor, dishole, discompute, &
                               disfloodseed, disflooded, diswalldonor, disreceiver)
 
                call setisdonor(flowdoms(d1, level, sps)%status(i1, j1, k1), &
                                cisdonor .or. disdonor)
 
                call setiswalldonor(flowdoms(d1, level, sps)%status(i1, j1, k1), &
                                    ciswalldonor .or. diswalldonor)
            end do
        end do completesends
 
        ! Complete the nonblocking sends.
 
        size = commpattern(level)%nProcRecv
        do i = 1, commpattern(level)%nProcRecv
            call mpi_waitany(size, sendrequests, index, mpistatus, ierr)
        end do
 
        deallocate (recvbuf, sendbuf)
 
    end subroutine exchangestatustranspose
 
    subroutine setupfringeglobalind(level, sps)
 
        use constants
        use blockpointers
        use communication
        use utils, only: echk
        implicit none
 
        ! This subroutine is used to record the global index of each of
        ! the donors for overset fringes. It has the same comm structure
        ! as  wOverset and flagInvalidDonors.
 
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: level, sps
 
        !
        !      Local variables.
        !
        integer :: size, procId, ierr, index
        integer, dimension(mpi_status_size) :: mpiStatus
 
        integer(kind=intType) :: nVar
        integer(kind=intType) :: i, j, k, ii, jj, iii, jjj, kkk, iFringe
        integer(kind=intType) :: d1, i1, j1, k1, d2, i2, j2, k2, ind
        integer(kind=intType), dimension(:), allocatable :: sendBufInt
        integer(kind=intType), dimension(:), allocatable :: recvBufInt
        logical :: invalid
        type(commtype), pointer :: commPattern
        type(internalcommtype), pointer :: internal
 
        commpattern => commpatternoverset(level, sps)
        internal => internaloverset(level, sps)
 
        ii = commpattern%nProcSend
        ii = commpattern%nsendCum(ii)
        jj = commpattern%nProcRecv
        jj = commpattern%nrecvCum(jj)
        nvar = 8
        allocate (sendbufint(ii * nvar), recvbufint(jj * nvar), stat=ierr)
 
        ! Send the variables. The data is first copied into
        ! the send buffer after which the buffer is sent asap.
 
        ii = 1
        sends: do i = 1, commpattern%nProcSend
 
            ! Store the processor id and the size of the message
            ! a bit easier.
 
            procid = commpattern%sendProc(i)
            size = nvar * commpattern%nsend(i)
 
            ! Copy the data in the correct part of the send buffer.
 
            jj = ii
            do j = 1, commpattern%nsend(i)
 
                ! Store the block id and the indices of the donor
                ! a bit easier.
 
                d1 = commpattern%sendList(i)%block(j)
                i1 = commpattern%sendList(i)%indices(j, 1)
                j1 = commpattern%sendList(i)%indices(j, 2)
                k1 = commpattern%sendList(i)%indices(j, 3)
 
                ! Loop over the 8 donors:
                do kkk = k1, k1 + 1
                    do jjj = j1, j1 + 1
                        do iii = i1, i1 + 1
                            sendbufint(jj) = flowdoms(d1, level, sps)%globalCell(iii, jjj, kkk)
                            jj = jj + 1
                        end do
                    end do
                end do
            end do
 
            ! Send the data.
 
            call mpi_isend(sendbufint(ii), size, adflow_integer, procid, &
                           procid, adflow_comm_world, sendrequests(i), &
                           ierr)
            call echk(ierr, __file__, __line__)
 
            ! Set ii to jj for the next processor.
 
            ii = jj
 
        end do sends
 
        ! Post the nonblocking receives.
 
        ii = 1
        receives: do i = 1, commpattern%nProcRecv
 
            ! Store the processor id and the size of the message
            ! a bit easier.
 
            procid = commpattern%recvProc(i)
            size = nvar * commpattern%nrecv(i)
 
            ! Post the receive.
 
            call mpi_irecv(recvbufint(ii), size, adflow_integer, procid, &
                           myid, adflow_comm_world, recvrequests(i), ierr)
            call echk(ierr, __file__, __line__)
 
            ! And update ii.
 
            ii = ii + size
 
        end do receives
 
        ! Do the local interpolation.
        localinterp: do i = 1, internal%ncopy
 
            ! Store the block and the indices of the donor a bit easier.
 
            d1 = internal%donorBlock(i)
            i1 = internal%donorIndices(i, 1)
            j1 = internal%donorIndices(i, 2)
            k1 = internal%donorIndices(i, 3)
 
            ! Idem for the halo's.
 
            d2 = internal%haloBlock(i)
            i2 = internal%haloIndices(i, 1)
            j2 = internal%haloIndices(i, 2)
            k2 = internal%haloIndices(i, 3)
 
            ! Loop over the 8 donors:
            ind = 0
            do kkk = k1, k1 + 1
                do jjj = j1, j1 + 1
                    do iii = i1, i1 + 1
                        ind = ind + 1
                        flowdoms(d2, level, sps)%gInd(ind, i2, j2, k2) = &
                            flowdoms(d1, level, sps)%globalCell(iii, jjj, kkk)
                    end do
                end do
            end do
        end do localinterp
 
        ! Complete the nonblocking receives in an arbitrary sequence and
        ! copy the variables from the buffer into the halo's.
 
        size = commpattern%nProcRecv
        completerecvs: do i = 1, commpattern%nProcRecv
 
            ! Complete any of the requests.
 
            call mpi_waitany(size, recvrequests, index, mpistatus, ierr)
            call echk(ierr, __file__, __line__)
 
            ! Copy the data just arrived in the halo's.
 
            ii = index
            jj = nvar * commpattern%nrecvCum(ii - 1)
            do j = 1, commpattern%nrecv(ii)
 
                ! Store the block and the indices of the halo a bit easier.
 
                d2 = commpattern%recvList(ii)%block(j)
                i2 = commpattern%recvList(ii)%indices(j, 1)
                j2 = commpattern%recvList(ii)%indices(j, 2)
                k2 = commpattern%recvList(ii)%indices(j, 3)
 
                ! Just set the 8 values
                do ind = 1, 8
                    flowdoms(d2, level, sps)%gInd(ind, i2, j2, k2) = &
                        recvbufint(jj + ind)
                end do
 
                jj = jj + 8
            end do
        end do completerecvs
 
        ! Complete the nonblocking sends.
 
        size = commpattern%nProcSend
        do i = 1, commpattern%nProcSend
            call mpi_waitany(size, sendrequests, index, mpistatus, ierr)
            call echk(ierr, __file__, __line__)
        end do
        deallocate (sendbufint, recvbufint)
 
    end subroutine setupfringeglobalind
 
    subroutine exchangesurfacedelta(zipperFamList, level, sps, commPattern, internal)
        !
        !       ExchangeSurfaceDelta exchanges surface delta to fill up halo
        !       surface cells from adjacent blocks.
        !
        use constants
        use blockpointers, onlY: ndom, flowdoms, nbocos, bctype, bcfaceid, bcdata, &
                                 ib, il, jb, jl, kb, kl
        use communication, only: commtype, internalcommtype, commpatterncell_1st, &
                                 internalcell_1st
        use utils, only: setpointers
        use haloexchange, only: whalo1to1realgeneric
        use sorting, only: faminlist
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in), dimension(:) :: zipperFamList
        integer(kind=intType), intent(in) :: level, sps
 
        type(commtype), dimension(*), intent(in) :: commPattern
        type(internalcommtype), dimension(*), intent(in) :: internal
 
        ! Local
        integer(kind=intType) :: i, j, k, ii, nn, mm
        real(kind=realtype), dimension(:), allocatable :: psave
        real(kind=realtype), dimension(:, :), pointer :: deltaptr
 
        ! Just cheat by exchangint pressure. saving Pressure, dumping deltaPtr into the pressure,
        ! exchanging that and then restoring the pressure
 
        do nn = 1, ndom
            call setpointers(nn, level, sps)
 
            ! Allocate pointer space for the integer flag communication
            allocate (flowdoms(nn, level, sps)%realCommVars(1)%var(1:ib + 1, 1:jb + 1, 1:kb + 1))
 
            ! Push the surface iblank back to the generic volume variable rVar1
            bocoloop: do mm = 1, nbocos
                faminclude: if (faminlist(bcdata(mm)%famID, zipperfamlist)) then
 
                    select case (bcfaceid(mm))
                    case (imin)
                        deltaptr => flowdoms(nn, level, sps)%realCommVars(1)%var(2 + 1, :, :)
                    case (imax)
                        deltaptr => flowdoms(nn, level, sps)%realCommVars(1)%var(il + 1, :, :)
                    case (jmin)
                        deltaptr => flowdoms(nn, level, sps)%realCommVars(1)%var(:, 2 + 1, :)
                    case (jmax)
                        deltaptr => flowdoms(nn, level, sps)%realCommVars(1)%var(:, jl + 1, :)
                    case (kmin)
                        deltaptr => flowdoms(nn, level, sps)%realCommVars(1)%var(:, :, 2 + 1)
                    case (kmax)
                        deltaptr => flowdoms(nn, level, sps)%realCommVars(1)%var(:, :, kl + 1)
                    end select
 
                    ! NO HALOS!
                    do j = bcdata(mm)%jnBeg + 1, bcdata(mm)%jnEnd
                        do i = bcdata(mm)%inBeg + 1, bcdata(mm)%inEnd
 
                            ! Remember to account for the pointer offset since
                            ! the iblank starts at zero
                            deltaptr(i + 1, j + 1) = bcdata(mm)%delta(i, j)
                        end do
                    end do
                end if faminclude
            end do bocoloop
        end do
 
        ! Exchange the variane
        call whalo1to1realgeneric(1, level, sps, commpatterncell_1st, internalcell_1st)
 
        ! Copy back out
        ii = 0
        do nn = 1, ndom
            call setpointers(nn, level, sps)
 
            ! Extract the surface iblank from the volume.
            bocoloop2: do mm = 1, nbocos
                faminclude2: if (faminlist(bcdata(mm)%famID, zipperfamlist)) then
 
                    select case (bcfaceid(mm))
                    case (imin)
                        deltaptr => flowdoms(nn, level, sps)%realCommVars(1)%var(2 + 1, :, :)
                    case (imax)
                        deltaptr => flowdoms(nn, level, sps)%realCommVars(1)%var(il + 1, :, :)
                    case (jmin)
                        deltaptr => flowdoms(nn, level, sps)%realCommVars(1)%var(:, 2 + 1, :)
                    case (jmax)
                        deltaptr => flowdoms(nn, level, sps)%realCommVars(1)%var(:, jl + 1, :)
                    case (kmin)
                        deltaptr => flowdoms(nn, level, sps)%realCommVars(1)%var(:, :, 2 + 1)
                    case (kmax)
                        deltaptr => flowdoms(nn, level, sps)%realCommVars(1)%var(:, :, jl + 1)
                    end select
 
                    ! INCLUDE THE HALOS!
                    do j = bcdata(mm)%jcBeg, bcdata(mm)%jcEnd
                        do i = bcdata(mm)%icBeg, bcdata(mm)%icEnd
 
                            ! Remember to account for the pointer offset since
                            ! the iblank starts at zero
                            bcdata(mm)%delta(i, j) = deltaptr(i + 1, j + 1)
                        end do
                    end do
                end if faminclude2
            end do bocoloop2
 
            ! Now deallocate this pointer
            deallocate (flowdoms(nn, level, sps)%realCommVars(1)%var)
        end do
    end subroutine exchangesurfacedelta
 
    subroutine exchangesurfaceiblanks(zipperFamList, level, sps, commPattern, internal)
        !
        !       ExchangeIblank exchanges the 1 to 1 internal halo's for the
        !       given level and sps instance.
        !
        use constants
        use blockpointers, only: ndom, ib, jb, kb, iblank, bcdata, &
                                 il, jl, kl, bcfaceid, nbocos, flowdoms, bctype
        use communication, only: commtype, internalcommtype
        use utils, only: setpointers
        use haloexchange, only: whalo1to1intgeneric
        use sorting, only: faminlist
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: level, sps
        integer(kind=intType), intent(in), dimension(:) :: zipperFamList
 
        type(commtype), dimension(*), intent(in) :: commPattern
        type(internalcommtype), dimension(*), intent(in) :: internal
 
        ! Local
        integer(kind=intType) :: i, j, k, ii, nn, mm
        integer(kind=intType), dimension(:), allocatable :: iBlankSave
        integer(kind=intType), dimension(:, :), pointer :: ibp
 
        ! Just cheat by saving iBlank iblank array, resusing itand
        ii = 0
        do nn = 1, ndom
            call setpointers(nn, level, sps)
            ii = ii + (ib + 1) * (jb + 1) * (kb + 1)
        end do
 
        allocate (iblanksave(ii))
        ii = 0
        do nn = 1, ndom
            call setpointers(nn, level, sps)
            do k = 0, kb
                do j = 0, jb
                    do i = 0, ib
                        ii = ii + 1
                        iblanksave(ii) = iblank(i, j, k)
                        ! The following algorithm uses the volume iblank array to communicate surface blocking
                        ! across block boundaries. However, for h-topology meshes, for example at a sharp
                        ! trailing edge this breaks, because the surface on the top an bottom of the shape are
                        ! not topologically adjacent in the block structure. If we leave the existing volume
                        ! blanking in place the correct values are communicated, if we zero it as done originally
                        ! the connection is broken. Therefore the following line is commented out.
                        !iblank(i,j,k) = 0 !commented out to fix issue with h-topology blocks on the zipper
 
                    end do
                end do
            end do
 
            ! Push the surface iblank back to the volume:
            bocoloop: do mm = 1, nbocos
                faminclude: if (faminlist(bcdata(mm)%famID, zipperfamlist)) then
 
                    select case (bcfaceid(mm))
                    case (imin)
                        ibp => iblank(2, :, :)
                    case (imax)
                        ibp => iblank(il, :, :)
                    case (jmin)
                        ibp => iblank(:, 2, :)
                    case (jmax)
                        ibp => iblank(:, jl, :)
                    case (kmin)
                        ibp => iblank(:, :, 2)
                    case (kmax)
                        ibp => iblank(:, :, kl)
                    end select
 
                    ! NO HALOS!
                    do j = bcdata(mm)%jnBeg + 1, bcdata(mm)%jnEnd
                        do i = bcdata(mm)%inBeg + 1, bcdata(mm)%inEnd
 
                            ! Remember to account for the pointer offset since
                            ! the iblank starts at zero
                            ibp(i + 1, j + 1) = bcdata(mm)%iBlank(i, j)
                        end do
                    end do
                end if faminclude
            end do bocoloop
        end do
 
        ! Exchange iblanks
        domainloop: do nn = 1, ndom
            flowdoms(nn, level, sps)%intCommVars(1)%var => &
                flowdoms(nn, level, sps)%iblank(:, :, :)
        end do domainloop
 
        ! Run the generic integer exchange
        call whalo1to1intgeneric(1, level, sps, commpattern, internal)
 
        ii = 0
        do nn = 1, ndom
            call setpointers(nn, level, sps)
 
            ! Extract the surface iblank from the volume.
            bocoloop2: do mm = 1, nbocos
                faminclude2: if (faminlist(bcdata(mm)%famID, zipperfamlist)) then
 
                    select case (bcfaceid(mm))
                    case (imin)
                        ibp => iblank(2, :, :)
                    case (imax)
                        ibp => iblank(il, :, :)
                    case (jmin)
                        ibp => iblank(:, 2, :)
                    case (jmax)
                        ibp => iblank(:, jl, :)
                    case (kmin)
                        ibp => iblank(:, :, 2)
                    case (kmax)
                        ibp => iblank(:, :, kl)
                    end select
 
                    ! INCLUDE THE HALOS!
                    do j = bcdata(mm)%jnBeg, bcdata(mm)%jnEnd + 1
                        do i = bcdata(mm)%inBeg, bcdata(mm)%inEnd + 1
                            ! Remember to account for the pointer offset since
                            ! the iblank starts at zero
                            bcdata(mm)%iBlank(i, j) = ibp(i + 1, j + 1)
                        end do
                    end do
                end if faminclude2
            end do bocoloop2
 
            ! Restore the saved array
            do k = 0, kb
                do j = 0, jb
                    do i = 0, ib
                        ii = ii + 1
                        iblank(i, j, k) = iblanksave(ii)
                    end do
                end do
            end do
        end do
        deallocate (iblanksave)
    end subroutine exchangesurfaceiblanks
 
    subroutine emptyoversetcomm(level, sps)
 
        !  Short cut function to make empty overset comm structure for
        !  problems that do not use overset meshes.
 
        use constants
        use communication, only: commpatternoverset, internaloverset
        implicit none
 
        ! Function
        integer(kind=intType), intent(in) :: level, sps
 
        ! Working
        integer(Kind=intType) :: nn, mm, ierr
 
        commpatternoverset(level, sps)%nProcRecv = 0
        allocate (commpatternoverset(level, sps)%recvProc(0))
        allocate (commpatternoverset(level, sps)%nRecv(0))
        allocate (commpatternoverset(level, sps)%recvList(0))
 
        commpatternoverset(level, sps)%nProcSend = 0
        allocate (commpatternoverset(level, sps)%sendProc(0))
        allocate (commpatternoverset(level, sps)%nSend(0))
        allocate (commpatternoverset(level, sps)%sendList(0))
 
        internaloverset(level, sps)%nCopy = 0
        allocate (internaloverset(level, sps)%donorBlock(0))
        allocate (internaloverset(level, sps)%donorIndices(0, 3))
        allocate (internaloverset(level, sps)%donorInterp(0, 8))
        allocate (internaloverset(level, sps)%haloBlock(0))
        allocate (internaloverset(level, sps)%haloIndices(0, 3))
 
    end subroutine emptyoversetcomm
 
    subroutine updateoversetconnectivity(level, sps)
 
        ! This subroutine updates the overset connectivity for a perturbed
        ! mesh. It does *not* completely redo the connectivity. Rather, a
        ! newton search on the existing donors are performed using the
        ! updated coordinates. This type of update is only applicable if the
        ! entire volume mesh is warped as one like with USMesh in IDWarp. This
        ! actually ends up being a fairly small correction most of the time,
        ! however, desipite looks to the contrary is actually quite fast to
        ! run.
 
        use constants
        use communication
        use blockpointers, only: ndom, il, jl, kl, xseed, flowdoms, x, ib, jb, kb, &
                                 ie, je, ke, fringes, scratch
        use haloexchange, only: whalo1to1realgeneric
        use oversetutilities, only: newtonupdate, fractoweights
        use utils, only: setpointers
 
        implicit none
 
        ! Input
        integer(kind=intType), intent(in) :: level, sps
        type(commtype), pointer :: commPattern
        type(internalcommtype), pointer :: internal
 
        ! Working
        integer(kind=intType) :: nn, ii, jj, ierr, i, j, k, d1, i1, j1, k1, d2, i2, j2, k2
        integer(kind=intType) :: size, procID, index, iii, jjj
        integer, dimension(mpi_status_size) :: mpiStatus
        real(kind=realtype) :: frac(3), frac0(3), xcen(3)
        integer(kind=intType), dimension(8), parameter :: indices = (/1, 2, 4, 3, 5, 6, 8, 7/)
 
        ! Set a tolerance for checking whether fractions are between 0 and 1
        real(kind=realtype) :: fractol = 1e-4
 
        ! Pointers to the overset comms to make it easier to read
        commpattern => commpatternoverset(level, sps)
        internal => internaloverset(level, sps)
 
        ! Step 1: Since we need to update donors for all cells including the
        ! donors for double halos, we must know the new cell center
        ! locations for the all of these receivers. Unfortunately, the
        ! double halos don't have coordinates, so we must first perform a
        ! (forward) block-to-block halo exchange to populate the xSeed
        ! values for all cells, including double halos.
 
        do nn = 1, ndom
            call setpointers(nn, level, sps)
 
            if (.not. associated(flowdoms(nn, level, sps)%xSeed)) then
                allocate (flowdoms(nn, level, sps)%XSeed(0:ib, 0:jb, 0:kb, 3))
            end if
            xseed => flowdoms(nn, level, sps)%xSeed
            xseed = zero
            do k = 2, kl
                do j = 2, jl
                    do i = 2, il
                        xseed(i, j, k, :) = eighth * ( &
                                            x(i - 1, j - 1, k - 1, :) + &
                                            x(i, j - 1, k - 1, :) + &
                                            x(i - 1, j, k - 1, :) + &
                                            x(i, j, k - 1, :) + &
                                            x(i - 1, j - 1, k, :) + &
                                            x(i, j - 1, k, :) + &
                                            x(i - 1, j, k, :) + &
                                            x(i, j, k, :)) !+ fringes(i,j,k)%offset
                    end do
                end do
            end do
        end do
 
        ! Exchange the xSeeds.
        do nn = 1, ndom
            flowdoms(nn, level, sps)%realCommVars(1)%var => flowdoms(nn, level, sps)%xSeed(:, :, :, 1)
            flowdoms(nn, level, sps)%realCommVars(2)%var => flowdoms(nn, level, sps)%xSeed(:, :, :, 2)
            flowdoms(nn, level, sps)%realCommVars(3)%var => flowdoms(nn, level, sps)%xSeed(:, :, :, 3)
        end do
 
        ! Run the (foward) generic halo exchange.
        call whalo1to1realgeneric(3, level, sps, commpatterncell_2nd, internalcell_2nd)
 
        ! Step 2: Next we need to communicate the xSeeds to their donor
        ! procs. This means running the overset exchange in REVERSE (ie from
        ! receiver to donor).  Most of this code will look like
        ! wOverset_b. We will runt he newtonUpdate code (below) on the fly
        ! as we receive the data, which should hide some of the comm time.
 
        ! Gather up the seeds into the *recv* buffer. Note we loop over
        ! nProcRECV here! After the buffer is assembled it is send off.
 
        jj = 1
        ii = 1
        recvs: do i = 1, commpattern%nProcRecv
 
            ! Store the processor id and the size of the message
            ! a bit easier.
 
            procid = commpattern%recvProc(i)
            size = 3 * commpattern%nrecv(i)
 
            ! Copy the data into the buffer
 
            do j = 1, commpattern%nrecv(i)
 
                ! Store the block and the indices to make code a bit easier to read
 
                d2 = commpattern%recvList(i)%block(j)
                i2 = commpattern%recvList(i)%indices(j, 1)
                j2 = commpattern%recvList(i)%indices(j, 2)
                k2 = commpattern%recvList(i)%indices(j, 3)
 
                ! Copy the xSeed
                recvbuffer(jj) = flowdoms(d2, level, sps)%xSeed(i2, j2, k2, 1)
                recvbuffer(jj + 1) = flowdoms(d2, level, sps)%xSeed(i2, j2, k2, 2)
                recvbuffer(jj + 2) = flowdoms(d2, level, sps)%xSeed(i2, j2, k2, 3)
                jj = jj + 3
            end do
 
            ! Send the data.
            call mpi_isend(recvbuffer(ii), size, adflow_real, procid, &
                           procid, adflow_comm_world, recvrequests(i), &
                           ierr)
 
            ! Set ii to jj for the next processor.
 
            ii = jj
 
        end do recvs
 
        ! Post the nonblocking receives.
 
        ii = 1
        sends: do i = 1, commpattern%nProcSend
 
            ! Store the processor id and the size of the message
            ! a bit easier.
 
            procid = commpattern%sendProc(i)
            size = 3 * commpattern%nsend(i)
 
            ! Post the receive.
 
            call mpi_irecv(sendbuffer(ii), size, adflow_real, procid, &
                           myid, adflow_comm_world, sendrequests(i), ierr)
 
            ! And update ii.
 
            ii = ii + size
 
        end do sends
 
        ! Do the local interpolation.
 
        localinterp: do i = 1, internal%ncopy
 
            ! Store the block and the indices of the donor a bit easier.
 
            d1 = internal%donorBlock(i)
            i1 = internal%donorIndices(i, 1)
            j1 = internal%donorIndices(i, 2)
            k1 = internal%donorIndices(i, 3)
 
            ! Idem for the halo's.
 
            d2 = internal%haloBlock(i)
            i2 = internal%haloIndices(i, 1)
            j2 = internal%haloIndices(i, 2)
            k2 = internal%haloIndices(i, 3)
 
            ! xCen is the '2'. This was the receiver, but since this is
            ! reverse, it's now the "input"
            xcen = flowdoms(d2, level, sps)%xSeed(i2, j2, k2, :)
 
            ! Store in the comm structure. We need it for the derivatives.
            internal%xCen(i, :) = xcen
 
            ! Do newton update
            frac0 = (/half, half, half/)
            call newtonupdate(xcen, &
                              flowdoms(d1, level, sps)%x(i1 - 1:i1 + 1, j1 - 1:j1 + 1, k1 - 1:k1 + 1, :), frac0, frac)
 
            ! Check if the fractions are between 0 and 1
            if (maxval(frac) > one + fractol .or. minval(frac) < zero - fractol) then
                print *, "Invalid overset connectivity update. Use 'frozen' or 'full' oversetUpdateMode instead."
                error stop
            end if
 
            ! Set the new weights
            call fractoweights(frac, internal%donorInterp(i, :))
        end do localinterp
 
        ! Complete the nonblocking receives in an arbitrary sequence and
        ! copy the variables from the buffer into the halo's.
 
        size = commpattern%nProcSend
        completesends: do i = 1, commpattern%nProcSend
 
            ! Complete any of the requests.
 
            call mpi_waitany(size, sendrequests, index, mpistatus, ierr)
 
            ii = index
 
            jj = 3 * commpattern%nsendCum(ii - 1)
            do j = 1, commpattern%nsend(ii)
 
                ! Store the block and the indices of the halo a bit easier.
 
                d2 = commpattern%sendList(ii)%block(j)
                i2 = commpattern%sendList(ii)%indices(j, 1)
                j2 = commpattern%sendList(ii)%indices(j, 2)
                k2 = commpattern%sendList(ii)%indices(j, 3)
 
                xcen = sendbuffer(jj + 1:jj + 3)
                jj = jj + 3
 
                ! Store in the comm structure. We need it for derivatives.
                commpattern%sendList(ii)%xCen(j, :) = xcen
 
                ! Compute new fraction
                frac0 = (/half, half, half/)
                call newtonupdate(xcen, &
                                  flowdoms(d2, level, sps)%x(i2 - 1:i2 + 1, j2 - 1:j2 + 1, k2 - 1:k2 + 1, :), &
                                  frac0, frac)
 
                ! Check if the fractions are between zero and one
                if (maxval(frac) > one + fractol .or. minval(frac) < zero - fractol) then
                    print *, "Invalid overset connectivity update. Use 'frozen' or 'full' oversetUpdateMode instead."
                    error stop
                end if
 
                ! Set the new weights
                call fractoweights(frac, commpattern%sendList(ii)%interp(j, :))
            end do
 
        end do completesends
 
        ! Complete the nonblocking sends.
 
        size = commpattern%nProcRecv
        do i = 1, commpattern%nProcRecv
            call mpi_waitany(size, recvrequests, index, mpistatus, ierr)
        end do
 
    end subroutine updateoversetconnectivity
#ifndef USE_COMPLEX
    subroutine updateoversetconnectivity_d(level, sps)
 
        ! Forward mode linearization of updateOversetConnectivity
 
        use constants
        use communication
        use blockpointers, only: ndom, il, jl, kl, xseed, flowdoms, x, ib, jb, kb, &
                                 ie, je, ke, fringes, scratch, flowdomsd, xd
        use haloexchange, only: whalo1to1realgeneric
        use oversetutilities_d, only: newtonupdate_d, fractoweights_d
        use utils, only: setpointers_d
 
        implicit none
 
        ! Input
        integer(kind=intType), intent(in) :: level, sps
        type(commtype), pointer :: commPattern
        type(internalcommtype), pointer :: internal
 
        ! Working
        integer(kind=intType) :: nn, ii, jj, ierr, i, j, k, d1, i1, j1, k1, d2, i2, j2, k2
        integer(kind=intType) :: size, procID, index, iii, jjj
        integer, dimension(mpi_status_size) :: mpiStatus
        real(kind=realtype) :: frac(3), fracd(3), frac0(3), xcen(3), xcend(3), weight(8)
        integer(kind=intType), dimension(8), parameter :: indices = (/1, 2, 4, 3, 5, 6, 8, 7/)
 
        ! Pointers to the overset comms to make it easier to read
        commpattern => commpatternoverset(level, sps)
        internal => internaloverset(level, sps)
 
        ! Step 1: Since we need to update donors for all cells including the
        ! donors for double halos, we must know the new cell center
        ! locations for the all of these receivers. Unfortunately, the
        ! double halos don't have coordinates, so we must first perform a
        ! (forward) block-to-block halo exchange to populate the xSeed
        ! values for all cells, including double halos.
 
        do nn = 1, ndom
            call setpointers_d(nn, level, sps)
 
            if (.not. associated(flowdoms(nn, level, sps)%xSeed)) then
                allocate (flowdoms(nn, level, sps)%XSeed(0:ib, 0:jb, 0:kb, 3))
            end if
            xseed => flowdoms(nn, level, sps)%xSeed
            xseed = zero
            do k = 2, kl
                do j = 2, jl
                    do i = 2, il
                        xseed(i, j, k, :) = eighth * ( &
                                            x(i - 1, j - 1, k - 1, :) + &
                                            x(i, j - 1, k - 1, :) + &
                                            x(i - 1, j, k - 1, :) + &
                                            x(i, j, k - 1, :) + &
                                            x(i - 1, j - 1, k, :) + &
                                            x(i, j - 1, k, :) + &
                                            x(i - 1, j, k, :) + &
                                            x(i, j, k, :)) !+ fringes(i,j,k)%offset
 
                        ! Offset is not active so the xSeed_d just has the x
                        ! part. Just dump the values into scratch so we don't
                        ! acllocate any additional memory.
                        scratch(i, j, k, 1:3) = eighth * ( &
                                                xd(i - 1, j - 1, k - 1, :) + &
                                                xd(i, j - 1, k - 1, :) + &
                                                xd(i - 1, j, k - 1, :) + &
                                                xd(i, j, k - 1, :) + &
                                                xd(i - 1, j - 1, k, :) + &
                                                xd(i, j - 1, k, :) + &
                                                xd(i - 1, j, k, :) + &
                                                xd(i, j, k, :))
                    end do
                end do
            end do
        end do
 
        ! Exchange the xSeeds.
        do nn = 1, ndom
            flowdoms(nn, level, sps)%realCommVars(1)%var => flowdoms(nn, level, sps)%xSeed(:, :, :, 1)
            flowdoms(nn, level, sps)%realCommVars(2)%var => flowdoms(nn, level, sps)%xSeed(:, :, :, 2)
            flowdoms(nn, level, sps)%realCommVars(3)%var => flowdoms(nn, level, sps)%xSeed(:, :, :, 3)
            flowdoms(nn, level, sps)%realCommVars(4)%var => flowdoms(nn, level, sps)%scratch(:, :, :, 1)
            flowdoms(nn, level, sps)%realCommVars(5)%var => flowdoms(nn, level, sps)%scratch(:, :, :, 2)
            flowdoms(nn, level, sps)%realCommVars(6)%var => flowdoms(nn, level, sps)%scratch(:, :, :, 3)
        end do
 
        ! Run the (foward) generic halo exchange.
        call whalo1to1realgeneric(6, level, sps, commpatterncell_2nd, internalcell_2nd)
 
        ! Step 2: Next we need to communicate the xSeeds to their donor
        ! procs. This means running the overset exchange in REVERSE (ie from
        ! receiver to donor).  Most of this code will look like
        ! wOverset_b. We will runt he newtonUpdate code (below) on the fly
        ! as we receive the data, which should hide some of the comm time.
 
        ! Gather up the seeds into the *recv* buffer. Note we loop over
        ! nProcRECV here! After the buffer is assembled it is send off.
 
        jj = 1
        ii = 1
        recvs: do i = 1, commpattern%nProcRecv
 
            ! Store the processor id and the size of the message
            ! a bit easier.
 
            procid = commpattern%recvProc(i)
            size = 6 * commpattern%nrecv(i)
 
            ! Copy the data into the buffer
 
            do j = 1, commpattern%nrecv(i)
 
                ! Store the block and the indices to make code a bit easier to read
 
                d2 = commpattern%recvList(i)%block(j)
                i2 = commpattern%recvList(i)%indices(j, 1)
                j2 = commpattern%recvList(i)%indices(j, 2)
                k2 = commpattern%recvList(i)%indices(j, 3)
 
                ! Copy the xSeed and it's derivative
                recvbuffer(jj) = flowdoms(d2, level, sps)%xSeed(i2, j2, k2, 1)
                recvbuffer(jj + 1) = flowdoms(d2, level, sps)%xSeed(i2, j2, k2, 2)
                recvbuffer(jj + 2) = flowdoms(d2, level, sps)%xSeed(i2, j2, k2, 3)
                recvbuffer(jj + 3) = flowdoms(d2, level, sps)%scratch(i2, j2, k2, 1)
                recvbuffer(jj + 4) = flowdoms(d2, level, sps)%scratch(i2, j2, k2, 2)
                recvbuffer(jj + 5) = flowdoms(d2, level, sps)%scratch(i2, j2, k2, 3)
 
                jj = jj + 6
            end do
 
            ! Send the data.
            call mpi_isend(recvbuffer(ii), size, adflow_real, procid, &
                           procid, adflow_comm_world, recvrequests(i), &
                           ierr)
 
            ! Set ii to jj for the next processor.
 
            ii = jj
 
        end do recvs
 
        ! Post the nonblocking receives.
 
        ii = 1
        sends: do i = 1, commpattern%nProcSend
 
            ! Store the processor id and the size of the message
            ! a bit easier.
 
            procid = commpattern%sendProc(i)
            size = 6 * commpattern%nsend(i)
 
            ! Post the receive.
 
            call mpi_irecv(sendbuffer(ii), size, adflow_real, procid, &
                           myid, adflow_comm_world, sendrequests(i), ierr)
 
            ! And update ii.
 
            ii = ii + size
 
        end do sends
 
        ! Do the local interpolation.
 
        localinterp: do i = 1, internal%ncopy
 
            ! Store the block and the indices of the donor a bit easier.
 
            d1 = internal%donorBlock(i)
            i1 = internal%donorIndices(i, 1)
            j1 = internal%donorIndices(i, 2)
            k1 = internal%donorIndices(i, 3)
 
            ! Idem for the halo's.
 
            d2 = internal%haloBlock(i)
            i2 = internal%haloIndices(i, 1)
            j2 = internal%haloIndices(i, 2)
            k2 = internal%haloIndices(i, 3)
 
            xcen = flowdoms(d2, level, sps)%xSeed(i2, j2, k2, :)
            xcend = flowdoms(d2, level, sps)%scratch(i2, j2, k2, 1:3)
            frac0 = (/half, half, half/)
            call newtonupdate_d(xcen, xcend, &
                                flowdoms(d1, level, sps)%x(i1 - 1:i1 + 1, j1 - 1:j1 + 1, k1 - 1:k1 + 1, :), &
                                flowdomsd(d1, level, sps)%x(i1 - 1:i1 + 1, j1 - 1:j1 + 1, k1 - 1:k1 + 1, :), &
                                frac0, frac, fracd)
 
            ! Set the new weights
            call fractoweights_d(frac, fracd, weight, internal%donorInterpd(i, :))
 
        end do localinterp
 
        ! Complete the nonblocking receives in an arbitrary sequence and
        ! copy the variables from the buffer into the halo's.
 
        size = commpattern%nProcSend
        completesends: do i = 1, commpattern%nProcSend
 
            ! Complete any of the requests.
 
            call mpi_waitany(size, sendrequests, index, mpistatus, ierr)
 
            ii = index
 
            jj = 6 * commpattern%nsendCum(ii - 1)
            do j = 1, commpattern%nsend(ii)
 
                ! Store the block and the indices of the halo a bit easier.
 
                d2 = commpattern%sendList(ii)%block(j)
                i2 = commpattern%sendList(ii)%indices(j, 1)
                j2 = commpattern%sendList(ii)%indices(j, 2)
                k2 = commpattern%sendList(ii)%indices(j, 3)
 
                xcen = sendbuffer(jj + 1:jj + 3)
                xcend = sendbuffer(jj + 4:jj + 6)
                jj = jj + 6
 
                ! Compute new fraction
                frac0 = (/half, half, half/)
                call newtonupdate_d(xcen, xcend, &
                                    flowdoms(d2, level, sps)%x(i2 - 1:i2 + 1, j2 - 1:j2 + 1, k2 - 1:k2 + 1, :), &
                                    flowdomsd(d2, level, sps)%x(i2 - 1:i2 + 1, j2 - 1:j2 + 1, k2 - 1:k2 + 1, :), &
                                    frac0, frac, fracd)
 
                ! Set the new weights
                call fractoweights_d(frac, fracd, weight, &
                                     commpattern%sendList(ii)%interpd(j, :))
            end do
 
        end do completesends
 
        ! Complete the nonblocking sends.
 
        size = commpattern%nProcRecv
        do i = 1, commpattern%nProcRecv
            call mpi_waitany(size, recvrequests, index, mpistatus, ierr)
        end do
 
    end subroutine updateoversetconnectivity_d
 
    subroutine updateoversetconnectivity_b(level, sps)
 
        ! Reverse  mode linearization of updateOversetConnectivity
 
        use constants
        use communication
        use blockpointers, only: ndom, il, jl, kl, xseed, flowdoms, x, ib, jb, kb, &
                                 ie, je, ke, fringes, scratch, flowdomsd, xd
        use haloexchange, only: whalo1to1realgeneric_b
        use oversetutilities_b, only: newtonupdate_b, fractoweights_b, newtonupdate, fractoweights
        use utils, only: setpointers_b
 
        implicit none
 
        ! Input
        integer(kind=intType), intent(in) :: level, sps
        type(commtype), pointer :: commPattern
        type(internalcommtype), pointer :: internal
 
        ! Working
        integer(kind=intType) :: nn, ii, jj, kk, ierr, i, j, k, d1, i1, j1, k1, d2, i2, j2, k2
        integer(kind=intType) :: size, procID, index, iii, jjj
        integer, dimension(mpi_status_size) :: mpiStatus
        real(kind=realtype) :: frac(3), fracd(3), frac0(3), xcen(3), xcend(3), weight(8), add(3)
        integer(kind=intType), dimension(8), parameter :: indices = (/1, 2, 4, 3, 5, 6, 8, 7/)
 
        ! Pointers to the overset comms to make it easier to read
        commpattern => commpatternoverset(level, sps)
        internal => internaloverset(level, sps)
 
        ! Zero out xSeedd (in scratch)
        do nn = 1, ndom
            flowdoms(nn, 1, sps)%scratch(:, :, :, 1:3) = zero
        end do
 
        ! In reverse the fist thing we must do is the compute the
        ! sensitivites of xCen and send it to the receiving
        ! processor. This comm pattern is the same as wOverset forward.
 
        ii = 1
        sends: do i = 1, commpattern%nProcSend
 
            ! Store the processor id and the size of the message
            ! a bit easier.
 
            procid = commpattern%sendProc(i)
            size = 3 * commpattern%nsend(i)
 
            ! Copy the data in the correct part of the send buffer.
 
            jj = ii
            do j = 1, commpattern%nsend(i)
 
                ! Store the block id and the indices of the donor
                ! a bit easier.
                d1 = commpattern%sendList(i)%block(j)
                i1 = commpattern%sendList(i)%indices(j, 1)
                j1 = commpattern%sendList(i)%indices(j, 2)
                k1 = commpattern%sendList(i)%indices(j, 3)
 
                ! -------- Recompute forward pass -------------
                xcen = commpattern%sendList(i)%xCen(j, :)
 
                ! Do newton update
                frac0 = (/half, half, half/)
                call newtonupdate(xcen, &
                                  flowdoms(d1, level, sps)%x(i1 - 1:i1 + 1, j1 - 1:j1 + 1, k1 - 1:k1 + 1, :), &
                                  frac0, frac)
 
                ! Set the new weights
                call fractoweights(frac, weight)
 
                ! ------------- Reverse pass -----------
 
                ! Transfer the weights back to the frac
                call fractoweights_b(frac, fracd, weight, commpattern%sendList(i)%interpd(j, :))
 
                ! Run the reverse newton update.  Note that we are
                ! accumulating into the local xd here in the newton_b call.
                frac0 = (/half, half, half/)
                call newtonupdate_b(xcen, xcend, &
                                    flowdoms(d1, level, sps)%x(i1 - 1:i1 + 1, j1 - 1:j1 + 1, k1 - 1:k1 + 1, :), &
                                    flowdomsd(d1, level, sps)%x(i1 - 1:i1 + 1, j1 - 1:j1 + 1, k1 - 1:k1 + 1, :), &
                                    frac0, frac, fracd)
 
                ! -------------------------------------
 
                ! We want to send xCend to the receiver
                sendbuffer(jj:jj + 2) = xcend
                jj = jj + 3
 
            end do
            ! Send the data.
 
            call mpi_isend(sendbuffer(ii), size, adflow_real, procid, &
                           procid, adflow_comm_world, sendrequests(i), &
                           ierr)
 
            ! Set ii to jj for the next processor.
 
            ii = jj
 
        end do sends
 
        ! Post the nonblocking receives.
 
        ii = 1
        receives: do i = 1, commpattern%nProcRecv
 
            ! Store the processor id and the size of the message
            ! a bit easier.
 
            procid = commpattern%recvProc(i)
            size = 3 * commpattern%nrecv(i)
 
            ! Post the receive.
 
            call mpi_irecv(recvbuffer(ii), size, adflow_real, procid, &
                           myid, adflow_comm_world, recvrequests(i), ierr)
 
            ! And update ii.
 
            ii = ii + size
 
        end do receives
 
        ! Do the local stuff while we're waiting
 
        localinterp: do i = 1, internal%ncopy
 
            ! Store the block and the indices of the donor a bit easier.
            d1 = internal%donorBlock(i)
            i1 = internal%donorIndices(i, 1)
            j1 = internal%donorIndices(i, 2)
            k1 = internal%donorIndices(i, 3)
 
            d2 = internal%haloBlock(i)
            i2 = internal%haloIndices(i, 1)
            j2 = internal%haloIndices(i, 2)
            k2 = internal%haloIndices(i, 3)
 
            ! -------- Recompute forward pass -------------
            xcen = internal%XCen(i, :)
 
            ! Do newton update
            frac0 = (/half, half, half/)
            call newtonupdate(xcen, &
                              flowdoms(d1, level, sps)%x(i1 - 1:i1 + 1, j1 - 1:j1 + 1, k1 - 1:k1 + 1, :), frac0, frac)
 
            ! Set the new weights
            call fractoweights(frac, weight)
 
            ! ------------- Reverse pass -----------
 
            ! Transfer the weights back to the frac
            call fractoweights_b(frac, fracd, weight, internal%Donorinterpd(i, :))
 
            ! Run the reverse newton update.  Note that we are
            ! accumulating into the local xd here in the newton_b call.
            frac0 = (/half, half, half/)
            call newtonupdate_b(xcen, xcend, &
                                flowdoms(d1, level, sps)%x(i1 - 1:i1 + 1, j1 - 1:j1 + 1, k1 - 1:k1 + 1, :), &
                                flowdomsd(d1, level, sps)%x(i1 - 1:i1 + 1, j1 - 1:j1 + 1, k1 - 1:k1 + 1, :), &
                                frac0, frac, fracd)
 
            ! Accumulate in to the xSeedd (which we're using scratch for)
            flowdoms(d2, level, sps)%scratch(i2, j2, k2, 1:3) = &
                flowdoms(d2, level, sps)%scratch(i2, j2, k2, 1:3) + xcend
 
        end do localinterp
 
        ! Complete the nonblocking receives in an arbitrary sequence and
        ! copy the variables from the buffer into the halo's.
 
        size = commpattern%nProcRecv
        completerecvs: do i = 1, commpattern%nProcRecv
 
            ! Complete any of the requests.
 
            call mpi_waitany(size, recvrequests, index, mpistatus, ierr)
 
            ! Copy the data just arrived in the halo's.
 
            ii = index
            jj = 3 * commpattern%nrecvCum(ii - 1)
            do j = 1, commpattern%nrecv(ii)
 
                ! Store the block and the indices of the halo a bit easier.
 
                d2 = commpattern%recvList(ii)%block(j)
                i2 = commpattern%recvList(ii)%indices(j, 1)
                j2 = commpattern%recvList(ii)%indices(j, 2)
                k2 = commpattern%recvList(ii)%indices(j, 3)
 
                flowdoms(d2, level, sps)%scratch(i2, j2, k2, 1:3) = &
                    flowdoms(d2, level, sps)%scratch(i2, j2, k2, 1:3) + recvbuffer(jj + 1:jj + 3)
 
                jj = jj + 3
            end do
        end do completerecvs
 
        ! Complete the nonblocking sends.
 
        size = commpattern%nProcSend
        do i = 1, commpattern%nProcSend
            call mpi_waitany(size, sendrequests, index, mpistatus, ierr)
        end do
 
        ! Now we have accumulated back as far as xSeedd (stored in
        ! scratch). We can now do the whalo1to1_b
 
        ! Exchange the xSeeds in reverse.
        do nn = 1, ndom
            flowdoms(nn, level, sps)%realCommVars(1)%var => flowdoms(nn, level, sps)%scratch(:, :, :, 1)
            flowdoms(nn, level, sps)%realCommVars(2)%var => flowdoms(nn, level, sps)%scratch(:, :, :, 2)
            flowdoms(nn, level, sps)%realCommVars(3)%var => flowdoms(nn, level, sps)%scratch(:, :, :, 3)
        end do
 
        call whalo1to1realgeneric_b(3, level, sps, commpatterncell_2nd, internalcell_2nd)
 
        ! Finaly we can push back to the local x
        do nn = 1, ndom
            call setpointers_b(nn, level, sps)
 
            do k = 2, kl
                do j = 2, jl
                    do i = 2, il
                        ! Add is accumulate seed for xSeed (stored in scratch )
                        add = eighth * scratch(i, j, k, 1:3)
                        do kk = k - 1, k
                            do jj = j - 1, j
                                do ii = i - 1, i
                                    xd(ii, jj, kk, :) = xd(ii, jj, kk, :) + add
                                end do
                            end do
                        end do
                    end do
                end do
            end do
        end do
 
    end subroutine updateoversetconnectivity_b
#endif
end module oversetcommutilities