oversetUtilities.F90

Go to the documentation of this file.

module oversetutilities
contains
#ifndef USE_TAPENADE
 
    subroutine tic(index)
        use constants
        use oversetdata, only: tstart
        implicit none
        integer(kind=intType), intent(in) :: index
        tstart(index) = mpi_wtime()
 
    end subroutine tic
 
    subroutine toc(index)
        use constants
        use oversetdata, only: tstart, oversettimes
        implicit none
        integer(kind=intType), intent(in) :: index
        oversettimes(index) = oversettimes(index) + mpi_wtime() - tstart(index)
    end subroutine toc
 
    subroutine unwindindex(index, il, jl, kl, i, j, k)
        ! Unwind a 1-based index based on the double halo size of 0:ib,
        ! 0:jb, 0:kb
        use constants
        implicit none
        integer(kind=intType), intent(in) :: index, il, jl, kl
        integer(kind=intType), intent(out) :: i, j, k
        integer(kind=intType) :: ID, isize, jsize, ksize
        id = index - 1
        isize = il + 3
        jsize = jl + 3
        ksize = kl + 3
        i = mod(id, isize)
        j = mod(id / isize, jsize)
        k = id / (isize * jsize)
 
    end subroutine unwindindex
 
    function windindex(i, j, k, il, jl, kl)
 
        use constants
        implicit none
        integer(kind=intType), intent(in) :: i, j, k, il, jl, kl
        integer(kind=intType) :: isize, jsize, ksize, windindex
 
        isize = il + 3
        jsize = jl + 3
        ksize = kl + 3
 
        windindex = k * isize * jsize + j * isize + i + 1
 
    end function windindex
 
    subroutine printoverlapmatrix(overlap)
 
        ! This is a debugging routine to print out the overlap matrix.
        use constants
        use communication, only: myid
        use oversetdata, only: csrmatrix
        implicit none
 
        ! Input/output
        type(csrmatrix), intent(in) :: overlap
 
        ! Working
        integer(kind=intType) :: i, jj
 
        if (myid == 0) then
            ! Now dump out who owns what:
            do i = 1, overlap%nrow
                write (*, "(a,I4, a)", advance='no') 'Row:', i, "   "
                do jj = overlap%rowPtr(i), overlap%rowPtr(i + 1) - 1
                    write (*, "(a,I2, a, es10.5)", advance='no') "(", overlap%colInd(jj), ")", overlap%data(jj)
                end do
                write (*, *) " "
            end do
 
            print *, '--------------------------------------'
            ! Now dump out who owns what:
            do i = 1, overlap%nRow
                write (*, "(a,I4, a)", advance='no') 'Row:', i, "   "
                do jj = overlap%rowPtr(i), overlap%rowPtr(i + 1) - 1
                    write (*, "(a,I2, a, I8)", advance='no') "(", overlap%colInd(jj), ")", int(overlap%assignedProc(jj))
                end do
                write (*, *) " "
            end do
        end if
    end subroutine printoverlapmatrix
 
    subroutine getcumulativeform(sizeArray, n, cumArray)
 
        use constants
        implicit none
 
        ! Input/Output
        integer(kind=intType), intent(in) :: n
        integer(kind=intType), dimension(n), intent(in) :: sizeArray
        integer(kind=intType), dimension(0:n), intent(out) :: cumArray
 
        ! Working
        integer(kind=intType) :: i
 
        cumarray(0) = 0
        do i = 1, n
            cumarray(i) = cumarray(i - 1) + sizearray(i)
        end do
 
    end subroutine getcumulativeform
 
    subroutine transposeoverlap(A, B)
 
        ! Create the matrix Create the matrix transpose.
        ! Inspired by: https://people.sc.fsu.edu/~jburkardt/f_src/sparsekit/sparsekit.f90
        use constants
        use oversetdata, only: csrmatrix
        implicit none
 
        ! Input/Output
        type(csrmatrix), intent(in) :: A
        type(csrmatrix), intent(inout) :: B
 
        ! Working
        integer(kind=intType) :: col, colp1, i, k, next
 
        ! A CSR matrix is the same as a CSC matrix of the transpose. So
        ! essentially the algorithm is convert A as a CSC matrix to B
        ! (a CSR matrix)
 
        ! Allocate space for everything in B
        b%nnz = a%nnz
        b%nRow = a%nCol
        b%nCol = a%nRow
        allocate (b%data(b%nnz), b%colInd(b%nnz), &
                  b%assignedProc(b%nnz), b%rowPtr(b%nRow + 1))
        b%allocated = .true.
        !  Compute lengths of rows of B (ie the columns of A)
 
        b%rowPtr = 0
 
        do i = 1, a%nRow
            do k = a%rowPTr(i), a%rowPtr(i + 1) - 1
                colp1 = a%colInd(k) + 1
                b%rowPtr(colp1) = b%rowPtr(colp1) + 1
 
            end do
        end do
        !
        !  Compute pointers from lengths.
        !
        b%rowPtr(1) = 1
        do i = 1, a%nRow
            b%rowPtr(i + 1) = b%rowPtr(i) + b%rowPtr(i + 1)
        end do
        !
        !  Do the actual copying.
        !
        do i = 1, a%nRow
            do k = a%rowPtr(i), a%rowPtr(i + 1) - 1
                col = a%colInd(k)
                next = b%rowPtr(col)
 
                b%data(next) = a%data(k)
                b%assignedProc(next) = a%assignedProc(k)
 
                b%colInd(next) = i
                b%rowPtr(col) = next + 1
            end do
        end do
 
        !  Reshift B%rowPtr
        do i = a%nRow, 1, -1
            b%rowPtr(i + 1) = b%rowPtr(i)
        end do
        b%rowPtr(1) = 1
 
    end subroutine transposeoverlap
 
    subroutine deallocatecsrmatrix(mat1)
 
        use constants
        use oversetdata, only: csrmatrix
        implicit none
 
        type(csrmatrix), intent(inout) :: mat1
 
        if (mat1%allocated) then
            deallocate ( &
                mat1%data, &
                mat1%colInd, &
                mat1%rowPtr, &
                mat1%assignedProc)
        end if
 
    end subroutine deallocatecsrmatrix
 
    subroutine computefringeprocarray(fringes, n, fringeProc, cumFringeProc, nFringeProc)
        ! Compute the breakpoints "cumFringeProc" for a list of sorted n
        ! fringes "fringes". nFringeProc is the total number of unique
        ! processors. fringeProc is the processor number for each section.
        use constants
        use block, only: fringetype
        use oversetdata, only: csrmatrix
        use communication, only: nproc
 
        implicit none
 
        ! Input/Output
        integer(kind=intType), intent(in) :: n
        type(fringetype), intent(in), dimension(n) :: fringes
        integer(kind=intType), intent(out) :: nFringeProc
        integer(kind=intType), intent(out) :: fringeProc(nProc), cumFringeProc(1:nProc + 1)
 
        ! Working
        integer(kind=intType) :: i, currentProc
 
        fringeproc = -1
        nfringeproc = 0
        cumfringeproc(1) = 1
        currentproc = -1
 
        do i = 1, n
            if (currentproc /= fringes(i)%donorProc) then
                nfringeproc = nfringeproc + 1
                cumfringeproc(nfringeproc) = i
                fringeproc(nfringeproc) = fringes(i)%donorProc
                currentproc = fringes(i)%donorProc
            end if
        end do
 
        ! Finally, the nFringeProc+1 entry is always n+1
        cumfringeproc(nfringeproc + 1) = n + 1
 
    end subroutine computefringeprocarray
 
    subroutine deallocateoblocks(oBlocks, n)
 
        ! This subroutine deallocates all data stores in a list of oBlocks
        use constants
        use adtbuild, only: destroyserialhex
        use oversetdata, only: oversetblock
        implicit none
 
        ! Input Params
        integer(kind=intType) :: n
        type(oversetblock), dimension(n), intent(inout) :: oBLocks
 
        ! Working Parameters
        integer(kind=intType) :: i
 
        do i = 1, n
 
            ! oBlock:
            if (oblocks(i)%allocated) then
                deallocate ( &
                    oblocks(i)%hexaConn, &
                    oblocks(i)%globalCell, &
                    oblocks(i)%nearWall, &
                    oblocks(i)%invalidDonor, &
                    oblocks(i)%qualDonor, &
                    oblocks(i)%xADT)
                if (allocated(oblocks(i)%rbuffer)) then
                    deallocate (oblocks(i)%rBuffer, &
                                oblocks(i)%iBuffer)
                end if
                call destroyserialhex(oblocks(i)%ADT)
            end if
        end do
    end subroutine deallocateoblocks
 
    subroutine deallocateofringes(oFringes, n)
 
        ! This subroutine deallocates all data stores in a list of oFringes
        use constants
        use oversetdata, only: oversetfringe
        implicit none
 
        ! Input Params
        integer(kind=intType) :: n
        type(oversetfringe), dimension(n), intent(inout) :: oFringes
 
        ! Working Parameters
        integer(kind=intType) :: i
 
        do i = 1, n
            if (allocated(ofringes(i)%x)) &
                deallocate (ofringes(i)%x)
 
            if (allocated(ofringes(i)%isWall)) &
                deallocate (ofringes(i)%isWall)
 
            if (allocated(ofringes(i)%xSeed)) &
                deallocate (ofringes(i)%xSeed)
 
            if (allocated(ofringes(i)%wallInd)) &
                deallocate (ofringes(i)%wallInd)
 
            if (allocated(ofringes(i)%rbuffer)) &
                deallocate (ofringes(i)%rBuffer)
 
            if (allocated(ofringes(i)%ibuffer)) &
                deallocate (ofringes(i)%iBuffer)
 
            if (associated(ofringes(i)%fringeIntBuffer)) &
                deallocate (ofringes(i)%fringeIntBuffer)
 
            if (associated(ofringes(i)%fringeRealBuffer)) &
                deallocate (ofringes(i)%fringeRealBuffer)
 
            ofringes(i)%allocated = .false.
        end do
    end subroutine deallocateofringes
 
    subroutine deallocateosurfs(oSurfs, n)
 
        ! This subroutine deallocates all data stores in a list of oSurfs
 
        use constants
        use adtbuild, only: destroyserialquad
        use oversetdata, only: oversetwall
        use kdtree2_module, only: kdtree2destroy
        implicit none
 
        ! Input Params
        integer(kind=intType) :: n
        type(oversetwall), dimension(n), intent(inout) :: oSurfs
 
        ! Working Parameters
        integer(kind=intType) :: i
 
        do i = 1, n
            if (osurfs(i)%allocated) then
                deallocate ( &
                    osurfs(i)%x, &
                    osurfs(i)%conn, &
                    osurfs(i)%iblank, &
                    osurfs(i)%cellPtr)
                call destroyserialquad(osurfs(i)%ADT)
                if (osurfs(i)%nNodes > 0) then
                    call kdtree2destroy(osurfs(i)%tree)
                end if
            end if
            osurfs(i)%allocated = .false.
        end do
    end subroutine deallocateosurfs
 
    subroutine wallsonblock(wallsPresent)
 
        use constants
        use blockpointers, only: nbkglobal
        use cgnsgrid, only: cgnsdoms
        implicit none
 
        logical, intent(out) :: wallsPresent
        integer(kind=intType) :: mm
        wallspresent = .false.
        ! Check THE ORIGINAL CGNS blocks for BCs, because the block may have
        ! been split.
        do mm = 1, cgnsdoms(nbkglobal)%nBocos
            if (cgnsdoms(nbkglobal)%bocoInfo(mm)%BCType == nswalladiabatic .or. &
                cgnsdoms(nbkglobal)%bocoInfo(mm)%BCType == nswallisothermal .or. &
                cgnsdoms(nbkglobal)%bocoInfo(mm)%BCType == eulerwall) then
                wallspresent = .true.
            end if
        end do
    end subroutine wallsonblock
 
    subroutine flagforcedrecv
 
        use constants
        use blockpointers, only: nx, ny, nz, ie, je, ke, bcdata, bcfaceid, nbocos, bctype, &
                                 forcedrecv, flowdoms, ndom, il, jl, kl, iblank, status
        use utils, only: setpointers
        use communication
        use haloexchange, only: whalo1to1intgeneric, whalo1to1intgeneric_b
        use stencils
        implicit none
 
        ! This is generic routine for filling up a 3D array of 1st level halos
        ! cells (1:ie, 1:je, 1:ke) indicating cells that are forced
        ! receivers. BlockPointers must have already been set.
 
        integer(kind=intType) :: nn, i, j, k, mm, iStart, iEnd, jStart, jEnd, kStart, kEnd
        integer(kind=intType) :: ii, jj, kk, i_stencil
        logical :: floodOrBlank, floodOrBlank2
        do nn = 1, ndom
            call setpointers(nn, 1, 1)
            forcedrecv = 0
            do mm = 1, nbocos
                ! Just record the ranges necessarvy and we'll add in a generic
                ! loop. Why is it the first three? Well, the first level of halos
                ! off of an overset outer bound is completely
                ! meaningless. Essentially we ignore those. So the outer two
                ! layers of cells are indices 2 and 3. Therefore the first 3 on
                ! either side need to be flagged as invalid.
 
                select case (bcfaceid(mm))
                case (imin)
                    istart = 1; iend = 3;
                    jstart = bcdata(mm)%inBeg + 1; jend = bcdata(mm)%inEnd
                    kstart = bcdata(mm)%jnBeg + 1; kend = bcdata(mm)%jnEnd
                case (imax)
                    istart = nx; iend = ie;
                    jstart = bcdata(mm)%inBeg + 1; jend = bcdata(mm)%inEnd
                    kstart = bcdata(mm)%jnBeg + 1; kend = bcdata(mm)%jnEnd
                case (jmin)
                    istart = bcdata(mm)%inBeg + 1; iend = bcdata(mm)%inEnd
                    jstart = 1; jend = 3;
                    kstart = bcdata(mm)%jnBeg + 1; kend = bcdata(mm)%jnEnd
                case (jmax)
                    istart = bcdata(mm)%inBeg + 1; iend = bcdata(mm)%inEnd
                    jstart = ny; jend = je;
                    kstart = bcdata(mm)%jnBeg + 1; kend = bcdata(mm)%jnEnd
                case (kmin)
                    istart = bcdata(mm)%inBeg + 1; iend = bcdata(mm)%inEnd
                    jstart = bcdata(mm)%jnBeg + 1; jend = bcdata(mm)%jnEnd
                    kstart = 1; kend = 3;
                case (kmax)
                    istart = bcdata(mm)%inBeg + 1; iend = bcdata(mm)%inEnd
                    jstart = bcdata(mm)%jnBeg + 1; jend = bcdata(mm)%jnEnd
                    kstart = nz; kend = ke;
                end select
 
                if (bctype(mm) == oversetouterbound) then
                    do k = kstart, kend
                        do j = jstart, jend
                            do i = istart, iend
                                forcedrecv(i, j, k) = 1
                            end do
                        end do
                    end do
                end if
            end do
 
            ! Add to the invalid donor list if it got flooded with iblank of -2 or -3:
            do k = 2, kl
                do j = 2, jl
                    do i = 2, il
                        ! Flooded or explictly blanked cell
                        floodorblank = isflooded(status(i, j, k)) .or. &
                                       isfloodseed(status(i, j, k)) .or. &
                                       iblank(i, j, k) == -4
                        if (floodorblank) then
                            stencilloop: do i_stencil = 1, n_visc_drdw
                                ii = visc_drdw_stencil(i_stencil, 1) + i
                                jj = visc_drdw_stencil(i_stencil, 2) + j
                                kk = visc_drdw_stencil(i_stencil, 3) + k
                                ! Flag as a forced reciver if it *isn't* flooded
                                ! or explictly blanked
 
                                floodorblank2 = isflooded(status(ii, jj, kk)) .or. &
                                                isfloodseed(status(ii, jj, kk)) .or. &
                                                iblank(ii, jj, kk) == -4
 
                                if (.not. floodorblank2) then
                                    forcedrecv(ii, jj, kk) = 1
                                end if
                            end do stencilloop
                        end if
                    end do
                end do
            end do
        end do
 
        ! Update the info across block boundaries
        domainloop: do nn = 1, ndom
            flowdoms(nn, 1, 1)%intCommVars(1)%var => &
                flowdoms(nn, 1, 1)%forcedRecv(:, :, :)
        end do domainloop
 
        ! Run the reverse halo exchange first. This is necessary if there
        ! is 1 cell wide block next to a overset outer bound like the following:
        !
        !         Blk1         Blk2
        !   ----+-----+------++-----+
        !       |     |      ||     |    <= this face has overset outer bound BC
        !   ----+-----+------++-----+
        !       |     |      ||     |
        !   ----+-----+------++-----+
        !       |     |      ||     |
        !   ----+-----+------++-----+
        !                    ^block boundary
        !
        ! So what happens, is blk2 (1 cell wide) sets the two layers of
        ! cells off of the BC as forced receivers. However, the second of
        ! those layers is a halo cell. Blk1 then never gets this
        ! information. as it clearly should. So what we have to do is a
        ! reverse halo exchange that takes halo information and combines
        ! it with real cell information. Essentially we will accumulate
        ! forcedRecv from the halo to the real cell. For this we run the
        ! generic reverse halo exchange.
 
        call whalo1to1intgeneric_b(1, 1, 1, commpatterncell_2nd, internalcell_2nd)
 
        ! Finally we now need to run the forward halo exchange to make
        ! sure any halos on other procs are set correctly that may be part of a stencil
        call whalo1to1intgeneric(1, 1, 1, commpatterncell_2nd, internalcell_2nd)
 
    end subroutine flagforcedrecv
 
    ! Utility function for unpacking/accessing the status variable
 
    function isdonor(i)
        use constants
        implicit none
        logical :: isdonor, ishole, iscompute, isfloodseed, isflooded, iswalldonor, isreceiver
        integer(kind=intType), intent(in) :: i
        call getstatus(i, isdonor, isdonor, iscompute, isfloodseed, isflooded, iswalldonor, isreceiver)
    end function isdonor
 
    function ishole(i)
        use constants
        implicit none
        logical :: isdonor, ishole, iscompute, isfloodseed, isflooded, iswalldonor, isreceiver
        integer(kind=intType), intent(in) :: i
        call getstatus(i, isdonor, ishole, iscompute, isfloodseed, isflooded, iswalldonor, isreceiver)
    end function ishole
 
    function iscompute(i)
        use constants
        implicit none
        logical :: isdonor, ishole, iscompute, isfloodseed, isflooded, iswalldonor, isreceiver
        integer(kind=intType), intent(in) :: i
        call getstatus(i, isdonor, ishole, iscompute, isfloodseed, isflooded, iswalldonor, isreceiver)
    end function iscompute
 
    function isfloodseed(i)
        use constants
        implicit none
        logical :: isdonor, ishole, iscompute, isfloodseed, isflooded, iswalldonor, isreceiver
        integer(kind=intType), intent(in) :: i
        call getstatus(i, isdonor, ishole, iscompute, isfloodseed, isflooded, iswalldonor, isreceiver)
    end function isfloodseed
 
    function isflooded(i)
        use constants
        implicit none
        logical :: isdonor, ishole, iscompute, isfloodseed, isflooded, iswalldonor, isreceiver
        integer(kind=intType), intent(in) :: i
        call getstatus(i, isdonor, ishole, iscompute, isfloodseed, isflooded, iswalldonor, isreceiver)
    end function isflooded
 
    function iswalldonor(i)
        use constants
        implicit none
        logical :: isdonor, ishole, iscompute, isfloodseed, isflooded, iswalldonor, isreceiver
        integer(kind=intType), intent(in) :: i
        call getstatus(i, isdonor, ishole, iscompute, isfloodseed, isflooded, iswalldonor, isreceiver)
    end function iswalldonor
 
    function isreceiver(i)
        use constants
        implicit none
        logical :: isdonor, ishole, iscompute, isfloodseed, isflooded, iswalldonor, isreceiver
        integer(kind=intType), intent(in) :: i
        call getstatus(i, isdonor, ishole, iscompute, isfloodseed, isflooded, iswalldonor, isreceiver)
    end function isreceiver
 
    subroutine setisdonor(i, flag)
        use constants
        implicit none
        integer(kind=intType), intent(inout) :: i
        logical :: isDonor, isHole, isCompute, isFloodSeed, isFlooded, isWallDonor, isReceiver, flag
        call getstatus(i, isdonor, ishole, iscompute, isfloodseed, isflooded, iswalldonor, isreceiver)
        call setstatus(i, flag, ishole, iscompute, isfloodseed, isflooded, iswalldonor, isreceiver)
    end subroutine setisdonor
 
    subroutine setishole(i, flag)
        use constants
        implicit none
        integer(kind=intType), intent(inout) :: i
        logical :: isDonor, isHole, isCompute, isFloodSeed, isFlooded, isWallDonor, isReceiver, flag
        call getstatus(i, isdonor, ishole, iscompute, isfloodseed, isflooded, iswalldonor, isreceiver)
        call setstatus(i, isdonor, flag, iscompute, isfloodseed, isflooded, iswalldonor, isreceiver)
    end subroutine setishole
 
    subroutine setiscompute(i, flag)
        use constants
        implicit none
        integer(kind=intType), intent(inout) :: i
        logical :: isDonor, isHole, isCompute, isFloodSeed, isFlooded, isWallDonor, isReceiver, flag
        call getstatus(i, isdonor, ishole, iscompute, isfloodseed, isflooded, iswalldonor, isreceiver)
        call setstatus(i, isdonor, ishole, flag, isfloodseed, isflooded, iswalldonor, isreceiver)
    end subroutine setiscompute
 
    subroutine setisfloodseed(i, flag)
        use constants
        implicit none
        integer(kind=intType), intent(inout) :: i
        logical :: isDonor, isHole, isCompute, isFloodSeed, isFlooded, isWallDonor, isReceiver, flag
        call getstatus(i, isdonor, ishole, iscompute, isfloodseed, isflooded, iswalldonor, isreceiver)
        call setstatus(i, isdonor, ishole, iscompute, flag, isflooded, iswalldonor, isreceiver)
    end subroutine setisfloodseed
 
    subroutine setisflooded(i, flag)
        use constants
        implicit none
        integer(kind=intType), intent(inout) :: i
        logical :: isDonor, isHole, isCompute, isFloodSeed, isFlooded, isWallDonor, isReceiver, flag
        call getstatus(i, isdonor, ishole, iscompute, isfloodseed, isflooded, iswalldonor, isreceiver)
        call setstatus(i, isdonor, ishole, iscompute, isfloodseed, flag, iswalldonor, isreceiver)
    end subroutine setisflooded
 
    subroutine setiswalldonor(i, flag)
        use constants
        implicit none
        integer(kind=intType), intent(inout) :: i
        logical :: isDonor, isHole, isCompute, isFloodSeed, isFlooded, isWallDonor, isReceiver, flag
        call getstatus(i, isdonor, ishole, iscompute, isfloodseed, isflooded, iswalldonor, isreceiver)
        call setstatus(i, isdonor, ishole, iscompute, isfloodseed, isflooded, flag, isreceiver)
    end subroutine setiswalldonor
 
    subroutine setisreceiver(i, flag)
        use constants
        implicit none
        integer(kind=intType), intent(inout) :: i
        logical :: isDonor, isHole, isCompute, isFloodSeed, isFlooded, isWallDonor, isReceiver, flag
        call getstatus(i, isdonor, ishole, iscompute, isfloodseed, isflooded, iswalldonor, isreceiver)
        call setstatus(i, isdonor, ishole, iscompute, isfloodseed, isflooded, iswalldonor, flag)
    end subroutine setisreceiver
 
    subroutine setstatus(i, isDonor, isHole, isCompute, isFloodSeed, isFlooded, isWallDonor, isReceiver)
 
        use constants
        implicit none
        integer(kind=intType), intent(out) :: i
        logical :: isDonor, isHole, isCompute, isFloodSeed, isFlooded, isWallDonor, isReceiver
        i = 0
 
        if (isdonor) i = i + 1
        if (ishole) i = i + 2
        if (iscompute) i = i + 4
        if (isfloodseed) i = i + 8
        if (isflooded) i = i + 16
        if (iswalldonor) i = i + 32
        if (isreceiver) i = i + 64
    end subroutine setstatus
 
    subroutine getstatus(i, isDonor, isHole, isCompute, isFloodSeed, isFlooded, isWallDonor, isReceiver)
 
        use constants
        implicit none
        logical :: isDonor, isHole, isCompute, isFloodSeed, isFlooded, isWallDonor, isReceiver
        integer(kind=intType) :: i, j
        j = i
 
        isdonor = .false.
        ishole = .false.
        iscompute = .false.
        isfloodseed = .false.
        isflooded = .false.
        iswalldonor = .false.
        isreceiver = .false.
 
        if (j / 64 > 0) then
            isreceiver = .true.
            j = j - 64
        end if
 
        if (j / 32 > 0) then
            iswalldonor = .true.
            j = j - 32
        end if
 
        if (j / 16 > 0) then
            isflooded = .true.
            j = j - 16
        end if
 
        if (j / 8 > 0) then
            isfloodseed = .true.
            j = j - 8
        end if
 
        if (j / 4 > 0) then
            iscompute = .true.
            j = j - 4
        end if
 
        if (j / 2 > 0) then
            ishole = .true.
            j = j - 2
        end if
 
        if (j / 1 > 0) then
            isdonor = .true.
            j = j - 1
        end if
    end subroutine getstatus
 
    !
    subroutine binsearchnodes(arr, searchNode, nn, searchInd)
 
        !  binSearchNodes does binary search for a node 'searchNode'
        !  in arr(1:nn) and returns index 'searchInd' where
        !  'searchNode' lies in arr. searchInd = -1 if not found.
 
        use constants
        implicit none
 
        ! Input parameters
        integer(kind=intType), intent(in) :: nn, searchNode
        integer(kind=intType), intent(out) :: searchInd
        integer(kind=intType), intent(in) :: arr(nn)
 
        ! Local variables
        integer(kind=intType) :: first, last, middle
 
        first = 1
        last = nn
 
        middle = (first + last) / 2
 
        do while (first <= last)
            if (arr(middle) < searchnode) then
                first = middle + 1
            else if (arr(middle) == searchnode) then
                searchind = middle
                exit
            else
                last = middle - 1
            end if
 
            middle = (first + last) / 2
        end do !while
 
        if (first > last) then
            searchind = -1
            print *, ' binSearchNode fails for searchNode ', searchnode
            stop
        end if
    end subroutine binsearchnodes
 
    subroutine binsearchpocketedgetype(arr, search, nn, searchInd)
 
        !  binSearchPocketEdgeType does binary searche for
        !  pocketEdgeType 'search' edge and returns index 'searchInd'
        !  where 'search' lies in arr.
 
        use constants
        use oversetdata ! cannot use only becuase of <= operator
        implicit none
 
        ! Input parameters
        integer(kind=intType), intent(in) :: nn
        integer(kind=intType), intent(out) :: searchInd
 
        type(pocketedge), intent(in) :: search
        type(pocketedge), dimension(*), intent(in) :: arr
 
        ! Local variables
        integer(kind=intType) :: first, last, middle
 
        first = 1
        last = nn
 
        middle = (first + last) / 2
 
        do while (first <= last)
            if (arr(middle) < search) then
                first = middle + 1
            else if (arr(middle) == search) then
                searchind = middle
                exit
            else
                last = middle - 1
            end if
 
            middle = (first + last) / 2
        end do !while
 
        if (first > last) then
            print *, ' binSearchPocketEdgeType fails for Edge with nodes ', &
                search%n1, search%n2
            stop
        end if
    end subroutine binsearchpocketedgetype
 
    !
    subroutine qsortedgetype(arr, nn)
        !
        !       qsortEdgeType sorts the given number of oversetString master
        !       Edges in increasing order based on the <= operator for this
        !       derived data type.
        !       (Generously copied from qsortFringeType.F90)
        !
        use constants
        use oversetdata ! cannot use only becuase of <= operator
        use utils, only: terminate
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: nn
 
        type(oversetedge), 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(oversetedge) :: 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("qsortEdgeType", &
                           "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("qsortEdgeType", &
                                       "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("qsortEdgeType", &
                                       "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("qsortEdgeType", &
                                       "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("qsortEdgeType", &
                                       "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("qsortEdgeType", &
                           "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("qsortEdgeType", &
                                   "Array is not sorted correctly")
            end do
        end if
 
    end subroutine qsortedgetype
 
    subroutine qsortfringetype(arr, nn, sortType)
        !
        !       qsortFringeListTy sorts the given number of fringes
        !       increasing order based on the <= operator for this derived
        !       data type.
        !
        use constants
        use block ! Cannot use-only becuase of <= operator
        use utils, only: terminate
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: nn, sortType
 
        type(fringetype), 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(fringetype) :: a, tmp
 
        integer(kind=intType), allocatable, dimension(:) :: stack
        integer(kind=intType), allocatable, dimension(:) :: tmpStack
 
        if (sorttype == sortbydonor) then
            fringesorttype = sortbydonor
        else if (sorttype == sortbyreceiver) then
            fringesorttype = sortbyreceiver
        else
            call terminate("qsortfringeType", &
                           "Uuknown sort type")
        end if
 
        ! Allocate the memory for stack.
 
        nstack = 100
        allocate (stack(nstack), stat=ierr)
        if (ierr /= 0) &
            call terminate("qsortfringeType", &
                           "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("qsortfringeType", &
                                       "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("qsortfringeType", &
                                       "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("qsortfringeType", &
                                       "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("qsortfringeType", &
                                       "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("qsortfringeType", &
                           "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("qsortfringeType", &
                                   "Array is not sorted correctly")
            end do
        end if
 
    end subroutine qsortfringetype
 
    subroutine addtofringelist(fringeList, n, fringe)
 
        ! Generic subroutine to add to a fringe "List". This isn't
        ! actually a list but rather an allocatable (pointer) array that
        ! is periodically resized as necessary. n is the current size
        ! which is automatically incremented by this routine. This
        ! operation occurs multiple times throughout the overset code.
 
        use constants
        use block, only: fringetype
 
        implicit none
 
        ! Input Params
        type(fringetype), dimension(:), pointer :: fringeList
        type(fringetype) :: fringe
        integer(kind=intType), intent(inout) :: n
 
        ! Working Paramters
        integer(kind=intType) :: fSize
        type(fringetype), dimension(:), pointer :: tmpFringePtr
        fsize = size(fringelist)
 
        ! Increment n for next item
        n = n + 1
 
        if (n > fsize) then
 
            ! Pointer to existing data:
            tmpfringeptr => fringelist
 
            ! Allocate new space
            allocate (fringelist(int(1.5 * fsize)))
 
            ! Copy exsitng values
            fringelist(1:fsize) = tmpfringeptr(1:fsize)
 
            ! Free original memory
            deallocate (tmpfringeptr)
 
        end if
 
        fringelist(n) = fringe
 
    end subroutine addtofringelist
 
    subroutine addtofringebuffer(intBuffer, realBuffer, n, fringe)
 
        ! Generic subroutine to add to a fringe "List". It isn't actually
        ! a list of fringe types but rather a real array and an int array.
 
        use constants
        use block, only: fringetype
 
        implicit none
 
        ! Input Params
        integer(kind=intType), dimension(:, :), pointer :: intBuffer
        real(kind=realtype), dimension(:, :), pointer :: realbuffer
        type(fringetype) :: fringe
        integer(kind=intType), intent(inout) :: n
 
        ! Working Paramters
        integer(kind=intType) :: fsize
        integer(kind=intType), dimension(:, :), pointer :: tmpint
        real(kind=realtype), dimension(:, :), pointer :: tmpreal
        fsize = size(intbuffer, 2)
 
        ! Increment n for next item
        n = n + 1
 
        if (n > fsize) then
 
            ! Pointers to existing data:
            tmpint => intbuffer
            tmpreal => realbuffer
 
            ! Allocate new space
            allocate (intbuffer(5, int(1.5 * fsize)))
            allocate (realbuffer(4, int(1.5 * fsize)))
 
            ! Copy exsitng values
            intbuffer(:, 1:fsize) = tmpint(:, 1:fsize)
            realbuffer(:, 1:fsize) = tmpreal(:, 1:fsize)
 
            ! Free original memory
            deallocate (tmpint, tmpreal)
 
        end if
 
        ! Now we can safely add the information:
        intbuffer(1, n) = fringe%donorProc
        intbuffer(2, n) = fringe%donorBlock
        intbuffer(3, n) = fringe%dIndex
        intbuffer(4, n) = fringe%myBlock
        intbuffer(5, n) = fringe%myIndex
 
        realbuffer(1:3, n) = fringe%donorFrac
        realbuffer(4, n) = fringe%quality
 
    end subroutine addtofringebuffer
 
    subroutine qsortpocketedgetype(arr, nn)
        !
        !       qsortPocketEdgeType sorts the given number of oversetString
        !       master Edges in increasing order based on the <= operator for
        !       this derived data type.
        !       (Generously copied from qsortFringeType.F90)
        !
        use constants
        use oversetdata ! Cannot use-only becuase of <= operator
        use utils, onlY: terminate
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: nn
 
        type(pocketedge), 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(pocketedge) :: 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("qsortEdgeType", &
                           "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("qsortEdgeType", &
                                       "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("qsortEdgeType", &
                                       "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("qsortEdgeType", &
                                       "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("qsortEdgeType", &
                                       "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("qsortEdgeType", &
                           "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("qsortEdgeType", &
                                   "Array is not sorted correctly")
            end do
        end if
 
    end subroutine qsortpocketedgetype
 
    subroutine checkoverset(level, sps, totalOrphans, lastCall)
 
        !
        !       CheckOverset checks the integrity of the overset connectivity
        !       and holes. For every comptue cell (iblank = 1) it checks that
        !       every cell in its stencil are not blanked. If even 1 cell is
        !      * found with an incomplete stencil it is a fatal error.
        use constants
        use blockpointers, only: il, jl, kl, iblank, flowdoms, ndom, orphans, &
                                 ibegor, jbegor, kbegor, nbkglobal, norphans
        use stencils, only: visc_drdw_stencil, n_visc_drdw
        use communication, only: myid, adflow_comm_world
        use utils, only: setpointers, echk
        use inputoverset, only: oversetdebugprint
        implicit none
 
        ! Input/Output
        integer(kind=intType), intent(in) :: level, sps
        integer(kind=intType), intent(out) :: totalOrphans
        logical, intent(in) :: lastCall
 
        ! Working
        integer(kind=intType) :: i, j, k, nn, ii, jj, kk, n, ierr
        integer(kind=intType) :: magic, localOrphans, i_stencil
        logical :: badCell
 
        ! This pass just lets the user know where the bad cells are.
        do nn = 1, ndom
            call setpointers(nn, level, sps)
 
            ! On the first pass count up the total number of orphans for this block
            n = 0
            do k = 2, kl
                do j = 2, jl
                    do i = 2, il
                        if (iblank(i, j, k) == 1) then
                            badcell = .false.
 
                            stencilloop: do i_stencil = 1, n_visc_drdw
                                ii = visc_drdw_stencil(i_stencil, 1) + i
                                jj = visc_drdw_stencil(i_stencil, 2) + j
                                kk = visc_drdw_stencil(i_stencil, 3) + k
 
                                if (.not. (iblank(ii, jj, kk) == 1 .or. iblank(ii, jj, kk) == -1)) then
                                    badcell = .true.
                                end if
                            end do stencilloop
                            if (badcell .and. lastcall .and. oversetdebugprint) then
                                print *, 'Error in connectivity at :', nbkglobal, i + ibegor, j + jbegor, k + kbegor
                            end if
                        end if
                    end do
                end do
            end do
        end do
 
        magic = 33
        localorphans = 0
        do nn = 1, ndom
            call setpointers(nn, level, sps)
 
            ! On the first pass count up the total number of orphans for this block
            n = 0
            do k = 2, kl
                do j = 2, jl
                    do i = 2, il
                        if (iblank(i, j, k) == 0 .or. iblank(i, j, k) == -2 .or. &
                            iblank(i, j, k) == -3 .or. iblank(i, j, k) == -4) then
 
                            stencilloop2: do i_stencil = 1, n_visc_drdw
                                ii = visc_drdw_stencil(i_stencil, 1) + i
                                jj = visc_drdw_stencil(i_stencil, 2) + j
                                kk = visc_drdw_stencil(i_stencil, 3) + k
 
                                if (ii >= 2 .and. jj >= 2 .and. kk >= 2 .and. &
                                    ii <= il .and. jj <= jl .and. kk <= kl) then
 
                                    if (iblank(ii, jj, kk) == 1) then
                                        ! This cell is an orphan:
                                        n = n + 1
                                    end if
                                end if
                            end do stencilloop2
                        end if
                    end do
                end do
            end do
 
            localorphans = localorphans + n
 
            ! Remove any existing orphans
            if (associated(flowdoms(nn, level, sps)%orphans)) then
                deallocate (flowdoms(nn, level, sps)%orphans)
            end if
            allocate (flowdoms(nn, level, sps)%orphans(3, n))
 
            ! Save the total number of orphans on this block
            flowdoms(nn, level, sps)%nOrphans = n
 
            ! Manual set information from blockPointers that would be set
            ! with setPointers()
            orphans => flowdoms(nn, level, sps)%orphans
            norphans = n
 
            ! On the first pass count up the total number of orphans for this block
            n = 0
            do k = 2, kl
                do j = 2, jl
                    do i = 2, il
                        if (iblank(i, j, k) == 0 .or. iblank(i, j, k) == -2 .or. &
                            iblank(i, j, k) == -3 .or. iblank(i, j, k) == -4) then
 
                            stencilloop3: do i_stencil = 1, n_visc_drdw
                                ii = visc_drdw_stencil(i_stencil, 1) + i
                                jj = visc_drdw_stencil(i_stencil, 2) + j
                                kk = visc_drdw_stencil(i_stencil, 3) + k
 
                                if (ii >= 2 .and. jj >= 2 .and. kk >= 2 .and. &
                                    ii <= il .and. jj <= jl .and. kk <= kl) then
 
                                    if (iblank(ii, jj, kk) == 1) then
                                        ! This cell is an orphan:
                                        n = n + 1
                                        orphans(:, n) = (/ii, jj, kk/)
                                        if (lastcall) then
                                            ! we can modify iBlankLast because this is the last checkOverset call.
                                            ! we set iBlankLast to -5 to mark orphan cells, this value will then
                                            ! be moved to iBlank after we are done with other loops.
                                            flowdoms(nn, level, sps)%iBlankLast(ii, jj, kk) = -5
                                        end if
                                    end if
                                end if
                            end do stencilloop3
                        end if
                    end do
                end do
            end do
        end do
 
        ! Determine the total number of overset orphans
        call mpi_allreduce(localorphans, totalorphans, 1, adflow_integer, mpi_sum, &
                           adflow_comm_world, ierr)
        call echk(ierr, __file__, __line__)
 
        if (myid == 0) then
            print *, 'Total number of orphans:', totalorphans
        end if
 
        ! if this is the last checkOverset call with lastCall = .True., then
        ! we can move the orphan information to iBlank because we have done all
        ! the checking required.
        if (lastcall .and. (totalorphans .gt. 0)) then
            ! loop over the cells and write the orphan information to iBlank
            do nn = 1, ndom
                call setpointers(nn, level, sps)
                do k = 2, kl
                    do j = 2, jl
                        do i = 2, il
                            if (flowdoms(nn, level, sps)%iblankLast(i, j, k) == -5) iblank(i, j, k) = -5
                        end do
                    end do
                end do
            end do
        end if
 
    end subroutine checkoverset
 
    !
    !       Finds quadratic uvw weights [-1:1] for interpolant point
    !       xSearch by solving for its co-ord in donor element defined by
    !       xElem using newton-raphson iteration.
    !       Fractions uvwQuadratic come with initial guess from ADT search
    !       used for linear interpolation.
 
    subroutine computequadraticweights(xSearch, xElem, uvwQuadratic)
 
        use constants
 
        implicit none
 
        ! Input variables
        real(kind=realtype), intent(in) :: xsearch(3), xelem(3, -1:1, -1:1, -1:1)
        real(kind=realtype), intent(inout) :: uvwquadratic(3)
 
        ! Working variables
        ! -----------------------------------------------------------------
        ! newton related
        integer(kind=intType) :: n, niter
        real(kind=realtype) :: resid, residtol
        real(kind=realtype) :: b(3, 4)
 
        ! others
        integer(kind=intType) :: j, l, iii, jjj, kkk
        real(kind=realtype) :: ff(27), shp(3, 3), psi(3)
        real(kind=realtype) :: dff(27, 3), dshp(3, 3, 3) !-> differentials wrt 3 psi dirs
        logical :: ok_flag
        ! -----------------------------------------------------------------
 
        ! Initialize newton parameters
        niter = 0
        resid = 1.0e10
        residtol = 1.0e-15
 
        ! Begin newton iterations to solve for uvw weights wrt quadratic element
 
        newton_loop: do while (niter < 10 .and. resid > residtol)
 
            !step 1: find weights
            !--------------------
 
            ! Initialize weights
            psi(1:3) = uvwquadratic(1:3)
 
            ! Precopute the FE shape functions for each j-direction
            do j = 1, 3
                shp(1, j) = half * psi(j) * (psi(j) - one)
                shp(2, j) = -(psi(j)**2 - 1)
                shp(3, j) = half * psi(j) * (psi(j) + one)
            end do
 
            ! These are the 27 quadratic weights
            ff(1) = shp(1, 1) * shp(1, 2) * shp(1, 3)
            ff(2) = shp(2, 1) * shp(1, 2) * shp(1, 3)
            ff(3) = shp(3, 1) * shp(1, 2) * shp(1, 3)
 
            ff(4) = shp(1, 1) * shp(2, 2) * shp(1, 3)
            ff(5) = shp(2, 1) * shp(2, 2) * shp(1, 3)
            ff(6) = shp(3, 1) * shp(2, 2) * shp(1, 3)
 
            ff(7) = shp(1, 1) * shp(3, 2) * shp(1, 3)
            ff(8) = shp(2, 1) * shp(3, 2) * shp(1, 3)
            ff(9) = shp(3, 1) * shp(3, 2) * shp(1, 3)
 
            ff(10) = shp(1, 1) * shp(1, 2) * shp(2, 3)
            ff(11) = shp(2, 1) * shp(1, 2) * shp(2, 3)
            ff(12) = shp(3, 1) * shp(1, 2) * shp(2, 3)
 
            ff(13) = shp(1, 1) * shp(2, 2) * shp(2, 3)
            ff(14) = shp(2, 1) * shp(2, 2) * shp(2, 3)
            ff(15) = shp(3, 1) * shp(2, 2) * shp(2, 3)
 
            ff(16) = shp(1, 1) * shp(3, 2) * shp(2, 3)
            ff(17) = shp(2, 1) * shp(3, 2) * shp(2, 3)
            ff(18) = shp(3, 1) * shp(3, 2) * shp(2, 3)
 
            ff(19) = shp(1, 1) * shp(1, 2) * shp(3, 3)
            ff(20) = shp(2, 1) * shp(1, 2) * shp(3, 3)
            ff(21) = shp(3, 1) * shp(1, 2) * shp(3, 3)
 
            ff(22) = shp(1, 1) * shp(2, 2) * shp(3, 3)
            ff(23) = shp(2, 1) * shp(2, 2) * shp(3, 3)
            ff(24) = shp(3, 1) * shp(2, 2) * shp(3, 3)
 
            ff(25) = shp(1, 1) * shp(3, 2) * shp(3, 3)
            ff(26) = shp(2, 1) * shp(3, 2) * shp(3, 3)
            ff(27) = shp(3, 1) * shp(3, 2) * shp(3, 3)
 
            !step 2: find differentials of weights
            !-------------------------------------
 
            ! Linearize the FE shape functions wrt each psi(j) direction
            ! Note: only derivatives wrt psi(j) for any j are non-zero, rest are zero
 
            dshp(:, :, :) = 0.d0
            do j = 1, 3
                dshp(1, j, j) = half * (psi(j) - one) + half * psi(j)
                dshp(2, j, j) = -2.d0 * psi(j)
                dshp(3, j, j) = half * (psi(j) + one) + half * psi(j)
            end do
 
            ! Linearize 27 quadratic weights wrt each psi(j) dir, build from dshp
 
            loop_psi_j: do j = 1, 3
                dff(1, j) = dshp(1, 1, j) * shp(1, 2) * shp(1, 3) &
                            + shp(1, 1) * dshp(1, 2, j) * shp(1, 3) &
                            + shp(1, 1) * shp(1, 2) * dshp(1, 3, j)
 
                dff(2, j) = dshp(2, 1, j) * shp(1, 2) * shp(1, 3) &
                            + shp(2, 1) * dshp(1, 2, j) * shp(1, 3) &
                            + shp(2, 1) * shp(1, 2) * dshp(1, 3, j)
 
                dff(3, j) = dshp(3, 1, j) * shp(1, 2) * shp(1, 3) &
                            + shp(3, 1) * dshp(1, 2, j) * shp(1, 3) &
                            + shp(3, 1) * shp(1, 2) * dshp(1, 3, j)
 
                dff(4, j) = dshp(1, 1, j) * shp(2, 2) * shp(1, 3) &
                            + shp(1, 1) * dshp(2, 2, j) * shp(1, 3) &
                            + shp(1, 1) * shp(2, 2) * dshp(1, 3, j)
 
                dff(5, j) = dshp(2, 1, j) * shp(2, 2) * shp(1, 3) &
                            + shp(2, 1) * dshp(2, 2, j) * shp(1, 3) &
                            + shp(2, 1) * shp(2, 2) * dshp(1, 3, j)
 
                dff(6, j) = dshp(3, 1, j) * shp(2, 2) * shp(1, 3) &
                            + shp(3, 1) * dshp(2, 2, j) * shp(1, 3) &
                            + shp(3, 1) * shp(2, 2) * dshp(1, 3, j)
 
                dff(7, j) = dshp(1, 1, j) * shp(3, 2) * shp(1, 3) &
                            + shp(1, 1) * dshp(3, 2, j) * shp(1, 3) &
                            + shp(1, 1) * shp(3, 2) * dshp(1, 3, j)
 
                dff(8, j) = dshp(2, 1, j) * shp(3, 2) * shp(1, 3) &
                            + shp(2, 1) * dshp(3, 2, j) * shp(1, 3) &
                            + shp(2, 1) * shp(3, 2) * dshp(1, 3, j)
 
                dff(9, j) = dshp(3, 1, j) * shp(3, 2) * shp(1, 3) &
                            + shp(3, 1) * dshp(3, 2, j) * shp(1, 3) &
                            + shp(3, 1) * shp(3, 2) * dshp(1, 3, j)
 
                dff(10, j) = dshp(1, 1, j) * shp(1, 2) * shp(2, 3) &
                             + shp(1, 1) * dshp(1, 2, j) * shp(2, 3) &
                             + shp(1, 1) * shp(1, 2) * dshp(2, 3, j)
 
                dff(11, j) = dshp(2, 1, j) * shp(1, 2) * shp(2, 3) &
                             + shp(2, 1) * dshp(1, 2, j) * shp(2, 3) &
                             + shp(2, 1) * shp(1, 2) * dshp(2, 3, j)
 
                dff(12, j) = dshp(3, 1, j) * shp(1, 2) * shp(2, 3) &
                             + shp(3, 1) * dshp(1, 2, j) * shp(2, 3) &
                             + shp(3, 1) * shp(1, 2) * dshp(2, 3, j)
 
                dff(13, j) = dshp(1, 1, j) * shp(2, 2) * shp(2, 3) &
                             + shp(1, 1) * dshp(2, 2, j) * shp(2, 3) &
                             + shp(1, 1) * shp(2, 2) * dshp(2, 3, j)
 
                dff(14, j) = dshp(2, 1, j) * shp(2, 2) * shp(2, 3) &
                             + shp(2, 1) * dshp(2, 2, j) * shp(2, 3) &
                             + shp(2, 1) * shp(2, 2) * dshp(2, 3, j)
 
                dff(15, j) = dshp(3, 1, j) * shp(2, 2) * shp(2, 3) &
                             + shp(3, 1) * dshp(2, 2, j) * shp(2, 3) &
                             + shp(3, 1) * shp(2, 2) * dshp(2, 3, j)
 
                dff(16, j) = dshp(1, 1, j) * shp(3, 2) * shp(2, 3) &
                             + shp(1, 1) * dshp(3, 2, j) * shp(2, 3) &
                             + shp(1, 1) * shp(3, 2) * dshp(2, 3, j)
 
                dff(17, j) = dshp(2, 1, j) * shp(3, 2) * shp(2, 3) &
                             + shp(2, 1) * dshp(3, 2, j) * shp(2, 3) &
                             + shp(2, 1) * shp(3, 2) * dshp(2, 3, j)
 
                dff(18, j) = dshp(3, 1, j) * shp(3, 2) * shp(2, 3) &
                             + shp(3, 1) * dshp(3, 2, j) * shp(2, 3) &
                             + shp(3, 1) * shp(3, 2) * dshp(2, 3, j)
 
                dff(19, j) = dshp(1, 1, j) * shp(1, 2) * shp(3, 3) &
                             + shp(1, 1) * dshp(1, 2, j) * shp(3, 3) &
                             + shp(1, 1) * shp(1, 2) * dshp(3, 3, j)
 
                dff(20, j) = dshp(2, 1, j) * shp(1, 2) * shp(3, 3) &
                             + shp(2, 1) * dshp(1, 2, j) * shp(3, 3) &
                             + shp(2, 1) * shp(1, 2) * dshp(3, 3, j)
 
                dff(21, j) = dshp(3, 1, j) * shp(1, 2) * shp(3, 3) &
                             + shp(3, 1) * dshp(1, 2, j) * shp(3, 3) &
                             + shp(3, 1) * shp(1, 2) * dshp(3, 3, j)
 
                dff(22, j) = dshp(1, 1, j) * shp(2, 2) * shp(3, 3) &
                             + shp(1, 1) * dshp(2, 2, j) * shp(3, 3) &
                             + shp(1, 1) * shp(2, 2) * dshp(3, 3, j)
 
                dff(23, j) = dshp(2, 1, j) * shp(2, 2) * shp(3, 3) &
                             + shp(2, 1) * dshp(2, 2, j) * shp(3, 3) &
                             + shp(2, 1) * shp(2, 2) * dshp(3, 3, j)
 
                dff(24, j) = dshp(3, 1, j) * shp(2, 2) * shp(3, 3) &
                             + shp(3, 1) * dshp(2, 2, j) * shp(3, 3) &
                             + shp(3, 1) * shp(2, 2) * dshp(3, 3, j)
 
                dff(25, j) = dshp(1, 1, j) * shp(3, 2) * shp(3, 3) &
                             + shp(1, 1) * dshp(3, 2, j) * shp(3, 3) &
                             + shp(1, 1) * shp(3, 2) * dshp(3, 3, j)
 
                dff(26, j) = dshp(2, 1, j) * shp(3, 2) * shp(3, 3) &
                             + shp(2, 1) * dshp(3, 2, j) * shp(3, 3) &
                             + shp(2, 1) * shp(3, 2) * dshp(3, 3, j)
 
                dff(27, j) = dshp(3, 1, j) * shp(3, 2) * shp(3, 3) &
                             + shp(3, 1) * dshp(3, 2, j) * shp(3, 3) &
                             + shp(3, 1) * shp(3, 2) * dshp(3, 3, j)
 
            end do loop_psi_j
 
            ! Step 3: construct Jacobian d(R(psi(:))/d(psi(:)) and residue R(psi(:))
            ! ----------------------------------------------------------------------
 
            ! loop over x,y,z dirs, stored row-wise
            loop_xyz: do n = 1, 3
 
                ! construct LHS -d(R(psi(:))/d(psi(:))
 
                ! loop over each psi dir j
                do j = 1, 3
 
                    ! loop over nodes
                    l = 0
                    b(n, j) = 0.d0
 
                    do kkk = -1, 1
                        do jjj = -1, 1
                            do iii = -1, 1
                                l = l + 1
                                b(n, j) = b(n, j) + dff(l, j) * xelem(n, iii, jjj, kkk)
                            end do
                        end do
                    end do
                end do !j
 
                ! construct RHS (Xp - sum_i(ff_i * X_i), i =1,27) for nth row (x,y or z)
                l = 0
                b(n, 4) = xsearch(n)
 
                ! loop over nodes
                do kkk = -1, 1
                    do jjj = -1, 1
                        do iii = -1, 1
                            l = l + 1
                            b(n, 4) = b(n, 4) - ff(l) * xelem(n, iii, jjj, kkk)
                        end do
                    end do
                end do
 
            end do loop_xyz
 
            ! Get d(uvwQuadratic) weights
            ! invert 3x3 matrix returns solution in b(:,4)
            call matrixinv3by3(b, ok_flag)
            if (.not. ok_flag) stop 'Can not invert B in computeQuadraticWeights'
 
            ! update uvwQuadratic weights
            do n = 1, 3
                uvwquadratic(n) = uvwquadratic(n) + b(n, 4)
 
                ! sanity check: if weights lie outside [-1:1] reset to 0.5
                if ((uvwquadratic(n) - 1) * (uvwquadratic(n) + 1) > 0.d0) uvwquadratic(n) = half
            end do
 
            ! L2-norm of residue
            resid = 0.d0
            do n = 1, 3
                resid = resid + b(n, 4)**2
            end do
 
            niter = niter + 1
 
        end do newton_loop
    end subroutine computequadraticweights
    !
    !      Plain invert of A to solve Ax=B
    !      a(3,3) --- matrix to be inverted
    !      a(:,4) --- contains the right hand side vector, also stores
    !                 the final solution
 
    subroutine matrixinv3by3(a, ok_flag)
 
        use constants
        implicit none
 
        ! Input variables
        real(kind=realtype), intent(inout) :: a(3, 4)
        logical, intent(out) :: ok_flag
 
        ! Working variables
        integer(kind=intType) :: n
        real(kind=realtype) :: det, cofactor(3, 3), ainv(3, 3), rin(3, 1), rout(3, 1), myeps
 
        myeps = 1.0e-10
        ainv = 0.d0
 
        det = a(1, 1) * a(2, 2) * a(3, 3) &
              - a(1, 1) * a(2, 3) * a(3, 2) &
              - a(1, 2) * a(2, 1) * a(3, 3) &
              + a(1, 2) * a(2, 3) * a(3, 1) &
              + a(1, 3) * a(2, 1) * a(3, 2) &
              - a(1, 3) * a(2, 2) * a(3, 1)
 
        if (abs(det) <= myeps) then
            ainv = 0.0d0
            ok_flag = .false.
            return
        end if
 
        cofactor(1, 1) = +(a(2, 2) * a(3, 3) - a(2, 3) * a(3, 2))
        cofactor(1, 2) = -(a(2, 1) * a(3, 3) - a(2, 3) * a(3, 1))
        cofactor(1, 3) = +(a(2, 1) * a(3, 2) - a(2, 2) * a(3, 1))
        cofactor(2, 1) = -(a(1, 2) * a(3, 3) - a(1, 3) * a(3, 2))
        cofactor(2, 2) = +(a(1, 1) * a(3, 3) - a(1, 3) * a(3, 1))
        cofactor(2, 3) = -(a(1, 1) * a(3, 2) - a(1, 2) * a(3, 1))
        cofactor(3, 1) = +(a(1, 2) * a(2, 3) - a(1, 3) * a(2, 2))
        cofactor(3, 2) = -(a(1, 1) * a(2, 3) - a(1, 3) * a(2, 1))
        cofactor(3, 3) = +(a(1, 1) * a(2, 2) - a(1, 2) * a(2, 1))
 
        ainv = transpose(cofactor) / det
 
        ok_flag = .true.
 
        do n = 1, 3
            rin(n, 1) = a(n, 4)
        end do
 
        !save solution on a(:,4)
        rout = matmul(ainv, rin)
 
        do n = 1, 3
            a(n, 4) = rout(n, 1)
        end do
 
    end subroutine matrixinv3by3
 
    subroutine fringereduction(level, sps)
 
        use constants
        use blockpointers, only: ndom, il, jl, kl, status, fringeptr
        use stencils, only: visc_drdw_stencil, n_visc_drdw
        use utils, only: setpointers
        implicit none
 
        ! Input/Output
        integer(kind=intType), intent(in) :: level, sps
 
        ! Working
        integer(kind=intType) :: i, j, k, nn, ii, jj, kk, i_stencil
        logical :: computeCellFound
 
        do nn = 1, ndom
            call setpointers(nn, level, sps)
 
            do k = 2, kl
                do j = 2, jl
                    do i = 2, il
 
                        ! Check if this cell is a fringe:
                        if (isreceiver(status(i, j, k))) then
 
                            computecellfound = .false.
 
                            stencilloop2: do i_stencil = 1, n_visc_drdw
                                ii = visc_drdw_stencil(i_stencil, 1) + i
                                jj = visc_drdw_stencil(i_stencil, 2) + j
                                kk = visc_drdw_stencil(i_stencil, 3) + k
 
                                if (iscompute(status(ii, jj, kk))) then
                                    ! This is a compute cell
                                    computecellfound = .true.
                                end if
                            end do stencilloop2
 
                            if (.not. computecellfound) then
                                ! This cell is a hole no compute cell
                                ! surrounding a fringe, we can hard iblank it.
                                call setishole(status(i, j, k), .true.)
                                call setiscompute(status(i, j, k), .false.)
                                call setisreceiver(status(i, j, k), .false.)
                                fringeptr(1, i, j, k) = 0
 
                            end if
                        end if
                    end do
                end do
            end do
        end do
 
    end subroutine fringereduction
 
    subroutine irregularcellcorrection(level, sps)
 
        use constants
        use blockpointers, only: ndom, il, jl, kl, status, ie, je, ke, forcedrecv
        use utils, only: setpointers
        implicit none
 
        ! Input/Output
        integer(kind=intType), intent(in) :: level, sps
 
        ! Working
        integer(kind=intType) :: i, j, k, nn
 
        do nn = 1, ndom
            call setpointers(nn, level, sps)
 
            do k = 2, kl
                do j = 2, jl
                    do i = 2, il
                        if (isdonor(status(i, j, k)) .and. &
                            isreceiver(status(i, j, k))) then
 
                            ! Clear the fringe
                            call setisdonor(status(i, j, k), .false.)
                            call setisreceiver(status(i, j, k), .false.)
                            call setiscompute(status(i, j, k), .true.)
                        end if
                    end do
                end do
            end do
        end do
 
    end subroutine irregularcellcorrection
 
    function checkoversetpresent()
 
        ! This routine determines if there are any overset boundaries
        ! present in the mesh.
 
        use constants
        use blockpointers, only: ndom, nbocos, bctype
        use communication, only: adflow_comm_world
        use utils, only: setpointers, echk
        implicit none
 
        ! Function
        logical :: checkoversetpresent, local
 
        ! Working
        integer(Kind=intType) :: nn, mm, ierr
 
        local = .false.
        do nn = 1, ndom
            call setpointers(nn, 1_inttype, 1_inttype)
 
            do mm = 1, nbocos
                if (bctype(mm) == oversetouterbound) then
                    local = .true.
                end if
            end do
        end do
 
        call mpi_allreduce(local, checkoversetpresent, 1, mpi_logical, mpi_lor, adflow_comm_world, ierr)
        call echk(ierr, __file__, __line__)
 
    end function checkoversetpresent
 
    subroutine setiblankarray(level, sps)
 
        use constants
        use block
        use blockpointers, only: ndom, il, jl, kl, status, iblank, flowdoms, ie, je, ke, forcedrecv
        use communication, only: myid, commpatterncell_2nd, internalcell_2nd, &
                                 adflow_comm_world
        use utils, only: setpointers, echk
        use haloexchange, only: whalo1to1intgeneric
        implicit none
 
        ! Input/Output
        integer(kind=intType), intent(in) :: level, sps
 
        ! Working
        integer(kind=intType) :: i, j, k, nn
        integer(kind=intType) :: nCompute, nFringe, nBlank, nFloodSeed, nFlooded, nExplicitBlanked
        integer(kind=intType) :: counts(6), ierr
        type(fringetype) :: fringe
        ncompute = 0
        nfringe = 0
        nblank = 0
        nfloodseed = 0
        nflooded = 0
        nexplicitblanked = 0
 
        do nn = 1, ndom
            call setpointers(nn, level, sps)
 
            do k = 2, kl
                do j = 2, jl
                    do i = 2, il
 
                        if (iblank(i, j, k) == -4) then
                            nexplicitblanked = nexplicitblanked + 1
 
                        else if (isreceiver(status(i, j, k))) then
                            iblank(i, j, k) = -1
                            nfringe = nfringe + 1
 
                        else if (isfloodseed(status(i, j, k))) then
                            iblank(i, j, k) = -3
                            nfloodseed = nfloodseed + 1
 
                        else if (isflooded(status(i, j, k))) then
                            iblank(i, j, k) = -2
                            nflooded = nflooded + 1
 
                        else if (ishole(status(i, j, k))) then
                            iblank(i, j, k) = 0
                            nblank = nblank + 1
 
                        else
                            ! We need to explictly make sure forced receivers
                            ! *NEVER EVER EVER EVER* get set as compute cells.
                            if (forcedrecv(i, j, k) .ne. 0) then
                                ! This is REALLY Bad. This cell must be a forced
                                ! receiver but never found a donor.
                                iblank(i, j, k) = 0
                                nblank = nblank + 1
                            else
                                ! Compute cell
                                ncompute = ncompute + 1
                                iblank(i, j, k) = 1
                            end if
                        end if
                    end do
                end do
            end do
        end do
 
        ! Update the iblank info.
        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, commpatterncell_2nd, internalcell_2nd)
 
        call mpi_reduce((/ncompute, nfringe, nblank, nflooded, nfloodseed, nexplicitblanked/), &
                        counts, 6, adflow_integer, mpi_sum, 0, adflow_comm_world, ierr)
        call echk(ierr, __file__, __line__)
 
        if (myid == 0) then
            print *, '+--------------------------------+'
            print *, '| Compute Cells           :', counts(1)
            print *, '| Fringe Cells            :', counts(2)
            print *, '| Blanked Cells           :', counts(3)
            print *, '| Explicitly Blanked Cells:', counts(6)
            print *, '| Flooded   Cells         :', counts(4)
            print *, '| FloodSeed Cells         :', counts(5)
            print *, '+--------------------------------+'
        end if
    end subroutine setiblankarray
 
    subroutine dumpiblank(level, sps)
 
        use constants
        use blockpointers, only: il, jl, kl, x, ndom, iblank
        use communication, only: myid
        use utils, only: setpointers
        implicit none
 
        ! Input/Output
        integer(kind=intType), intent(in) :: level, sps
 
        ! Working
        integer(kind=intType) :: i, j, k, nn
        real(kind=realtype) :: xp(3)
        character(80) :: fileName
 
        write (filename, "(a,I2.2,a)") "proc_", myid, ".dat"
        open (unit=19, file=trim(filename), form='formatted')
 
        do nn = 1, ndom
            call setpointers(nn, level, sps)
            do k = 2, kl
                do j = 2, jl
                    do i = 2, il
 
                        ! Compute the cell center:
                        xp = 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, :))
 
                        write (19, "(ES18.10, ES18.10, ES18.10, I3)") xp(1), xp(2), xp(3), iblank(i, j, k)
                    end do
                end do
            end do
        end do
 
        close (19)
 
    end subroutine dumpiblank
 
    subroutine getworkarray(overlap, work)
 
        use constants
        use communication, only: myid
        use oversetdata, only: csrmatrix, ndomtotal
        implicit none
 
        ! Input/Output
        type(csrmatrix), intent(in) :: overlap
        integer(kind=intType), dimension(:, :), allocatable :: work
 
        ! Local variables
        integer(kind=intType) :: nWork, jj, iDom, jDom
 
        nwork = 0
        do jj = 1, overlap%nnz
            if (overlap%assignedProc(jj) == myid) then
                nwork = nwork + 1
            end if
        end do
        allocate (work(4, nwork))
 
        nwork = 0
        do idom = 1, ndomtotal
            do jj = overlap%rowPtr(idom), overlap%rowPtr(idom + 1) - 1
                jdom = overlap%colInd(jj)
                if (overlap%assignedProc(jj) == myid) then
                    nwork = nwork + 1
                    work(1, nwork) = idom
                    work(2, nwork) = jdom
                    work(3, nwork) = jj
                    work(4, nwork) = 0
                end if
            end do
        end do
    end subroutine getworkarray
 
    subroutine writeoversetwall(oWall, fName)
 
        ! debug routine to dumb an owall to a file
        use constants
        use oversetdata
        implicit none
        type(oversetwall) :: oWall
        character(len=*), intent(in) :: fName
 
        integer(kind=intType) :: i, iDim
        open (unit=19, file=trim(fname), form='formatted')
        write (19, *) 'Variables = "CoordinateX" "CoordianteY" "CoordinateZ"'
        write (19, *) "Zone"
        write (19, *) "Nodes = ", owall%nNodes, " Elements= ", owall%nCells, " ZONETYPE=FEQUADRILATERAL"
        write (19, *) "DATAPACKING=BLOCK"
 
        do idim = 1, 3
            do i = 1, owall%nNodes
                write (19, *) owall%x(idim, i)
            end do
        end do
 
        do i = 1, owall%nCells
            write (19, *) owall%conn(1, i), owall%conn(2, i), owall%conn(3, i), owall%conn(4, i)
        end do
        close (19)
    end subroutine writeoversetwall
 
    subroutine getoversetiblank(blkList, n)
 
        ! This routine gathers a list of the iblank status of every
        ! compute cell in the original CGNS ordering. It then returns the
        ! full list on the root processor. Therefore, this routine is not
        ! memory or computation scalable. It is only meant to be used a
        ! debugging tool to ensure that the overset connectivity as
        ! computed with given parallel decomposition is the same as a
        ! different parallel decomposition.
 
        use constants
        use blockpointers, only: ndom, il, jl, kl, ibegor, jbegor, kbegor, &
                                 nbkglobal, iblank
        use cgnsgrid, only: cgnsdoms, cgnsndom
        use communication, only: adflow_comm_world, myid
        use inputtimespectral, only: ntimeintervalsspectral
        use utils, only: setpointers, echk, terminate
