userSurfaceIntegrations.F90

Go to the documentation of this file.

module usersurfaceintegrations
 
    use constants
    use communication, only: commtype, internalcommtype
    use usersurfaceintegrationdata
contains
 
    subroutine integrateusersurfaces(localValues, famList, sps)
 
        use constants
        use block, onlY: flowdoms, ndom
        use flowvarrefstate, only: pref, rhoref, pref, timeref, lref, tref
        use communication, only: myid, adflow_comm_world
        use utils, only: echk, mynorm2
        use flowutils, only: computeptot, computettot
        use sorting, only: faminlist
        use zipperintegrations, only: flowintegrationzipper
        use utils, only: terminate
        implicit none
 
        ! Input Parameters
        real(kind=realtype), dimension(nLocalValues), intent(inout) :: localvalues
        integer(kind=intType), dimension(:), intent(in) :: famList
        integer(kind=intType), intent(in) :: sps
 
        ! Working parameters
        integer(kind=intType) :: iSurf, i, j, k, jj, ierr, nn, iDim, nPts
        real(kind=realtype), dimension(:), allocatable :: recvbuffer1, recvbuffer2
        real(kind=realtype), dimension(:, :), allocatable :: vars
        integer(kind=intType), dimension(:), allocatable :: fams
        logical, dimension(:), allocatable :: ptValid
        type(userintsurf), pointer :: surf
 
        if (nuserintsurfs == 0) then
            return ! Nothing to do
        end if
 
        ! Set the pointers for the required communication variables
        domainloop: do nn = 1, ndom
            if (flowdoms(nn, 1, sps)%addGridVelocities) then
                call terminate("userSurfaceIntegrations", "Cannot use user-supplied surface integrations"&
                     &"on with moving grids")
            end if
 
            flowdoms(nn, 1, sps)%realCommVars(irho)%var => flowdoms(nn, 1, sps)%w(:, :, :, irho)
            flowdoms(nn, 1, sps)%realCommVars(ivx)%var => flowdoms(nn, 1, sps)%w(:, :, :, ivx)
            flowdoms(nn, 1, sps)%realCommVars(ivy)%var => flowdoms(nn, 1, sps)%w(:, :, :, ivy)
            flowdoms(nn, 1, sps)%realCommVars(ivz)%var => flowdoms(nn, 1, sps)%w(:, :, :, ivz)
            flowdoms(nn, 1, sps)%realCommVars(izippflowp)%var => flowdoms(nn, 1, sps)%P(:, :, :)
            flowdoms(nn, 1, sps)%realCommVars(izippflowgamma)%var => flowdoms(nn, 1, sps)%gamma(:, :, :)
            ! flowDoms(nn, 1, sps)%realCommVars(iZippFlowSface)%var => Not Implemented
 
            flowdoms(nn, 1, sps)%realCommVars(izippflowx)%var => flowdoms(nn, 1, sps)%x(:, :, :, 1)
            flowdoms(nn, 1, sps)%realCommVars(izippflowy)%var => flowdoms(nn, 1, sps)%x(:, :, :, 2)
            flowdoms(nn, 1, sps)%realCommVars(izippflowz)%var => flowdoms(nn, 1, sps)%x(:, :, :, 3)
        end do domainloop
 
        masterloop: do isurf = 1, nuserintsurfs
 
            ! Pointer for easier reading
            surf => userintsurfs(isurf)
 
            ! We will make a short-cut here: By definition user supplied
            ! surfaces have a fixed family, we won't do anything if we
            ! are not told to deal with this surface.
 
            faminclude: if (faminlist(surf%famID, famlist)) then
                npts = size(surf%pts, 2)
 
                ! Communicate the face values and the nodal values
                if (myid == 0) then
                    allocate (recvbuffer1(6 * npts), recvbuffer2(3 * npts))
                else
                    allocate (recvbuffer1(0), recvbuffer2(0))
                end if
 
                call commuserintegrationsurfacevars(recvbuffer1, irho, izippflowgamma, surf%flowComm)
                call commuserintegrationsurfacevars(recvbuffer2, izippflowx, izippflowz, surf%nodeComm)
 
                ! *Finally* we can do the actual integrations
                if (myid == 0) then
 
                    ! Allocate some temporary data needed to supply to the
                    ! zipper integration routine.
                    allocate (ptvalid(npts), vars(npts, nzippflowcomm), fams(size(surf%conn, 2)))
 
                    ! Initialize ptValid to True. If we find that it isn't,
                    ! we'll permenantly set that point to false. This could
                    ! come from either the node or the flow comms.
 
                    ptvalid = .true.
                    ! Prepare for the "zipper" integration call. We have to
                    ! re-order the data according to the "inv" array in each
                    ! of the two comms.
                    do i = 1, npts
 
                        ! Flow Variables
                        j = surf%flowComm%inv(i)
                        vars(j, irho:izippflowgamma) = recvbuffer1(6 * (i - 1) + irho:6 * (i - 1) + izippflowgamma)
 
                        if (.not. surf%flowComm%valid(i)) then
                            ptvalid(j) = .false.
                        end if
 
                        ! Sface is not implemented. To correctly do this,
                        ! interpolate the three components of 's', do the dot
                        ! product with the local normal to get the sFace value.
                        vars(j, izippflowsface) = zero
 
                        ! Node Comm Values
                        j = surf%nodeComm%inv(i)
                        vars(j, izippflowx:izippflowz) = recvbuffer2(3 * i - 2:3 * i)
 
                        ! The additional pt-valid array
                        if (.not. surf%nodeComm%valid(i)) then
                            ptvalid(j) = .false.
                        end if
                    end do
 
                    ! The family array is all the same value:
                    fams = surf%famID
                    ! Perform the actual integration
                    call flowintegrationzipper(surf%isInflow, surf%conn, fams, vars, localvalues, famlist, sps, ptvalid)
                    deallocate (ptvalid, vars, fams)
                end if
                deallocate (recvbuffer1, recvbuffer2)
            end if faminclude
        end do masterloop
    end subroutine integrateusersurfaces