#include <petsc/finclude/petsc.h>
        use petsc
        implicit none
 
        ! Input/Output
        integer(kind=intType), intent(in) :: n
        integer(kind=intType), dimension(n), intent(out) :: blkList
 
        ! Working
 
        ! Working parameters
        integer(kind=intType) :: i, j, k, ierr, l, nx_cg, ny_cg, nz_cg
        integer(kind=intType) :: ii, indx, indy, indz, nn, cgnsInd
        integer(kind=intType), allocatable, dimension(:) :: cellOffset
        real(kind=realtype), dimension(:), pointer :: localptr
        real(kind=realtype) :: vals(5)
        vec cgnsvec
 
        ! This routine cannot be used in timespectral mode
        if (ntimeintervalsspectral > 1) then
            call terminate('getOversetIBlank', 'This routine can only be used '&
                 &'with 1 spectral instance')
        end if
 
        allocate (celloffset(cgnsndom + 1))
        celloffset(1) = 0
        do nn = 1, cgnsndom
            celloffset(nn + 1) = celloffset(nn) + &
                                 cgnsdoms(nn)%nx * cgnsdoms(nn)%ny * cgnsdoms(nn)%nz
        end do
 
        ! Create the CGNSVector
        if (myid == 0) then
            call veccreatempi(adflow_comm_world, celloffset(cgnsndom + 1) * 5, &
                              petsc_determine, cgnsvec, ierr)
            call echk(ierr, __file__, __line__)
 
        else
            call veccreatempi(adflow_comm_world, 0, petsc_determine, cgnsvec, ierr)
            call echk(ierr, __file__, __line__)
        end if
 
        call vecsetblocksize(cgnsvec, 5, ierr)
        call echk(ierr, __file__, __line__)
        ii = 0
        do nn = 1, ndom
            call setpointers(nn, 1, 1)
            do k = 2, kl
                do j = 2, jl
                    do i = 2, il
 
                        nx_cg = cgnsdoms(nbkglobal)%nx
                        ny_cg = cgnsdoms(nbkglobal)%ny
                        nz_cg = cgnsdoms(nbkglobal)%nz
 
                        indx = ibegor + i - 2
                        indy = jbegor + j - 2
                        indz = kbegor + k - 2
 
                        ! cgnsInd is zero-based
                        cgnsind = celloffset(nbkglobal) + &
                                  (indz - 1) * ny_cg * nx_cg + &
                                  (indy - 1) * nx_cg + &
                                  (indx - 1)
                        vals(1) = real(iblank(i, j, k))
                        if (vals(1) <= -2) then
                            vals(1) = 0
                        end if
                        vals(2) = real(nbkglobal)
                        vals(3) = real(indx)
                        vals(4) = real(indy)
                        vals(5) = real(indz)
                        call vecsetvaluesblocked(cgnsvec, 1, (/cgnsind/), vals, insert_values, ierr)
                        call echk(ierr, __file__, __line__)
                    end do
                end do
            end do
        end do
 
        call vecassemblybegin(cgnsvec, ierr)
        call echk(ierr, __file__, __line__)
 
        call vecassemblyend(cgnsvec, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Get the local vector pointer. Only the root proc actually has
        ! values.
        call vecgetarrayf90(cgnsvec, localptr, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Convert back to integer.
        do i = 1, size(localptr)
            blklist(i) = int(localptr(i))
        end do
 
        call vecrestorearrayf90(cgnsvec, localptr, ierr)
        call echk(ierr, __file__, __line__)
 
    end subroutine getoversetiblank
 
#endif
 
    ! --------------------------------------------------
    !           Tapenade Routine BELOW this point
    ! --------------------------------------------------
 
    subroutine fractoweights(frac, weights)
        use constants
        implicit none
        real(kind=realtype), intent(in), dimension(3) :: frac
        real(kind=realtype), intent(out), dimension(8) :: weights
 
        weights(1) = (one - frac(1)) * (one - frac(2)) * (one - frac(3))
        weights(2) = (frac(1)) * (one - frac(2)) * (one - frac(3))
        weights(3) = (one - frac(1)) * (frac(2)) * (one - frac(3))
        weights(4) = (frac(1)) * (frac(2)) * (one - frac(3))
        weights(5) = (one - frac(1)) * (one - frac(2)) * (frac(3))
        weights(6) = (frac(1)) * (one - frac(2)) * (frac(3))
        weights(7) = (one - frac(1)) * (frac(2)) * (frac(3))
        weights(8) = (frac(1)) * (frac(2)) * (frac(3))
    end subroutine fractoweights
 
    subroutine fractoweights2(frac, weights)
        use constants
        implicit none
        real(kind=realtype), intent(in), dimension(3) :: frac
        real(kind=realtype), intent(out), dimension(8) :: weights
 
        weights(1) = (one - frac(1)) * (one - frac(2)) * (one - frac(3))
        weights(2) = (frac(1)) * (one - frac(2)) * (one - frac(3))
        weights(3) = (frac(1)) * (frac(2)) * (one - frac(3))
        weights(4) = (one - frac(1)) * (frac(2)) * (one - frac(3))
 
        weights(5) = (one - frac(1)) * (one - frac(2)) * (frac(3))
        weights(6) = (frac(1)) * (one - frac(2)) * (frac(3))
        weights(7) = (frac(1)) * (frac(2)) * (frac(3))
        weights(8) = (one - frac(1)) * (frac(2)) * (frac(3))
 
    end subroutine fractoweights2
 
    subroutine newtonupdate(xCen, blk, frac0, frac)
 
        ! This routine performs the newton update to recompute the new
        ! "frac" (u,v,w) for the point xCen. The actual search is performed
        ! on the the dual cell formed by the cell centers of the 3x3x3 block
        ! of primal nodes. This routine is AD'd with tapenade in both
        ! forward and reverse.
 
        use constants
        implicit none
 
        ! Input
        real(kind=realtype), dimension(3), intent(in) :: xcen
        real(kind=realtype), dimension(3, 3, 3, 3), intent(in) :: blk
        real(kind=realtype), dimension(3), intent(in) :: frac0
        ! Output
        real(kind=realtype), dimension(3), intent(out) :: frac
 
        ! Working
        real(kind=realtype), dimension(3, 1:8) :: xn
        real(kind=realtype) :: u, v, w, uv, uw, vw, wvu, du, dv, dw
        real(kind=realtype) :: a11, a12, a13, a21, a22, a23, a31, a32, a33, val
        real(kind=realtype) :: f(3), x(3)
        integer(kind=intType), dimension(8), parameter :: indices = [1, 2, 4, 3, 5, 6, 8, 7]
        integer(kind=intType) :: i, j, k, ii, ll
        real(kind=realtype), parameter :: adteps = 1.e-25_realtype
        real(kind=realtype), parameter :: thresconv = 1.e-10_realtype
 
        ! Compute the cell center locations for the 8 nodes describing the
        ! dual cell. Note that this must be counter-clockwise ordering.
 
        ii = 0
        do k = 1, 2
            do j = 1, 2
                do i = 1, 2
                    ii = ii + 1
                    xn(:, indices(ii)) = eighth * ( &
                                         blk(i, j, k, :) + &
                                         blk(i + 1, j, k, :) + &
                                         blk(i, j + 1, k, :) + &
                                         blk(i + 1, j + 1, k, :) + &
                                         blk(i, j, k + 1, :) + &
                                         blk(i + 1, j, k + 1, :) + &
                                         blk(i, j + 1, k + 1, :) + &
                                         blk(i + 1, j + 1, k + 1, :))
                end do
            end do
        end do
 
        ! Compute the coordinates relative to node 1.
 
        do i = 2, 8
            xn(:, i) = xn(:, i) - xn(:, 1)
        end do
 
        ! Compute the location of our seach point relative to the first node.
        x = xcen - xn(:, 1)
 
        ! Modify the coordinates of node 3, 6, 8 and 7 such that
        ! they correspond to the weights of the u*v, u*w, v*w and
        ! u*v*w term in the transformation respectively.
 
        xn(1, 7) = xn(1, 7) + xn(1, 2) + xn(1, 4) + xn(1, 5) &
                   - xn(1, 3) - xn(1, 6) - xn(1, 8)
        xn(2, 7) = xn(2, 7) + xn(2, 2) + xn(2, 4) + xn(2, 5) &
                   - xn(2, 3) - xn(2, 6) - xn(2, 8)
        xn(3, 7) = xn(3, 7) + xn(3, 2) + xn(3, 4) + xn(3, 5) &
                   - xn(3, 3) - xn(3, 6) - xn(3, 8)
 
        xn(1, 3) = xn(1, 3) - xn(1, 2) - xn(1, 4)
        xn(2, 3) = xn(2, 3) - xn(2, 2) - xn(2, 4)
        xn(3, 3) = xn(3, 3) - xn(3, 2) - xn(3, 4)
 
        xn(1, 6) = xn(1, 6) - xn(1, 2) - xn(1, 5)
        xn(2, 6) = xn(2, 6) - xn(2, 2) - xn(2, 5)
        xn(3, 6) = xn(3, 6) - xn(3, 2) - xn(3, 5)
 
        xn(1, 8) = xn(1, 8) - xn(1, 4) - xn(1, 5)
        xn(2, 8) = xn(2, 8) - xn(2, 4) - xn(2, 5)
        xn(3, 8) = xn(3, 8) - xn(3, 4) - xn(3, 5)
 
        ! Set the starting values of u, v and w based on our previous values
 
        u = frac0(1); v = frac0(2); w = frac0(3);
        ! The Newton algorithm to determine the parametric
        ! weights u, v and w for the given coordinate.
 
        newtonhexa: do ll = 1, 15
 
            ! Compute the RHS.
 
            uv = u * v; uw = u * w; vw = v * w; wvu = u * v * w
 
            f(1) = xn(1, 2) * u + xn(1, 4) * v + xn(1, 5) * w &
                   + xn(1, 3) * uv + xn(1, 6) * uw + xn(1, 8) * vw &
                   + xn(1, 7) * wvu - x(1)
            f(2) = xn(2, 2) * u + xn(2, 4) * v + xn(2, 5) * w &
                   + xn(2, 3) * uv + xn(2, 6) * uw + xn(2, 8) * vw &
                   + xn(2, 7) * wvu - x(2)
            f(3) = xn(3, 2) * u + xn(3, 4) * v + xn(3, 5) * w &
                   + xn(3, 3) * uv + xn(3, 6) * uw + xn(3, 8) * vw &
                   + xn(3, 7) * wvu - x(3)
 
            ! Compute the Jacobian.
 
            a11 = xn(1, 2) + xn(1, 3) * v + xn(1, 6) * w + xn(1, 7) * vw
            a12 = xn(1, 4) + xn(1, 3) * u + xn(1, 8) * w + xn(1, 7) * uw
            a13 = xn(1, 5) + xn(1, 6) * u + xn(1, 8) * v + xn(1, 7) * uv
 
            a21 = xn(2, 2) + xn(2, 3) * v + xn(2, 6) * w + xn(2, 7) * vw
            a22 = xn(2, 4) + xn(2, 3) * u + xn(2, 8) * w + xn(2, 7) * uw
            a23 = xn(2, 5) + xn(2, 6) * u + xn(2, 8) * v + xn(2, 7) * uv
 
            a31 = xn(3, 2) + xn(3, 3) * v + xn(3, 6) * w + xn(3, 7) * vw
            a32 = xn(3, 4) + xn(3, 3) * u + xn(3, 8) * w + xn(3, 7) * uw
            a33 = xn(3, 5) + xn(3, 6) * u + xn(3, 8) * v + xn(3, 7) * uv
 
            ! Compute the determinant. Make sure that it is not zero
            ! and invert the value. The cut off is needed to be able
            ! to handle exceptional cases for degenerate elements.
 
            val = a11 * (a22 * a33 - a32 * a23) + a21 * (a13 * a32 - a12 * a33) &
                  + a31 * (a12 * a23 - a13 * a22)
            val = sign(one, val) / max(abs(val), adteps)
 
            ! Compute the new values of u, v and w.
 
            du = val * ((a22 * a33 - a23 * a32) * f(1) &
                        + (a13 * a32 - a12 * a33) * f(2) &
                        + (a12 * a23 - a13 * a22) * f(3))
            dv = val * ((a23 * a31 - a21 * a33) * f(1) &
                        + (a11 * a33 - a13 * a31) * f(2) &
                        + (a13 * a21 - a11 * a23) * f(3))
            dw = val * ((a21 * a32 - a22 * a31) * f(1) &
                        + (a12 * a31 - a11 * a32) * f(2) &
                        + (a11 * a22 - a12 * a21) * f(3))
 
            u = u - du; v = v - dv; w = w - dw
 
            ! Exit the loop if the update of the parametric
            ! weights is below the threshold
 
            val = sqrt(du * du + dv * dv + dw * dw)
            if (val <= thresconv) then
                exit newtonhexa
            end if
 
        end do newtonhexa
 
        ! We would *like* that all solutions fall inside the hexa, but we
        ! can't be picky here since we are not changing the donors. So
        ! whatever the u,v,w is we have to accept. Even if it is greater than
        ! 1 or less than zero, it shouldn't be by much.
 
        frac(1) = u
        frac(2) = v
        frac(3) = w
 
    end subroutine newtonupdate
 
end module oversetutilities