#ifndef USE_COMPLEX
    subroutine integrateusersurfaces_d(localValues, localValuesd, famList, sps)
 
        use constants
        use block, onlY: flowdoms, flowdomsd, ndom
        use flowvarrefstate, only: pref, rhoref, pref, timeref, lref, tref
        use communication, only: myid, adflow_comm_world
        use utils, only: echk, mynorm2
        use flowutils, only: computeptot, computettot
        use sorting, only: faminlist
        use zipperintegrations_d, only: flowintegrationzipper_d
        use utils, only: terminate
        implicit none
 
        ! Input Parameters
        real(kind=realtype), dimension(nLocalValues), intent(inout) :: localvalues, localvaluesd
        integer(kind=intType), dimension(:), intent(in) :: famList
        integer(kind=intType), intent(in) :: sps
 
        ! Working parameters
        integer(kind=intType) :: iSurf, i, j, k, jj, ierr, nn, iDim, nPts
        real(kind=realtype), dimension(:), allocatable :: recvbuffer1, recvbuffer2
        real(kind=realtype), dimension(:), allocatable :: recvbuffer1d, recvbuffer2d
        real(kind=realtype), dimension(:, :), allocatable :: vars, varsd
        integer(kind=intType), dimension(:), allocatable :: fams
        logical, dimension(:), allocatable :: ptValid
        type(userintsurf), pointer :: surf
 
        if (nuserintsurfs == 0) then
            return ! Nothing to do
        end if
 
        ! Set the pointers for the required communication variables
        domainloop: do nn = 1, ndom
            if (flowdoms(nn, 1, sps)%addGridVelocities) then
                call terminate("userSurfaceIntegrations", "Cannot use user-supplied surface integrations"&
                     &"on with moving grids")
            end if
 
            flowdoms(nn, 1, sps)%realCommVars(irho)%var => flowdoms(nn, 1, sps)%w(:, :, :, irho)
            flowdoms(nn, 1, sps)%realCommVars(ivx)%var => flowdoms(nn, 1, sps)%w(:, :, :, ivx)
            flowdoms(nn, 1, sps)%realCommVars(ivy)%var => flowdoms(nn, 1, sps)%w(:, :, :, ivy)
            flowdoms(nn, 1, sps)%realCommVars(ivz)%var => flowdoms(nn, 1, sps)%w(:, :, :, ivz)
            flowdoms(nn, 1, sps)%realCommVars(izippflowp)%var => flowdoms(nn, 1, sps)%P(:, :, :)
            flowdoms(nn, 1, sps)%realCommVars(izippflowgamma)%var => flowdoms(nn, 1, sps)%gamma(:, :, :)
            ! flowDoms(nn, 1, sps)%realCommVars(iZippFlowSface)%var => Not Implemented
 
            flowdoms(nn, 1, sps)%realCommVars(izippflowx)%var => flowdoms(nn, 1, sps)%x(:, :, :, 1)
            flowdoms(nn, 1, sps)%realCommVars(izippflowy)%var => flowdoms(nn, 1, sps)%x(:, :, :, 2)
            flowdoms(nn, 1, sps)%realCommVars(izippflowz)%var => flowdoms(nn, 1, sps)%x(:, :, :, 3)
 
            flowdoms(nn, 1, sps)%realCommVars(irho + nzippflowcomm)%var => flowdomsd(nn, 1, sps)%w(:, :, :, irho)
            flowdoms(nn, 1, sps)%realCommVars(ivx + nzippflowcomm)%var => flowdomsd(nn, 1, sps)%w(:, :, :, ivx)
            flowdoms(nn, 1, sps)%realCommVars(ivy + nzippflowcomm)%var => flowdomsd(nn, 1, sps)%w(:, :, :, ivy)
            flowdoms(nn, 1, sps)%realCommVars(ivz + nzippflowcomm)%var => flowdomsd(nn, 1, sps)%w(:, :, :, ivz)
            flowdoms(nn, 1, sps)%realCommVars(izippflowp + nzippflowcomm)%var => flowdomsd(nn, 1, sps)%P(:, :, :)
            flowdoms(nn, 1, sps)%realCommVars(izippflowgamma + nzippflowcomm)%var => &
                flowdomsd(nn, 1, sps)%gamma(:, :, :)
            ! flowDoms(nn, 1, sps)%realCommVars(iZippFlowSface+nZippFlowComm)%var => Not Implemented
 
            flowdoms(nn, 1, sps)%realCommVars(izippflowx + nzippflowcomm)%var => flowdomsd(nn, 1, sps)%x(:, :, :, 1)
            flowdoms(nn, 1, sps)%realCommVars(izippflowy + nzippflowcomm)%var => flowdomsd(nn, 1, sps)%x(:, :, :, 2)
            flowdoms(nn, 1, sps)%realCommVars(izippflowz + nzippflowcomm)%var => flowdomsd(nn, 1, sps)%x(:, :, :, 3)
 
        end do domainloop
 
        masterloop: do isurf = 1, nuserintsurfs
 
            ! Pointer for easier reading
            surf => userintsurfs(isurf)
 
            ! We will make a short-cut here: By definition user supplied
            ! surfaces have a fixed family, we won't do anything if we
            ! are not told to deal with this surface.
 
            faminclude: if (faminlist(surf%famID, famlist)) then
                npts = size(surf%pts, 2)
 
                ! Communicate the face values and the nodal values
                if (myid == 0) then
                    allocate (recvbuffer1(6 * npts), recvbuffer2(3 * npts))
                    allocate (recvbuffer1d(6 * npts), recvbuffer2d(3 * npts))
                else
                    allocate (recvbuffer1(0), recvbuffer2(0))
                    allocate (recvbuffer1d(0), recvbuffer2d(0))
                end if
 
                call commuserintegrationsurfacevars_d(recvbuffer1, recvbuffer1d, irho, izippflowgamma, surf%flowComm)
                call commuserintegrationsurfacevars_d(recvbuffer2, recvbuffer2d, izippflowx, izippflowz, surf%nodeComm)
 
                ! *Finally* we can do the actual integrations
                if (myid == 0) then
 
                    ! Allocate some temporary data needed to supply to the
                    ! zipper integration routine.
                    allocate (ptvalid(npts), vars(npts, nzippflowcomm), &
                              varsd(npts, nzippflowcomm), fams(size(surf%conn, 2)))
 
                    ! Initialize ptValid to True. If we find that it isn't,
                    ! we'll permenantly set that point to false. This could
                    ! come from either the node or the flow comms.
 
                    ptvalid = .true.
 
                    ! Prepare for the "zipper" integration call. We have to
                    ! re-order the data according to the "inv" array in each
                    ! of the two comms.
                    do i = 1, npts
 
                        ! Flow Variables
                        j = surf%flowComm%inv(i)
 
                        vars(j, irho:izippflowgamma) = recvbuffer1(6 * (i - 1) + irho:6 * (i - 1) + izippflowgamma)
                        varsd(j, irho:izippflowgamma) = recvbuffer1d(6 * (i - 1) + irho:6 * (i - 1) + izippflowgamma)
 
                        if (.not. surf%flowComm%valid(i)) then
                            ptvalid(j) = .false.
                        end if
 
                        ! Sface is not implemented. To correctly do this,
                        ! interpolate the three components of 's', do the dot
                        ! product with the local normal to get the sFace value.
                        vars(j, izippflowsface) = zero
                        varsd(j, izippflowsface) = zero
 
                        ! Node Comm Values
                        j = surf%nodeComm%inv(i)
                        vars(j, izippflowx:izippflowz) = recvbuffer2(3 * i - 2:3 * i)
                        varsd(j, izippflowx:izippflowz) = recvbuffer2d(3 * i - 2:3 * i)
 
                        ! The additional pt-valid array
                        if (.not. surf%nodeComm%valid(i)) then
                            ptvalid(j) = .false.
                        end if
                    end do
 
                    ! The family array is all the same value:
                    fams = surf%famID
 
                    ! Perform the actual integration
                    call flowintegrationzipper_d(surf%isInflow, surf%conn, fams, vars, &
                                                 varsd, localvalues, localvaluesd, &
                                                 famlist, sps, ptvalid)
                    deallocate (ptvalid, vars, varsd, fams)
                end if
                deallocate (recvbuffer1, recvbuffer2, recvbuffer1d, recvbuffer2d)
            end if faminclude
        end do masterloop
    end subroutine integrateusersurfaces_d
 
    subroutine integrateusersurfaces_b(localValues, localValuesd, famList, sps)
 
        use constants
        use block, onlY: flowdoms, flowdomsd, ndom
        use flowvarrefstate, only: pref, rhoref, pref, timeref, lref, tref
        use communication, only: myid, adflow_comm_world
        use utils, only: echk, mynorm2
        use flowutils, only: computeptot, computettot
        use sorting, only: faminlist
        use zipperintegrations_b, only: flowintegrationzipper_b
        use utils, only: terminate
        implicit none
 
        ! Input Parameters
        real(kind=realtype), dimension(nLocalValues), intent(inout) :: localvalues, localvaluesd
        integer(kind=intType), dimension(:), intent(in) :: famList
        integer(kind=intType), intent(in) :: sps
 
        ! Working parameters
        integer(kind=intType) :: iSurf, i, j, k, jj, ierr, nn, iDim, nPts
        real(kind=realtype), dimension(:), allocatable :: recvbuffer1, recvbuffer2
        real(kind=realtype), dimension(:), allocatable :: recvbuffer1d, recvbuffer2d
        real(kind=realtype), dimension(:, :), allocatable :: vars, varsd
        integer(kind=intType), dimension(:), allocatable :: fams
        logical, dimension(:), allocatable :: ptValid
        type(userintsurf), pointer :: surf
 
        if (nuserintsurfs == 0) then
            return ! Nothing to do
        end if
 
        ! Run the foward mode code pass:
        call integrateusersurfaces(localvalues, famlist, sps)
 
        ! Set the pointers for the required communication variables
        domainloop: do nn = 1, ndom
            if (flowdoms(nn, 1, sps)%addGridVelocities) then
                call terminate("userSurfaceIntegrations", "Cannot use user-supplied surface integrations"&
                     &"on with moving grids")
            end if
 
            flowdoms(nn, 1, sps)%realCommVars(irho)%var => flowdoms(nn, 1, sps)%w(:, :, :, irho)
            flowdoms(nn, 1, sps)%realCommVars(ivx)%var => flowdoms(nn, 1, sps)%w(:, :, :, ivx)
            flowdoms(nn, 1, sps)%realCommVars(ivy)%var => flowdoms(nn, 1, sps)%w(:, :, :, ivy)
            flowdoms(nn, 1, sps)%realCommVars(ivz)%var => flowdoms(nn, 1, sps)%w(:, :, :, ivz)
            flowdoms(nn, 1, sps)%realCommVars(izippflowp)%var => flowdoms(nn, 1, sps)%P(:, :, :)
            flowdoms(nn, 1, sps)%realCommVars(izippflowgamma)%var => flowdoms(nn, 1, sps)%gamma(:, :, :)
            ! flowDoms(nn, 1, sps)%realCommVars(iZippFlowSface)%var => Not Implemented
 
            flowdoms(nn, 1, sps)%realCommVars(izippflowx)%var => flowdoms(nn, 1, sps)%x(:, :, :, 1)
            flowdoms(nn, 1, sps)%realCommVars(izippflowy)%var => flowdoms(nn, 1, sps)%x(:, :, :, 2)
            flowdoms(nn, 1, sps)%realCommVars(izippflowz)%var => flowdoms(nn, 1, sps)%x(:, :, :, 3)
 
            flowdoms(nn, 1, sps)%realCommVars(irho + nzippflowcomm)%var => flowdomsd(nn, 1, sps)%w(:, :, :, irho)
            flowdoms(nn, 1, sps)%realCommVars(ivx + nzippflowcomm)%var => flowdomsd(nn, 1, sps)%w(:, :, :, ivx)
            flowdoms(nn, 1, sps)%realCommVars(ivy + nzippflowcomm)%var => flowdomsd(nn, 1, sps)%w(:, :, :, ivy)
            flowdoms(nn, 1, sps)%realCommVars(ivz + nzippflowcomm)%var => flowdomsd(nn, 1, sps)%w(:, :, :, ivz)
            flowdoms(nn, 1, sps)%realCommVars(izippflowp + nzippflowcomm)%var => flowdomsd(nn, 1, sps)%P(:, :, :)
            flowdoms(nn, 1, sps)%realCommVars(izippflowgamma + nzippflowcomm)%var => &
                flowdomsd(nn, 1, sps)%gamma(:, :, :)
            ! flowDoms(nn, 1, sps)%realCommVars(iZippFlowSface+nZippFlowComm)%var => Not Implemented
 
            flowdoms(nn, 1, sps)%realCommVars(izippflowx + nzippflowcomm)%var => flowdomsd(nn, 1, sps)%x(:, :, :, 1)
            flowdoms(nn, 1, sps)%realCommVars(izippflowy + nzippflowcomm)%var => flowdomsd(nn, 1, sps)%x(:, :, :, 2)
            flowdoms(nn, 1, sps)%realCommVars(izippflowz + nzippflowcomm)%var => flowdomsd(nn, 1, sps)%x(:, :, :, 3)
 
        end do domainloop
 
        masterloop: do isurf = 1, nuserintsurfs
 
            ! Pointer for easier reading
            surf => userintsurfs(isurf)
 
            ! We will make a short-cut here: By definition user supplied
            ! surfaces have a fixed family, we won't do anything if we
            ! are not told to deal with this surface.
 
            faminclude: if (faminlist(surf%famID, famlist)) then
                npts = size(surf%pts, 2)
 
                ! Communicate the face values and the nodal values
                if (myid == 0) then
                    allocate (recvbuffer1(6 * npts), recvbuffer2(3 * npts))
                    allocate (recvbuffer1d(6 * npts), recvbuffer2d(3 * npts))
                else
                    allocate (recvbuffer1(0), recvbuffer2(0))
                    allocate (recvbuffer1d(0), recvbuffer2d(0))
                end if
 
                call commuserintegrationsurfacevars(recvbuffer1, irho, izippflowgamma, surf%flowComm)
                call commuserintegrationsurfacevars(recvbuffer2, izippflowx, izippflowz, surf%nodeComm)
 
                ! *Finally* we can do the actual integrations
                if (myid == 0) then
 
                    ! Allocate some temporary data needed to supply to the
                    ! zipper integration routine.
                    allocate (ptvalid(npts), vars(npts, nzippflowcomm), &
                              varsd(npts, nzippflowcomm), fams(size(surf%conn, 2)))
 
                    ! Initialize ptValid to True. If we find that it isn't,
                    ! we'll permenantly set that point to false. This could
                    ! come from either the node or the flow comms.
 
                    ptvalid = .true.
 
                    varsd = zero
                    ! Prepare for the "zipper" integration call. We have to
                    ! re-order the data according to the "inv" array in each
                    ! of the two comms.
                    do i = 1, npts
 
                        ! Flow Variables
                        j = surf%flowComm%inv(i)
                        vars(j, irho:izippflowgamma) = recvbuffer1(6 * (i - 1) + irho:6 * (i - 1) + izippflowgamma)
 
                        if (.not. surf%flowComm%valid(i)) then
                            ptvalid(j) = .false.
                        end if
 
                        ! Sface is not implemented. To correctly do this,
                        ! interpolate the three components of 's', do the dot
                        ! product with the local normal to get the sFace value.
                        vars(j, izippflowsface) = zero
 
                        ! Node Comm Values
                        j = surf%nodeComm%inv(i)
                        vars(j, izippflowx:izippflowz) = recvbuffer2(3 * i - 2:3 * i)
 
                        ! The additional pt-valid array
                        if (.not. surf%nodeComm%valid(i)) then
                            ptvalid(j) = .false.
                        end if
                    end do
 
                    ! The family array is all the same value:
                    fams = surf%famID
 
                    ! Perform the actual (reverse) integration
                    call flowintegrationzipper_b(surf%isInflow, surf%conn, fams, vars, &
                                                 varsd, localvalues, localvaluesd, &
                                                 famlist, sps, ptvalid)
 
                    ! Accumulate into the receive buffers
                    recvbuffer1d = zero
                    recvbuffer2d = zero
                    do i = 1, npts
 
                        ! Accumulte back to the buffer
                        j = surf%flowComm%inv(i)
 
                        recvbuffer1d(6 * (i - 1) + irho:6 * (i - 1) + izippflowgamma) = &
                            recvbuffer1d(6 * (i - 1) + irho:6 * (i - 1) + izippflowgamma) + &
                            varsd(j, irho:izippflowgamma)
 
                        ! Sface is not implemented. No reverse seed. Just zero
 
                        ! Node Comm Values
                        j = surf%nodeComm%inv(i)
                        recvbuffer2d(3 * i - 2:3 * i) = recvbuffer2d(3 * i - 2:3 * i) + varsd(j, izippflowx:izippflowz)
                    end do
 
                    deallocate (ptvalid, vars, varsd, fams)
                end if
 
                ! Finish the reverse scatter
                call commuserintegrationsurfacevars_b(recvbuffer1, recvbuffer1d, irho, izippflowgamma, surf%flowComm)
                call commuserintegrationsurfacevars_b(recvbuffer2, recvbuffer2d, izippflowx, izippflowz, surf%nodeComm)
 
                deallocate (recvbuffer1, recvbuffer2, recvbuffer1d, recvbuffer2d)
            end if faminclude
        end do masterloop
    end subroutine integrateusersurfaces_b
 
#endif
 
    subroutine addintegrationsurface(pts, conn, famName, famID, isInflow, nPts, nConn)
        ! Add a user-supplied integration surface.
 
        use communication, only: myid
        use constants
 
        implicit none
 
        ! Input variables
        integer(kind=intType), intent(in) :: nPts, nConn, famID
        real(kind=realtype), dimension(3, nPts), intent(in) :: pts
        integer(kind=intType), dimension(3, nConn), intent(in) :: conn
        logical, intent(in) :: isInflow
        character(len=*) :: famName
        type(userintsurf), pointer :: surf
 
        ! Not really much to do here...we just have to save the data
        ! into the data structure untilly we actual have to do the
        ! search.
        nuserintsurfs = nuserintsurfs + 1
        if (nuserintsurfs > nuserintsurfsmax) then
            print *, "Error: Exceeded the maximum number of user-supplied "&
                 &"integration slices. Increase nUserIntSurfsMax"
            stop
        end if
 
        surf => userintsurfs(nuserintsurfs)
        if (myid == 0) then
            allocate (surf%pts(3, npts), surf%conn(3, nconn))
            surf%pts = pts
            surf%conn = conn
        end if
        surf%famName = famname
        surf%famID = famid
        surf%isInflow = isinflow
    end subroutine addintegrationsurface
 
    subroutine buildvolumeadts(oBlocks, useDual)
 
        ! This builds volume ADTs for the the owned blocks. It will build
        ! either the dual mesh or the primal mesh depening on the flag
        ! useDual.
 
        use constants
        use oversetdata, only: oversetblock
        use blockpointers, only: ndom, x, ie, je, ke, il, jl, kl, vol, ib, jb, kb, &
                                 iblank, bcdata, nbocos, bcfaceid, bctype
        use adtbuild, only: buildserialhex
        use utils, only: setpointers, echk
        implicit none
 
        ! Input/Output Parameters
        type(oversetblock), dimension(:), target :: oBlocks
        logical :: useDual
 
        ! Working Parameters
        integer(kind=intType) :: nInterpol, nn, i, j, k, iii, jjj, kkk
        integer(kind=intType) :: iStart, jStart, kStart, iEnd, jEnd, kEnd
        integer(kind=intType) :: ii, jj, kk, mm, nADT, nHexa, planeOffset
        type(oversetblock), pointer :: oBlock
 
        ninterpol = 1 ! we get the ADT to compute the interpolated volume for us.
 
        domainloop: do nn = 1, ndom
 
            call setpointers(nn, 1, 1)
            oblock => oblocks(nn)
 
            primalordual: if (usedual) then
 
                ! Now setup the data for the ADT
                nhexa = il * jl * kl
                nadt = ie * je * ke
                oblock%il = il
                oblock%jl = jl
                oblock%kl = kl
 
                allocate (oblock%xADT(3, nadt), oblock%hexaConn(8, nhexa), &
                          oblock%qualDonor(1, nadt))
                ! Fill up the xADT using cell centers (dual mesh)
                mm = 0
                do k = 1, ke
                    do j = 1, je
                        do i = 1, ie
                            mm = mm + 1
                            oblock%xADT(:, mm) = 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, :))
                            oblock%qualDonor(1, mm) = vol(i, j, k)
                        end do
                    end do
                end do
 
                do mm = 1, nbocos
                    ! We need to make sure that the interpolation does not
                    ! use a halo behind an overset outer bound. This could
                    ! happen because the last cell is interpolated (-1) and
                    ! the iblank on the halos are still 1. Just set the
                    ! quality very high so it is not accepted.
 
                    select case (bcfaceid(mm))
                    case (imin)
                        istart = 1; iend = 1;
                        jstart = bcdata(mm)%icBeg; jend = bcdata(mm)%icEnd
                        kstart = bcdata(mm)%jcBeg; kend = bcdata(mm)%jcEnd
                    case (imax)
                        istart = ie; iend = ie;
                        jstart = bcdata(mm)%icBeg; jend = bcdata(mm)%icEnd
                        kstart = bcdata(mm)%jcBeg; kend = bcdata(mm)%jcEnd
                    case (jmin)
                        istart = bcdata(mm)%icBeg; iend = bcdata(mm)%icEnd
                        jstart = 1; jend = 1
                        kstart = bcdata(mm)%jcBeg; kend = bcdata(mm)%jcEnd
                    case (jmax)
                        istart = bcdata(mm)%icBeg; iend = bcdata(mm)%icEnd
                        jstart = je; jend = je;
                        kstart = bcdata(mm)%jcBeg; kend = bcdata(mm)%jcEnd
                    case (kmin)
                        istart = bcdata(mm)%icBeg; iend = bcdata(mm)%icEnd
                        jstart = bcdata(mm)%jcBeg; jend = bcdata(mm)%jcEnd
                        kstart = 1; kend = 1;
                    case (kmax)
                        istart = bcdata(mm)%icBeg; iend = bcdata(mm)%icEnd
                        jstart = bcdata(mm)%jcBeg; jend = bcdata(mm)%jcEnd
                        kstart = ke; kend = ke;
                    end select
 
                    if (bctype(mm) == oversetouterbound) then
                        do k = kstart, kend
                            do j = jstart, jend
                                do i = istart, iend
                                    ! recompute the index
                                    kk = (k - 1) * ie * je + (j - 1) * ie + i
                                    oblock%qualDonor(1, kk) = large
                                end do
                            end do
                        end do
                    end if
                end do
 
                mm = 0
                ! These are the 'elements' of the dual mesh.
                planeoffset = ie * je
                do k = 2, ke
                    do j = 2, je
                        do i = 2, ie
                            mm = mm + 1
                            oblock%hexaConn(1, mm) = (k - 2) * planeoffset + (j - 2) * ie + (i - 2) + 1
                            oblock%hexaConn(2, mm) = oblock%hexaConn(1, mm) + 1
                            oblock%hexaConn(3, mm) = oblock%hexaConn(2, mm) + ie
                            oblock%hexaConn(4, mm) = oblock%hexaConn(3, mm) - 1
 
                            oblock%hexaConn(5, mm) = oblock%hexaConn(1, mm) + planeoffset
                            oblock%hexaConn(6, mm) = oblock%hexaConn(2, mm) + planeoffset
                            oblock%hexaConn(7, mm) = oblock%hexaConn(3, mm) + planeoffset
                            oblock%hexaConn(8, mm) = oblock%hexaConn(4, mm) + planeoffset
                        end do
                    end do
                end do
            else
                ! Note that we will be including the halo primal cells. This
                ! should slightly increase robusness for viscous off-wall
                ! spacing. This means the primal mesh has 1 MORE node/cell
                ! in each direction.
 
                ! Now setup the data for the ADT
                nhexa = ie * je * ke
                nadt = ib * jb * kb
                oblock%il = ie
                oblock%jl = je
                oblock%kl = ke
 
                allocate (oblock%xADT(3, nadt), oblock%hexaConn(8, nhexa), &
                          oblock%qualDonor(1, nadt))
 
                oblock%qualDonor = zero
                ! Fill up the xADT using the primal nodes
                mm = 0
                do k = 0, ke
                    do j = 0, je
                        do i = 0, ie
                            mm = mm + 1
                            oblock%xADT(:, mm) = x(i, j, k, :)
 
                            ! Since we don't have all 8 volumes surrounding the
                            ! halo nodes, clip the volumes to be between 0 and ib etc.
                            do iii = 0, 1
                                do jjj = 0, 1
                                    do kkk = 0, 1
                                        ii = min(max(0, iii + i), ib)
                                        jj = min(max(0, jjj + j), jb)
                                        kk = min(max(0, kkk + k), kb)
 
                                        oblock%qualDonor(1, mm) = oblock%qualDonor(1, mm) + &
                                                                  vol(ii, jj, kk)
                                    end do
                                end do
                            end do
 
                            ! Dividing by 8 isn't strictly necessary but we'll
                            ! do it anyway.
                            oblock%qualDonor(1, mm) = oblock%qualDonor(1, mm) * eighth
                        end do
                    end do
                end do
 
                mm = 0
                ! These are the 'elements' of the dual mesh.
                planeoffset = ib * jb
                do k = 1, ke
                    do j = 1, je
                        do i = 1, ie
                            mm = mm + 1
                            oblock%hexaConn(1, mm) = (k - 1) * planeoffset + (j - 1) * ib + (i - 1) + 1
                            oblock%hexaConn(2, mm) = oblock%hexaConn(1, mm) + 1
                            oblock%hexaConn(3, mm) = oblock%hexaConn(2, mm) + ib
                            oblock%hexaConn(4, mm) = oblock%hexaConn(3, mm) - 1
 
                            oblock%hexaConn(5, mm) = oblock%hexaConn(1, mm) + planeoffset
                            oblock%hexaConn(6, mm) = oblock%hexaConn(2, mm) + planeoffset
                            oblock%hexaConn(7, mm) = oblock%hexaConn(3, mm) + planeoffset
                            oblock%hexaConn(8, mm) = oblock%hexaConn(4, mm) + planeoffset
                        end do
                    end do
                end do
            end if primalordual
 
            ! Call the custom build routine -- Serial only, only Hexa volumes,
            ! we supply our own ADT Type
 
            call buildserialhex(nhexa, nadt, oblock%xADT, oblock%hexaConn, oblock%ADT)
        end do domainloop
 
    end subroutine buildvolumeadts
 
    subroutine performinterpolation(pts, oBlocks, useDual, comm)
 
        ! This routine performs the actual searches for the slices. It is
        ! generic in the sense that it will search an arbtitrary set of
        ! points on either the primal or dual meshes. The final required
        ! communication data is then written into the supplied comm.
 
        use constants
        use block, only: interppttype
        use communication, only: adflow_comm_world, myid, nproc
        use oversetdata, only: oversetblock
        use blockpointers, only: ndom, x, ie, je, ke, il, jl, kl, x, iblank, vol
        use adtlocalsearch, only: mindistancetreesearchsinglepoint, &
                                  containmenttreesearchsinglepoint
        use adtutils, only: stack
        use adtdata, only: adtbboxtargettype
        use utils, only: setpointers, mynorm2, echk
        use inputoverset, only: oversetprojtol
        use oversetutilities, only: fractoweights2, getcumulativeform
 
        implicit none
 
        ! Input parameters:
        real(kind=realtype), dimension(:, :), intent(in) :: pts
        type(userintsurf) :: surf
        type(oversetblock), dimension(:), target, intent(in) :: oBlocks
        logical, intent(in) :: useDual
        !type(oversetSurf), dimension(:), intent(in) :: oSurfs
        type(usersurfcommtype) :: comm
 
        ! Working parameters
        type(oversetblock), pointer :: oBlock
        type(interppttype), dimension(:), allocatable :: surfFringes
 
        integer(Kind=intType) :: i, j, k, ii, jj, kk, iii, jjj, kkk, nn, mm
        integer(kind=intType) :: iSurf, ierr, nInterpol, iProc
        integer(kind=intType) :: nHexa, nAdt, planeOffset, elemID, nPts
        real(kind=realtype) :: xc(4), weight(8)
        integer mpiStatus(MPI_STATUS_SIZE)
 
        real(kind=realtype) :: uvw(5), uvw2(5), donorqual, xcheck(3)
        integer(kind=intType) :: intInfo(3), intInfo2(3)
        logical :: failed, invalid
 
        integer(kind=intType), dimension(:, :), allocatable :: donorInfo, intSend
        integer(kind=intType), dimension(:), allocatable :: procSizes
        real(kind=realtype), dimension(:, :), allocatable :: donorfrac, realsend
 
        ! Variables we have to pass the ADT search routine
        integer(kind=intType), dimension(:), pointer :: BB
        type(adtbboxtargettype), dimension(:), pointer :: BB2
        integer(kind=intType), dimension(:), pointer :: frontLeaves
        integer(kind=intType), dimension(:), pointer :: frontLeavesNew
 
        ! Data for the search
        allocate (bb(20), frontleaves(25), frontleavesnew(25), stack(100), bb2(20))
 
        npts = size(pts, 2)
 
        ! Allocate donor information arrays
        allocate (donorfrac(4, npts), donorinfo(5, npts))
        donorinfo = -1
        donorfrac(4, :) = large
        ninterpol = 1
        domainsearch: do nn = 1, ndom
            oblock => oblocks(nn)
            call setpointers(nn, 1, 1)
 
            ! Search the supplied pts one at a time
            elemloop: do i = 1, npts
 
                xc(1:3) = pts(:, i)
 
                ! Call the standard tree search
                call containmenttreesearchsinglepoint(oblock%ADT, xc, intinfo, uvw, &
                                                      oblock%qualDonor, ninterpol, bb, &
                                                      frontleaves, frontleavesnew, failed)
 
                ! Make sure this point is not garbage.
                if (intinfo(1) >= 0) then
                    call fractoweights2(uvw(1:3), weight)
                    xcheck = zero
                    do j = 1, 8
                        xcheck = xcheck + weight(j) * oblock%xADT(:, oblock%hexaConn(j, intinfo(3)))
                    end do
 
                    if (mynorm2(xcheck - xc(1:3)) > oversetprojtol) then
                        failed = .true.
                    end if
                end if
 
                if (intinfo(1) >= 0 .and. failed) then
                    ! we "found" a point but it is garbage. Do the failsafe search
                    xc(4) = large
                    call mindistancetreesearchsinglepoint(oblock%ADT, xc, intinfo, uvw, &
                                                          oblock%qualDonor, ninterpol, bb2, frontleaves, frontleavesnew)
 
                    ! Check this one:
                    call fractoweights2(uvw(1:3), weight)
                    xcheck = zero
                    do j = 1, 8
                        xcheck = xcheck + weight(j) * oblock%xADT(:, oblock%hexaConn(j, intinfo(3)))
                    end do
 
                    ! Since this is the last line of defence, relax the tolerance a bit
                    if (mynorm2(xcheck - xc(1:3)) > 100 * oversetprojtol) then
                        ! This fringe has not found a donor
                        intinfo(1) = -1
                    else
                        ! This one has now passed.
 
                        ! Important! uvw(4) is the distance squared for this search
                        ! not interpolated value
                        uvw(4) = uvw(5)
                    end if
                end if
 
                elemfound: if (intinfo(1) >= 0) then
 
                    ! Donor and block and index information for this donor.
                    donorqual = uvw(4)
                    elemid = intinfo(3) - 1 ! Make it zero based
 
                    ! The dual mesh needs an offset of 1 becuase it only used
                    ! 1:ie values. This is not necessary for the primal.
                    if (usedual) then
                        ii = mod(elemid, oblock%il) + 1
                        jj = mod(elemid / oblock%il, oblock%jl) + 1
                        kk = elemid / (oblock%il * oblock%jl) + 1
                    else
                        ii = mod(elemid, oblock%il)
                        jj = mod(elemid / oblock%il, oblock%jl)
                        kk = elemid / (oblock%il * oblock%jl)
                    end if
                    ! Rememebr donorFrac(4, i) is the current best quality
                    if (donorqual < donorfrac(4, i)) then
 
                        invalid = .false.
 
                        ! For the dual mesh search, we have to make sure the
                        ! potential donors are valid. Such a check is not
                        ! necessary for the primal search since all nodes are
                        ! considered valid.
 
                        if (usedual) then
                            ! Check if the point is invalid. We can do this
                            ! with i-blank array. We can only accept compute
                            ! cells (iblank=1) or interpolated
                            ! cells(iblank=-1).
                            do kkk = 0, 1
                                do jjj = 0, 1
                                    do iii = 0, 1
                                        if (.not. (iblank(ii + iii, jj + jjj, kk + kkk) == 1 .or. &
                                                   iblank(ii + iii, jj + jjj, kk + kkk) == -1)) then
                                            invalid = .true.
                                        end if
                                    end do
                                end do
                            end do
                        end if
 
                        if (.not. invalid) then
 
                            ! Set the quality of the donor to the one we
                            ! just found. Save the rest of the necessary
                            ! information.
                            donorinfo(1, i) = myid
                            donorinfo(2, i) = nn
                            donorinfo(3, i) = ii
                            donorinfo(4, i) = jj
                            donorinfo(5, i) = kk
                            donorfrac(1:3, i) = uvw(1:3)
                            donorfrac(4, i) = donorqual
                        end if
                    end if
                end if elemfound
            end do elemloop
        end do domainsearch
 
        ! Next count up the number of valid donors we've found and compact
        ! the info back to that length.
        if (myid /= 0) then
            j = 0
            do i = 1, npts
                if (donorinfo(1, i) /= -1) then
                    j = j + 1
                end if
            end do
            allocate (intsend(6, j), realsend(4, i))
            if (j > 0) then
                j = 0
                do i = 1, npts
                    if (donorinfo(1, i) /= -1) then
                        j = j + 1
                        intsend(1:5, j) = donorinfo(:, i)
                        intsend(6, j) = i
                        realsend(:, j) = donorfrac(:, i)
                    end if
                end do
            end if
        else
            ! On the root proc, use intSend and realSend as the receiver
            ! buffer. These can be at most nPts sized.
            allocate (intsend(6, npts), realsend(4, npts))
        end if
 
        ! Gather up the sizes (j) to the root processor so he know who to
        ! expect data from.
        allocate (procsizes(0:nproc - 1))
 
        call mpi_gather(j, 1, adflow_integer, procsizes, 1, &
                        adflow_integer, 0, adflow_comm_world, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Next all the procs need to send all the information back to the
        ! root processor where we will determine the proper donors for
        ! each of the cells
 
        ! All procs except root fire off their data.
        if (myid >= 1) then
            if (j > 0) then
                call mpi_send(intsend, j * 6, adflow_integer, 0, myid, &
                              adflow_comm_world, ierr)
                call echk(ierr, __file__, __line__)
 
                call mpi_send(realsend, j * 4, adflow_real, 0, myid, &
                              adflow_comm_world, ierr)
                call echk(ierr, __file__, __line__)
            end if
        end if
 
        ! And the root processor recieves it...
        if (myid == 0) then
            do iproc = 1, nproc - 1
                ! Determine if this proc has sent anything:
                if (procsizes(iproc) /= 0) then
 
                    call mpi_recv(intsend, 6 * npts, adflow_integer, iproc, iproc, &
                                  adflow_comm_world, mpistatus, ierr)
                    call echk(ierr, __file__, __line__)
 
                    call mpi_recv(realsend, 4 * npts, adflow_real, iproc, iproc, &
                                  adflow_comm_world, mpistatus, ierr)
                    call echk(ierr, __file__, __line__)
 
                    ! Now process the data (intSend and realSend) that we
                    ! just received. We don't need to check the status for
                    ! the sizes becuase we already know the sizes from the
                    ! initial gather we did.
 
                    do i = 1, procsizes(iproc)
                        ii = intsend(6, i)
 
                        if (realsend(4, i) < donorfrac(4, ii)) then
                            ! The incoming quality is better. Accept it.
                            donorinfo(1:5, ii) = intsend(1:5, i)
                            donorfrac(:, ii) = realsend(:, i)
                        end if
                    end do
                end if
            end do
 
            ! To make this easier, convert the information we have to a
            ! 'fringeType' array so we can use the pre-existing sorting
            ! routine.
            allocate (surffringes(npts))
            do i = 1, npts
                surffringes(i)%donorProc = donorinfo(1, i)
                surffringes(i)%donorBlock = donorinfo(2, i)
                surffringes(i)%dI = donorinfo(3, i)
                surffringes(i)%dJ = donorinfo(4, i)
                surffringes(i)%dK = donorinfo(5, i)
                surffringes(i)%donorFrac = donorfrac(1:3, i)
                ! Use the myBlock attribute to keep track of the original
                ! index. When we sort the fringes, they will no longer be
                ! in the same order
                surffringes(i)%myBlock = i
            end do
 
            ! Perform the actual sort.
            call qsortinterppttype(surffringes, npts)
 
            ! We will reuse-proc sizes to now mean the number of elements
            ! that the processor *actually* has to send. We will include
            ! the root proc itself in the calc becuase that will tell us
            ! the size of the internal comm structure.
 
            procsizes = 0
            allocate (comm%valid(npts))
            comm%valid = .true.
            do i = 1, npts
                if (surffringes(i)%donorProc < 0) then
                    ! We dont have a donor. Flag this point as invalid
                    comm%valid(i) = .false.
                end if
 
                ! Dump the points without donors on the root proc by making
                ! sure j is at least 0 for the root proc. These will just
                ! simply be ignored during the comm.
                j = max(surffringes(i)%donorProc, 0)
                procsizes(j) = procsizes(j) + 1
            end do
        end if
 
        ! Simply broadcast out the the proc sizes back to everyone so all
        ! processors know if they are to receive anything back.
        call mpi_bcast(procsizes, nproc, adflow_integer, 0, adflow_comm_world, ierr)
        call echk(ierr, __file__, __line__)
 
        ! We can now save some of the final comm data required on the
        ! root proc for this surf.
        allocate (comm%procSizes(0:nproc - 1), comm%procDisps(0:nproc))
 
        ! Copy over procSizes and generate the cumulative form of
        ! the size array, procDisps
        comm%procSizes = procsizes
        call getcumulativeform(comm%procSizes, nproc, comm%procDisps)
 
        ! Record the elemInverse which is necessary to index into
        ! the original conn array.
        if (myid == 0) then
            allocate (comm%Inv(npts))
            do i = 1, npts
                comm%Inv(i) = surffringes(i)%myBlock
            end do
        end if
 
        ! Now we can send out the final donor information to the
        ! processors that must supply it.
        comm%nDonor = procsizes(myid)
        allocate (comm%frac(3, comm%nDonor), comm%donorInfo(4, comm%nDonor))
 
        if (myid >= 1) then
            if (comm%nDonor > 0) then
                ! We are responible for at least 1 donor. We have to make
                ! use of the intSend and realSend buffers again (which
                ! are guaranteed to be big enough). The reason we can't
                ! dump the data in directlyis that intSend and realSend
                ! have a different leading index than we need on the
                ! final data structure.
 
                call mpi_recv(intsend, 6 * comm%nDonor, adflow_integer, 0, myid, &
                              adflow_comm_world, mpistatus, ierr)
                call echk(ierr, __file__, __line__)
 
                call mpi_recv(realsend, 4 * comm%nDonor, adflow_real, 0, myid, &
                              adflow_comm_world, mpistatus, ierr)
                call echk(ierr, __file__, __line__)
 
                ! Copy into final structure
                do i = 1, comm%nDonor
                    comm%donorInfo(:, i) = intsend(1:4, i)
                    comm%frac(:, i) = realsend(1:3, i)
                end do
            end if
        else
            ! We are the root processor.
            if (comm%nDonor > 0) then
                ! We need to copy out our donor info on the root proc if we have any
                do i = comm%procDisps(myid) + 1, comm%procDisps(myid + 1)
                    comm%donorInfo(1, i) = surffringes(i)%donorBlock
                    comm%donorInfo(2, i) = surffringes(i)%dI
                    comm%donorInfo(3, i) = surffringes(i)%dJ
                    comm%donorInfo(4, i) = surffringes(i)%dK
                    comm%frac(1:3, i) = surffringes(i)%donorFrac
                end do
            end if
 
            ! Now loop over the rest of the procs and send out the info we
            ! need. We have to temporarily copy the data back out of
            ! fringes to the intSend and realSend arrays
            do iproc = 1, nproc - 1
 
                if (comm%procSizes(iproc) > 0) then
                    ! Have something to send here:
                    j = 0
                    do i = comm%procDisps(iproc) + 1, comm%procDisps(iproc + 1)
                        j = j + 1
 
                        intsend(1, j) = surffringes(i)%donorBlock
                        intsend(2, j) = surffringes(i)%dI
                        intsend(3, j) = surffringes(i)%dJ
                        intsend(4, j) = surffringes(i)%dK
                        realsend(1:3, j) = surffringes(i)%donorFrac
                    end do
 
                    call mpi_send(intsend, j * 6, adflow_integer, iproc, iproc, &
                                  adflow_comm_world, ierr)
                    call echk(ierr, __file__, __line__)
 
                    call mpi_send(realsend, j * 4, adflow_real, iproc, iproc, &
                                  adflow_comm_world, ierr)
                    call echk(ierr, __file__, __line__)
                end if
            end do
 
            ! Deallocate data allocatd only on root proc
            deallocate (surffringes)
        end if
 
        ! Nuke rest of allocated on all procs
        deallocate (intsend, realsend, procsizes, donorinfo, donorfrac)
        deallocate (bb, frontleaves, frontleavesnew, stack, bb2)
    end subroutine performinterpolation
 
    subroutine interpolateintegrationsurfaces
 
        ! This routine performs the actual searches for the slices. We
        ! reuse much of the same machinery as is used in the overset code.
 
        use constants
        use communication, only: adflow_comm_world, myid
        use oversetdata, only: oversetblock
        use blockpointers, only: ndom, ie, je, ke, il, jl, kl
        use adtbuild, only: buildserialhex, destroyserialhex
        use utils, only: setpointers, echk
 
        implicit none
 
        ! Working parameters
        type(oversetblock), dimension(nDom), target :: oBlocks
        type(userintsurf), pointer :: surf
 
        integer(Kind=intType) :: iSurf, ii, i, nn, nPts, ierr
        real(kind=realtype), dimension(:, :), allocatable :: pts
        logical :: useDual
        if (nuserintsurfs == 0) then
            return
        end if
        primaldualloop: do ii = 1, 2
            if (ii == 1) then
                usedual = .true.
            else
                usedual = .false.
            end if
 
            call buildvolumeadts(oblocks, usedual)
 
            masterloop: do isurf = 1, nuserintsurfs
 
                surf => userintsurfs(isurf)
 
                if (myid == 0) then
 
                    ! We are interpolating the nodal values for both the
                    ! nodes and the solution variables.
                    npts = size(surf%pts, 2)
                    allocate (pts(3, npts))
                    nodeloop: do i = 1, npts
                        pts(:, i) = surf%pts(:, i)
                    end do nodeloop
                end if
 
                ! Send the number of points back to all procs:
                call mpi_bcast(npts, 1, adflow_integer, 0, adflow_comm_world, ierr)
                call echk(ierr, __file__, __line__)
 
                ! All other procs except the root allocate space and receive
                ! the pt array.
                if (myid /= 0) then
                    allocate (pts(3, npts))
                end if
 
                call mpi_bcast(pts, 3 * npts, adflow_real, 0, adflow_comm_world, ierr)
                call echk(ierr, __file__, __line__)
 
                ! Call the actual interpolation routine
                if (ii == 1) then
                    call performinterpolation(pts, oblocks, .true., surf%flowComm)
                else
                    call performinterpolation(pts, oblocks, .false., surf%nodeComm)
                end if
 
                deallocate (pts)
            end do masterloop
 
            ! Destroy the ADT Data and allocated values
            do nn = 1, ndom
                call destroyserialhex(oblocks(nn)%ADT)
                deallocate (oblocks(nn)%xADT, oblocks(nn)%hexaConn, oblocks(nn)%qualDonor)
            end do
        end do primaldualloop
    end subroutine interpolateintegrationsurfaces
 
    subroutine commuserintegrationsurfacevars(recvBuffer, varStart, varEnd, comm)
 
        use constants
        use block, onlY: flowdoms, ndom
        use communication, only: myid, adflow_comm_world
        use utils, only: echk
        use oversetutilities, only: fractoweights
 
        implicit none
 
        ! Input/Output
        real(kind=realtype), dimension(:) :: recvbuffer
        integer(kind=intType), intent(in) :: varStart, varEnd
        type(usersurfcommtype) :: comm
 
        ! Working
        real(kind=realtype), dimension(:), allocatable :: sendbuffer
        integer(Kind=intType) :: d1, i1, j1, k1, jj, k, nvar, i, ierr
        real(kind=realtype), dimension(8) :: weight
 
        ! The number of variables we are transferring:
        nvar = varend - varstart + 1
 
        ! We assume that the pointers to the realCommVars have already been set.
 
        allocate (sendbuffer(nvar * comm%nDonor))
        ! initialize the sendBuffer to an arbitrary value. If the
        ! integration tries to use this value, you know something is
        ! wrong.
        sendbuffer = -99999
        ! First generate the interpolated data necessary
        jj = 0
        donorloop: do i = 1, comm%nDonor
            ! Convert the frac to weights
            call fractoweights(comm%frac(:, i), weight)
 
            ! Block and indices for easier reading. The +1 is due to the
            ! pointer offset on realCommVars.
 
            d1 = comm%donorInfo(1, i) ! Block Index
            i1 = comm%donorInfo(2, i) + 1 ! donor I index
            j1 = comm%donorInfo(3, i) + 1 ! donor J index
            k1 = comm%donorInfo(4, i) + 1 ! donor K index
 
            ! We are interpolating nVar variables
            do k = varstart, varend
                jj = jj + 1
                if (d1 > 0) then ! Is this pt valid?
                    sendbuffer(jj) = &
                        weight(1) * flowdoms(d1, 1, 1)%realCommVars(k)%var(i1, j1, k1) + &
                        weight(2) * flowdoms(d1, 1, 1)%realCommVars(k)%var(i1 + 1, j1, k1) + &
                        weight(3) * flowdoms(d1, 1, 1)%realCommVars(k)%var(i1, j1 + 1, k1) + &
                        weight(4) * flowdoms(d1, 1, 1)%realCommVars(k)%var(i1 + 1, j1 + 1, k1) + &
                        weight(5) * flowdoms(d1, 1, 1)%realCommVars(k)%var(i1, j1, k1 + 1) + &
                        weight(6) * flowdoms(d1, 1, 1)%realCommVars(k)%var(i1 + 1, j1, k1 + 1) + &
                        weight(7) * flowdoms(d1, 1, 1)%realCommVars(k)%var(i1, j1 + 1, k1 + 1) + &
                        weight(8) * flowdoms(d1, 1, 1)%realCommVars(k)%var(i1 + 1, j1 + 1, k1 + 1)
                end if
            end do
        end do donorloop
 
        ! Now we can do an mpi_gatherv to the root proc:
        call mpi_gatherv(sendbuffer, nvar * comm%nDonor, adflow_real, recvbuffer, &
                         nvar * comm%procSizes, nvar * comm%procDisps, adflow_real, 0, adflow_comm_world, ierr)
        call echk(ierr, __file__, __line__)
        deallocate (sendbuffer)
 
    end subroutine commuserintegrationsurfacevars
 
    subroutine commuserintegrationsurfacevars_d(recvBuffer, recvBufferd, varStart, varEnd, comm)
 
        use constants
        use block, onlY: flowdoms, ndom
        use communication, only: myid, adflow_comm_world
        use utils, only: echk
        use oversetutilities, only: fractoweights
 
        implicit none
 
        ! Input/Output
        real(kind=realtype), dimension(:) :: recvbuffer, recvbufferd
        integer(kind=intType), intent(in) :: varStart, varEnd
        type(usersurfcommtype) :: comm
 
        ! Working
        real(kind=realtype), dimension(:), allocatable :: sendbuffer, sendbufferd
        integer(Kind=intType) :: d1, i1, j1, k1, jj, k, nvar, i, ierr
        real(kind=realtype), dimension(8) :: weight
 
        ! The number of variables we are transferring:
        nvar = varend - varstart + 1
 
        ! We assume that the pointers to the realCommVars have already been set.
 
        allocate (sendbuffer(nvar * comm%nDonor), &
                  sendbufferd(nvar * comm%nDonor))
 
        ! First generate the interpolated data necessary
        jj = 0
        donorloop: do i = 1, comm%nDonor
            ! Convert the frac to weights
            call fractoweights(comm%frac(:, i), weight)
 
            ! Block and indices for easier reading. The +1 is due to the
            ! pointer offset on realCommVars.
 
            d1 = comm%donorInfo(1, i) ! Block Index
            i1 = comm%donorInfo(2, i) + 1 ! donor I index
            j1 = comm%donorInfo(3, i) + 1 ! donor J index
            k1 = comm%donorInfo(4, i) + 1 ! donor K index
 
            ! We are interpolating nVar variables
            do k = varstart, varend
                jj = jj + 1
                if (d1 > 0) then ! Is this pt valid?
                    sendbuffer(jj) = &
                        weight(1) * flowdoms(d1, 1, 1)%realCommVars(k)%var(i1, j1, k1) + &
                        weight(2) * flowdoms(d1, 1, 1)%realCommVars(k)%var(i1 + 1, j1, k1) + &
                        weight(3) * flowdoms(d1, 1, 1)%realCommVars(k)%var(i1, j1 + 1, k1) + &
                        weight(4) * flowdoms(d1, 1, 1)%realCommVars(k)%var(i1 + 1, j1 + 1, k1) + &
                        weight(5) * flowdoms(d1, 1, 1)%realCommVars(k)%var(i1, j1, k1 + 1) + &
                        weight(6) * flowdoms(d1, 1, 1)%realCommVars(k)%var(i1 + 1, j1, k1 + 1) + &
                        weight(7) * flowdoms(d1, 1, 1)%realCommVars(k)%var(i1, j1 + 1, k1 + 1) + &
                        weight(8) * flowdoms(d1, 1, 1)%realCommVars(k)%var(i1 + 1, j1 + 1, k1 + 1)
                    sendbufferd(jj) = &
                        weight(1) * flowdoms(d1, 1, 1)%realCommVars(k + nzippflowcomm)%var(i1, j1, k1) + &
                        weight(2) * flowdoms(d1, 1, 1)%realCommVars(k + nzippflowcomm)%var(i1 + 1, j1, k1) + &
                        weight(3) * flowdoms(d1, 1, 1)%realCommVars(k + nzippflowcomm)%var(i1, j1 + 1, k1) + &
                        weight(4) * flowdoms(d1, 1, 1)%realCommVars(k + nzippflowcomm)%var(i1 + 1, j1 + 1, k1) + &
                        weight(5) * flowdoms(d1, 1, 1)%realCommVars(k + nzippflowcomm)%var(i1, j1, k1 + 1) + &
                        weight(6) * flowdoms(d1, 1, 1)%realCommVars(k + nzippflowcomm)%var(i1 + 1, j1, k1 + 1) + &
                        weight(7) * flowdoms(d1, 1, 1)%realCommVars(k + nzippflowcomm)%var(i1, j1 + 1, k1 + 1) + &
                        weight(8) * flowdoms(d1, 1, 1)%realCommVars(k + nzippflowcomm)%var(i1 + 1, j1 + 1, k1 + 1)
                end if
            end do
        end do donorloop
 
        ! Now we can do an mpi_gatherv to the root proc:
        call mpi_gatherv(sendbuffer, nvar * comm%nDonor, adflow_real, recvbuffer, &
                         nvar * comm%procSizes, nvar * comm%procDisps, adflow_real, 0, adflow_comm_world, ierr)
        call echk(ierr, __file__, __line__)
 
        call mpi_gatherv(sendbufferd, nvar * comm%nDonor, adflow_real, recvbufferd, &
                         nvar * comm%procSizes, nvar * comm%procDisps, adflow_real, 0, adflow_comm_world, ierr)
        call echk(ierr, __file__, __line__)
 
        deallocate (sendbuffer, sendbufferd)
 
    end subroutine commuserintegrationsurfacevars_d
 
    subroutine commuserintegrationsurfacevars_b(recvBuffer, recvBufferd, varStart, varEnd, comm)
 
        use constants
        use block, onlY: flowdoms, ndom
        use communication, only: myid, adflow_comm_world
        use utils, only: echk
        use oversetutilities, only: fractoweights
 
        implicit none
 
        ! Input/Output
        real(kind=realtype), dimension(:) :: recvbuffer, recvbufferd
        integer(kind=intType), intent(in) :: varStart, varEnd
        type(usersurfcommtype) :: comm
 
        ! Working
        real(kind=realtype), dimension(:), allocatable :: sendbuffer, sendbufferd
        integer(Kind=intType) :: d1, i1, j1, k1, jj, k, nvar, i, ierr
        real(kind=realtype), dimension(8) :: weight
 
        ! The number of variables we are transferring:
        nvar = varend - varstart + 1
 
        allocate (sendbufferd(nvar * comm%nDonor))
 
        ! Adjoint of a gatherv is a scatterv. Flip the send/recv relative
        ! to the forward call.
        call mpi_scatterv(recvbufferd, nvar * comm%procSizes, nvar * comm%procDisps, adflow_real, &
                          sendbufferd, nvar * comm%nDonor, adflow_real, 0, adflow_comm_world, ierr)
        call echk(ierr, __file__, __line__)
 
        ! First generate the interpolated data necessary
        jj = 0
        donorloop: do i = 1, comm%nDonor
            ! Convert the frac to weights
            call fractoweights(comm%frac(:, i), weight)
 
            ! Block and indices for easier reading. The +1 is due to the
            ! pointer offset on realCommVars.
 
            d1 = comm%donorInfo(1, i) ! Block Index
            i1 = comm%donorInfo(2, i) + 1 ! donor I index
            j1 = comm%donorInfo(3, i) + 1 ! donor J index
            k1 = comm%donorInfo(4, i) + 1 ! donor K index
 
            ! We are interpolating nVar variables
            do k = varstart, varend
                jj = jj + 1
                if (d1 > 0) then ! Is this pt valid?
 
                    ! Accumulate back onto the derivative variables.
                    flowdoms(d1, 1, 1)%realCommVars(k + nzippflowcomm)%var(i1, j1, k1) = &
                        flowdoms(d1, 1, 1)%realCommVars(k + nzippflowcomm)%var(i1, j1, k1) + &
                        weight(1) * sendbufferd(jj)
 
                    flowdoms(d1, 1, 1)%realCommVars(k + nzippflowcomm)%var(i1 + 1, j1, k1) = &
                        flowdoms(d1, 1, 1)%realCommVars(k + nzippflowcomm)%var(i1 + 1, j1, k1) + &
                        weight(2) * sendbufferd(jj)
 
                    flowdoms(d1, 1, 1)%realCommVars(k + nzippflowcomm)%var(i1, j1 + 1, k1) = &
                        flowdoms(d1, 1, 1)%realCommVars(k + nzippflowcomm)%var(i1, j1 + 1, k1) + &
                        weight(3) * sendbufferd(jj)
 
                    flowdoms(d1, 1, 1)%realCommVars(k + nzippflowcomm)%var(i1 + 1, j1 + 1, k1) = &
                        flowdoms(d1, 1, 1)%realCommVars(k + nzippflowcomm)%var(i1 + 1, j1 + 1, k1) + &
                        weight(4) * sendbufferd(jj)
 
                    flowdoms(d1, 1, 1)%realCommVars(k + nzippflowcomm)%var(i1, j1, k1 + 1) = &
                        flowdoms(d1, 1, 1)%realCommVars(k + nzippflowcomm)%var(i1, j1, k1 + 1) + &
                        weight(5) * sendbufferd(jj)
 
                    flowdoms(d1, 1, 1)%realCommVars(k + nzippflowcomm)%var(i1 + 1, j1, k1 + 1) = &
                        flowdoms(d1, 1, 1)%realCommVars(k + nzippflowcomm)%var(i1 + 1, j1, k1 + 1) + &
                        weight(6) * sendbufferd(jj)
 
                    flowdoms(d1, 1, 1)%realCommVars(k + nzippflowcomm)%var(i1, j1 + 1, k1 + 1) = &
                        flowdoms(d1, 1, 1)%realCommVars(k + nzippflowcomm)%var(i1, j1 + 1, k1 + 1) + &
                        weight(7) * sendbufferd(jj)
 
                    flowdoms(d1, 1, 1)%realCommVars(k + nzippflowcomm)%var(i1 + 1, j1 + 1, k1 + 1) = &
                        flowdoms(d1, 1, 1)%realCommVars(k + nzippflowcomm)%var(i1 + 1, j1 + 1, k1 + 1) + &
                        weight(8) * sendbufferd(jj)
                end if
            end do
        end do donorloop
 
        deallocate (sendbufferd)
 
    end subroutine commuserintegrationsurfacevars_b
 
    subroutine qsortinterppttype(arr, nn)
 
        use constants
        use block ! Cannot use-only becuase of <= operator
        use utils, only: terminate
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: nn
 
        type(interppttype), dimension(*), intent(inout) :: arr
        !
        !      Local variables.
        !
        integer(kind=intType), parameter :: m = 7
 
        integer(kind=intType) :: nStack
        integer(kind=intType) :: i, j, k, r, l, jStack, ii
 
        integer :: ierr
 
        type(interppttype) :: a, tmp
 
        integer(kind=intType), allocatable, dimension(:) :: stack
        integer(kind=intType), allocatable, dimension(:) :: tmpStack
 
        ! Allocate the memory for stack.
 
        nstack = 100
        allocate (stack(nstack), stat=ierr)
        if (ierr /= 0) &
            call terminate("qsortinterpPt", &
                           "Memory allocation failure for stack")
 
        ! Initialize the variables that control the sorting.
 
        jstack = 0
        l = 1
        r = nn
 
        ! Start of the algorithm
 
        do
 
            ! Check for the size of the subarray.
 
            if ((r - l) < m) then
 
                ! Perform insertion sort
 
                do j = l + 1, r
                    a = arr(j)
                    do i = (j - 1), l, -1
                        if (arr(i) <= a) exit
                        arr(i + 1) = arr(i)
                    end do
                    arr(i + 1) = a
                end do
 
                ! In case there are no more elements on the stack, exit from
                ! the outermost do-loop. Algorithm has finished.
 
                if (jstack == 0) exit
 
                ! Pop stack and begin a new round of partitioning.
 
                r = stack(jstack)
                l = stack(jstack - 1)
                jstack = jstack - 2
 
            else
 
                ! Subarray is larger than the threshold for a linear sort.
                ! Choose median of left, center and right elements as
                ! partitioning element a.
                ! Also rearrange so that (l) <= (l+1) <= (r).
 
                k = (l + r) / 2
                tmp = arr(k)             ! Swap the elements
                arr(k) = arr(l + 1)           ! k and l+1.
                arr(l + 1) = tmp
 
                if (arr(r) < arr(l)) then
                    tmp = arr(l)             ! Swap the elements
                    arr(l) = arr(r)             ! r and l.
                    arr(r) = tmp
                end if
 
                if (arr(r) < arr(l + 1)) then
                    tmp = arr(l + 1)         ! Swap the elements
                    arr(l + 1) = arr(r)           ! r and l+1.
                    arr(r) = tmp
                end if
 
                if (arr(l + 1) < arr(l)) then
                    tmp = arr(l + 1)         ! Swap the elements
                    arr(l + 1) = arr(l)           ! l and l+1.
                    arr(l) = tmp
                end if
 
                ! Initialize the pointers for partitioning.
 
                i = l + 1
                j = r
                a = arr(l + 1)
 
                ! The innermost loop
 
                do
 
                    ! Scan up to find element >= a.
                    do
                        i = i + 1
                        if (a <= arr(i)) exit
                    end do
 
                    ! Scan down to find element <= a.
                    do
                        j = j - 1
                        if (arr(j) <= a) exit
                    end do
 
                    ! Exit the loop in case the pointers i and j crossed.
 
                    if (j < i) exit
 
                    ! Swap the element i and j.
 
                    tmp = arr(i)
                    arr(i) = arr(j)
                    arr(j) = tmp
                end do
 
                ! Swap the entries j and l+1. Remember that a equals
                ! arr(l+1).
 
                arr(l + 1) = arr(j)
                arr(j) = a
 
                ! Push pointers to larger subarray on stack,
                ! process smaller subarray immediately.
 
                jstack = jstack + 2
                if (jstack > nstack) then
 
                    ! Storage of the stack is too small. Reallocate.
 
                    allocate (tmpstack(nstack), stat=ierr)
                    if (ierr /= 0) &
                        call terminate("qsortinterpPt", &
                                       "Memory allocation error for tmpStack")
                    tmpstack = stack
 
                    ! Free the memory of stack, store the old value of nStack
                    ! in tmp and increase nStack.
 
                    deallocate (stack, stat=ierr)
                    if (ierr /= 0) &
                        call terminate("qsortinterpPt", &
                                       "Deallocation error for stack")
                    ii = nstack
                    nstack = nstack + 100
 
                    ! Allocate the memory for stack and copy the old values
                    ! from tmpStack.
 
                    allocate (stack(nstack), stat=ierr)
                    if (ierr /= 0) &
                        call terminate("qsortinterpPt", &
                                       "Memory reallocation error for stack")
                    stack(1:ii) = tmpstack(1:ii)
 
                    ! And finally release the memory of tmpStack.
 
                    deallocate (tmpstack, stat=ierr)
                    if (ierr /= 0) &
                        call terminate("qsortinterpPt", &
                                       "Deallocation error for tmpStack")
                end if
 
                if ((r - i + 1) >= (j - l)) then
                    stack(jstack) = r
                    r = j - 1
                    stack(jstack - 1) = j
                else
                    stack(jstack) = j - 1
                    stack(jstack - 1) = l
                    l = j
                end if
 
            end if
        end do
 
        ! Release the memory of stack.
 
        deallocate (stack, stat=ierr)
        if (ierr /= 0) &
            call terminate("qsortinterpPt", &
                           "Deallocation error for stack")
 
        ! Check in debug mode whether the array is really sorted.
 
        if (debug) then
            do i = 1, (nn - 1)
                if (arr(i + 1) < arr(i)) &
                    call terminate("qsortinterpPt", &
                                   "Array is not sorted correctly")
            end do
        end if
 
    end subroutine qsortinterppttype
 
end module usersurfaceintegrations