gridChecking.F90

Go to the documentation of this file.

module gridchecking
contains
    subroutine checkfaces
        !
        !       checkFaces determines whether or not a boundary condition or
        !       a connectivity has been specified for all block faces. If this
        !       is not the case, the corresponding blocks are printed and the
        !       code terminates.
        !
        use constants
        use blockpointers, only: ndom, nbkglobal
        use cgnsgrid, only: cgnsndom, cgnsdoms
        use communication, only: adflow_comm_world, myid, nproc
        use inputtimespectral, only: ntimeintervalsspectral
        use utils, only: setpointers, terminate
        use partitionmod, only: sortbadentities
        use commonformats, only: strings
        implicit none
        !
        !      Local variables.
        !
        integer :: ierr, size
 
        integer, dimension(nProc) :: recvcounts, displs
 
        integer(kind=intType) :: mm, nn, sps, multiple
        integer(kind=intType) :: nBad, faceID, nBadGlobal
 
        integer(kind=intType), dimension(nProc) :: counts
        integer(kind=intType), &
            dimension(4, nDom*nTimeIntervalsSpectral) :: bad
 
        integer(kind=intType), dimension(:, :), allocatable :: badGlobal
 
        real(kind=realtype) :: dummy(1)
 
        logical :: blockIsBad
 
        character(len=7) :: intString
 
        !       Determine the local bad blocks.
        !
        ! Loop over the number of spectral solutions to be checked and
        ! the number of local blocks.
 
        nbad = 0
        do sps = 1, ntimeintervalsspectral
            do nn = 1, ndom
 
                ! Check if the current block is okay.
 
                call setpointers(nn, 1_inttype, sps)
                call checkfacesblock(blockisbad, faceid, multiple)
 
                ! If the block is bad, update nBad and store the info in bad.
 
                if (blockisbad) then
                    nbad = nbad + 1
                    bad(1, nbad) = nbkglobal
                    bad(2, nbad) = faceid
                    bad(3, nbad) = multiple
                    bad(4, nbad) = sps
                end if
 
            end do
        end do
        !
        !       Determine the global number of bad blocks and gather this
        !       information.
        !
        ! Determine the number of bad blocks per processor.
 
        call mpi_allgather(nbad, 1, adflow_integer, counts, 1, &
                           adflow_integer, adflow_comm_world, ierr)
 
        ! Determine the global number of bad blocks and the arrays
        ! recvcounts and displs needed for the call to allgatherv.
 
        nbadglobal = counts(1)
        recvcounts(1) = 4 * counts(1)
        displs(1) = 0
 
        do nn = 2, nproc
            nbadglobal = nbadglobal + counts(nn)
            recvcounts(nn) = 4 * counts(nn)
            displs(nn) = displs(nn - 1) + recvcounts(nn - 1)
        end do
 
        ! Allocate the memory to store all the bad blocks.
 
        allocate (badglobal(4, nbadglobal), stat=ierr)
        if (ierr /= 0) &
            call terminate("checkFaces", &
                           "Memory allocation failure for badGlobal")
 
        ! Gather the data.
 
        size = 4 * nbad
        call mpi_allgatherv(bad, size, adflow_integer, badglobal, &
                            recvcounts, displs, adflow_integer, &
                            adflow_comm_world, ierr)
 
        ! Sort the bad blocks and get rid of the multiple entries.
        ! The last argument is .false. to indicate that only the
        ! integers must be sorted; dummy is passed for consistency.
 
        call sortbadentities(nbadglobal, badglobal, dummy, .false.)
 
        ! If bad blocks are present, print them to stdout and terminate.
 
        testbadpresent: if (nbadglobal > 0) then
 
            ! The data is only written by processor 0.
 
            testrootproc: if (myid == 0) then
 
                ! Write a header.
 
                write (intstring, "(i6)") nbadglobal
                intstring = adjustl(intstring)
 
                print "(a)", "#"
                print strings, "# Found ", trim(intstring), " blocks for which the boundary or connectivity information"
                print "(a)", "# is not correct. Here is the list of bad blocks"
                print "(a)", "#"
 
                ! Write the bad blocks.
 
                do mm = 1, nbadglobal
 
                    ! Abbreviate the contents of badGlobal a bit easier.
 
                    nn = badglobal(1, mm)
                    faceid = badglobal(2, mm)
                    multiple = badglobal(3, mm)
                    sps = badglobal(4, mm)
 
                    ! If multiple spectral solutions are read, write the info
                    ! about it. Otherwise just start the line with # sign.
 
                    if (ntimeintervalsspectral > 1) then
                        write (intstring, "(i6)") nbadglobal
                        intstring = adjustl(intstring)
                        write (*, strings, advance="no") "# Spectral grid ", trim(intstring), ", "
                    else
                        write (*, "(a)", advance="no") "# "
                    end if
 
                    ! Write the error message, depending on the value of
                    ! faceID and multiple.
 
                    if (multiple == 0) then
 
                        select case (faceid)
                        case (imin)
                            print strings, "Zone ", trim(cgnsdoms(nn)%zoneName), ": iMin block face not fully described"
                        case (imax)
                            print strings, "Zone ", trim(cgnsdoms(nn)%zoneName), ": iMax block face not fully described"
                        case (jmin)
                            print strings, "Zone ", trim(cgnsdoms(nn)%zoneName), ": jMin block face not fully described"
                        case (jmax)
                            print strings, "Zone ", trim(cgnsdoms(nn)%zoneName), ": jMax block face not fully described"
                        case (kmin)
                            print strings, "Zone ", trim(cgnsdoms(nn)%zoneName), ": kMin block face not fully described"
                        case (kmax)
                            print strings, "Zone ", trim(cgnsdoms(nn)%zoneName), ": kMax block face not fully described"
                        case default
                            print strings, "Zone ", trim(cgnsdoms(nn)%zoneName), &
                                ": Multiple block faces not fully described"
                        end select
 
                    else
 
                        select case (faceid)
                        case (imin)
                            print strings, "Zone ", trim(cgnsdoms(nn)%zoneName), ": iMin block face multiple described"
                        case (imax)
                            print strings, "Zone ", trim(cgnsdoms(nn)%zoneName), ": iMax block face multiple described"
                        case (jmin)
                            print strings, "Zone ", trim(cgnsdoms(nn)%zoneName), ": jMin block face multiple described"
                        case (jmax)
                            print strings, "Zone ", trim(cgnsdoms(nn)%zoneName), ": jMax block face multiple described"
                        case (kmin)
                            print strings, "Zone ", trim(cgnsdoms(nn)%zoneName), ": kMin block face multiple described"
                        case (kmax)
                            print strings, "Zone ", trim(cgnsdoms(nn)%zoneName), ": kMax block face multiple described"
                        case default
                            print strings, "Zone ", trim(cgnsdoms(nn)%zoneName), &
                                ": Multiple block faces multiple described"
                        end select
 
                    end if
 
                end do
 
                ! Terminate.
 
                call terminate("checkFaces", &
                               "Wrong block boundary info found")
 
            end if testrootproc
 
            ! The other processors wait to get killed.
 
            call mpi_barrier(adflow_comm_world, ierr)
 
        end if testbadpresent
 
        ! Deallocate the memory of badGlobal.
 
        deallocate (badglobal, stat=ierr)
        if (ierr /= 0) &
            call terminate("checkFaces", &
                           "Deallocation failure for badGlobal")
 
    end subroutine checkfaces
 
    !      ==================================================================
 
    subroutine checkfacesblock(blockIsBad, faceID, multiple)
        !
        !       checkFacesBlock checks if for the currently active block all
        !       the necessary boundary and connectivity info is specified.
        !       If not, blockIsBad is set to .true. and the block face ID
        !       which is bad, is returned as well.
        !
        use constants
        use blockpointers
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(out) :: faceID, multiple
        logical, intent(out) :: blockIsBad
        !
        !      Local variables.
        !
        integer, dimension(2:jl, 2:kl), target :: iMinFace, iMaxFace
        integer, dimension(2:il, 2:kl), target :: jMinFace, jMaxFace
        integer, dimension(2:il, 2:jl), target :: kMinFace, kMaxFace
 
        integer, dimension(:, :), pointer :: face
 
        integer(kind=intType) :: nn, i, j, k, nBad
        integer(kind=intType) :: istart, iEnd, jStart, jEnd
 
        logical :: iMinBad, iMaxBad, jMinBad, jMaxBad
        logical :: kMinBad, kMaxBad
 
        ! Initialize iMinFace, etc to 0. These arrays will contain the
        ! number of specified connectivities and BC's for each face on
        ! the block boundary.
 
        iminface = 0; imaxface = 0
        jminface = 0; jmaxface = 0
        kminface = 0; kmaxface = 0
 
        ! Loop over the number of subfaces of this block.
 
        do nn = 1, nsubface
 
            ! Determine the block face on which the subface is located and
            ! set a couple of variables accordingly. Note that the nodal
            ! ranges are used to determine the correct range of faces,
            ! because the actual cell range can contain halo cells.
 
            select case (bcfaceid(nn))
            case (imin)
                face => iminface; istart = min(jnbeg(nn), jnend(nn)) + 1
                iend = max(jnbeg(nn), jnend(nn))
                jstart = min(knbeg(nn), knend(nn)) + 1
                jend = max(knbeg(nn), knend(nn))
            case (imax)
                face => imaxface; istart = min(jnbeg(nn), jnend(nn)) + 1
                iend = max(jnbeg(nn), jnend(nn))
                jstart = min(knbeg(nn), knend(nn)) + 1
                jend = max(knbeg(nn), knend(nn))
            case (jmin)
                face => jminface; istart = min(inbeg(nn), inend(nn)) + 1
                iend = max(inbeg(nn), inend(nn))
                jstart = min(knbeg(nn), knend(nn)) + 1
                jend = max(knbeg(nn), knend(nn))
            case (jmax)
                face => jmaxface; istart = min(inbeg(nn), inend(nn)) + 1
                iend = max(inbeg(nn), inend(nn))
                jstart = min(knbeg(nn), knend(nn)) + 1
                jend = max(knbeg(nn), knend(nn))
            case (kmin)
                face => kminface; istart = min(inbeg(nn), inend(nn)) + 1
                iend = max(inbeg(nn), inend(nn))
                jstart = min(jnbeg(nn), jnend(nn)) + 1
                jend = max(jnbeg(nn), jnend(nn))
            case (kmax)
                face => kmaxface; istart = min(inbeg(nn), inend(nn)) + 1
                iend = max(inbeg(nn), inend(nn))
                jstart = min(jnbeg(nn), jnend(nn)) + 1
                jend = max(jnbeg(nn), jnend(nn))
            end select
 
            ! Loop over the faces and update their counter.
 
            do j = jstart, jend
                do i = istart, iend
                    face(i, j) = face(i, j) + 1
                end do
            end do
 
        end do
 
        ! Determine the bad block faces.
 
        multiple = 0
        iminbad = .false.
        imaxbad = .false.
 
        ! iMin and iMax face.
 
        do k = 2, kl
            do j = 2, jl
                if (iminface(j, k) /= 1) then
                    multiple = iminface(j, k)
                    faceid = imin
                    iminbad = .true.
                end if
 
                if (imaxface(j, k) /= 1) then
                    multiple = imaxface(j, k)
                    faceid = imax
                    imaxbad = .true.
                end if
            end do
        end do
 
        ! jMin and jMax face.
 
        jminbad = .false.
        jmaxbad = .false.
 
        do k = 2, kl
            do i = 2, il
                if (jminface(i, k) /= 1) then
                    multiple = jminface(i, k)
                    faceid = jmin
                    jminbad = .true.
                end if
 
                if (jmaxface(i, k) /= 1) then
                    multiple = jmaxface(i, k)
                    faceid = jmax
                    jmaxbad = .true.
                end if
            end do
        end do
 
        ! kMin and kMax face.
 
        kminbad = .false.
        kmaxbad = .false.
 
        do j = 2, jl
            do i = 2, il
                if (kminface(i, j) /= 1) then
                    multiple = kminface(i, j)
                    faceid = kmin
                    kminbad = .true.
                end if
 
                if (kmaxface(i, j) /= 1) then
                    multiple = kmaxface(i, j)
                    faceid = kmax
                    kmaxbad = .true.
                end if
            end do
        end do
 
        ! Determine the number of bad block faces.
 
        nbad = 0
        if (iminbad) nbad = nbad + 1
        if (imaxbad) nbad = nbad + 1
        if (jminbad) nbad = nbad + 1
        if (jmaxbad) nbad = nbad + 1
        if (kminbad) nbad = nbad + 1
        if (kmaxbad) nbad = nbad + 1
 
        ! Set blockIsBad if bad faces are present and correct the face
        ! id to something weird if multiple bad faces are present.
 
        blockisbad = .false.
        if (nbad > 0) blockisbad = .true.
        if (nbad > 1) faceid = huge(faceid)
 
    end subroutine checkfacesblock
    subroutine check1to1subfaces
        !
        !       check1to1Subfaces checks if the 1 to 1 internal subfaces,
        !       including the periodic ones, match up to a certain tolerance.
        !       If not, a warning will be printed. The computation is not
        !       returnFaild, because sometimes gaps are introduced on purpose,
        !       e.g. near a wing tip in an H-topology in spanwise direction.
        !
        use constants
        use blockpointers, only: ndom, nbocos, inbeg, inend, jnbeg, &
                                 jnend, knbeg, knend, n1to1, nbkglobal, neighproc, dinend, &
                                 dinbeg, djnbeg, djnend, dknbeg, dknend, x, cgnssubface, &
                                 neighblock, l1, l2, l3
        use cgnsgrid, only: cgnsdoms, cgnsndom
        use communication, only: myid, adflow_comm_world, nproc, sendrequests, &
                                 recvrequests
        use inputphysics
        use inputtimespectral, only: ntimeintervalsspectral
        use utils, only: delta, setpointers, terminate
        use partitionmod, only: sortbadentities
        use commonformats, only: strings
        implicit none
        !
        !      Local variables.
        !
        integer :: ierr, procId, size
 
        integer, dimension(mpi_status_size) :: mpiStatus
        integer, dimension(nProc) :: sizeMessage
        integer, dimension(nProc) :: recvcounts, displs
 
        integer(kind=intType) :: i, j, k, ll1, ll2, ll3, ic, jc, kc
        integer(kind=intType) :: sps, nn, mm, ll, ii, proc, nFCheck
        integer(kind=intType) :: stepI, stepJ, stepK
        integer(kind=intType) :: nMessagesSend, nMessagesReceive
        integer(kind=intType) :: nBad, nBadGlobal
 
        integer(kind=intType), dimension(3, 3) :: trMat
 
        integer(kind=intType), dimension(0:nProc) :: nFSend, nCoor
        integer(kind=intType), dimension(nProc) :: nFCount, nCCount
 
        integer(kind=intType), dimension(:, :), allocatable :: intBuf
        integer(kind=intType), dimension(:, :), allocatable :: intRecv
        integer(kind=intType), dimension(:, :), allocatable :: badSubfaces
        integer(kind=intType), dimension(:, :), allocatable :: badGlobal
 
        real(kind=realtype), dimension(:), allocatable :: baddist
        real(kind=realtype), dimension(:), allocatable :: baddistglobal
        real(kind=realtype), dimension(:, :), allocatable :: realbuf
        real(kind=realtype), dimension(:, :), allocatable :: realrecv
 
        character(len=7) :: intString
        character(len=maxStringLen) :: devFormat, spectralDevFormat
        !
        !       Determine the local number of faces that must be sent to other
        !       processors, including to myself. Also determine the number of
        !       faces I have to check.
        !
        nfcheck = 0
        nfsend = 0
        ncoor = 0
 
        ! Loop over the local blocks and its subfaces of the number of
        ! spectral solutions to be checked.
 
        do sps = 1, ntimeintervalsspectral
            do nn = 1, ndom
 
                call setpointers(nn, 1_inttype, sps)
 
                ! Loop over the 1 to 1 subfaces.
 
                do mm = 1, n1to1
 
                    ! Add the offset of nBocos to mm, such that the entries
                    ! in the arrays corresponds to this 1 to 1 subface.
 
                    ll = mm + nbocos
 
                    ! Update nFSend and nCoor of the correct processor and
                    ! update nFCheck.
 
                    ii = (abs(inend(ll) - inbeg(ll)) + 1) &
                         * (abs(jnend(ll) - jnbeg(ll)) + 1) &
                         * (abs(knend(ll) - knbeg(ll)) + 1)
                    ll = neighproc(ll) + 1
 
                    nfsend(ll) = nfsend(ll) + 1
                    ncoor(ll) = ncoor(ll) + ii
 
                    nfcheck = nfcheck + 1
 
                end do
            end do
        end do
 
        ! Put nFSend and nCoor in cumulative storage format. Note
        ! that nFSend and nCoor start at index 0. Store the starting
        ! value in nFCount and nCCount respectively.
 
        do nn = 1, nproc
            nfcount(nn) = nfsend(nn - 1)
            nccount(nn) = ncoor(nn - 1)
 
            nfsend(nn) = nfsend(nn) + nfsend(nn - 1)
            ncoor(nn) = ncoor(nn) + ncoor(nn - 1)
        end do
        !
        !       Determine the integer and real buffers to store the subface
        !       information to be communicated.
        !
        ! Allocate the memory for the integer and real buffers to store
        ! the information of the subfaces to be communicated.
 
        nn = nfsend(nproc)
        mm = ncoor(nproc)
 
        allocate (intbuf(10, nn), realbuf(3, mm), stat=ierr)
        if (ierr /= 0) &
             call terminate("check1to1Subfaces", &
             "Memory allocation failure for intBuf and &
             &realBuf")
 
        ! Repeat the loop over the number of local blocks and its subfaces
        ! of the number of spectral solutions to be checked.
 
        do sps = 1, ntimeintervalsspectral
            do nn = 1, ndom
 
                call setpointers(nn, 1_inttype, sps)
 
                ! Loop over the 1 to 1 subfaces.
 
                do mm = 1, n1to1
 
                    ! Add the offset of nBocos to mm, such that the entries
                    ! in the arrays corresponds to this 1 to 1 subface.
 
                    ll = mm + nbocos
 
                    ! Store the donor range, the local block number, the
                    ! spectral solution and global block ID in intBuf.
 
                    proc = neighproc(ll) + 1
                    nfcount(proc) = nfcount(proc) + 1
                    ii = nfcount(proc)
 
                    intbuf(1, ii) = dinbeg(ll)
                    intbuf(2, ii) = dinend(ll)
                    intbuf(3, ii) = djnbeg(ll)
                    intbuf(4, ii) = djnend(ll)
                    intbuf(5, ii) = dknbeg(ll)
                    intbuf(6, ii) = dknend(ll)
 
                    intbuf(7, ii) = neighblock(ll)
                    intbuf(8, ii) = sps
 
                    intbuf(9, ii) = nbkglobal
                    intbuf(10, ii) = cgnssubface(ll)
 
                    ! Determine whether the subface has positive or negative
                    ! running indices. To be sure that the correct sequence is
                    ! stored in realBuf, the loop must be performed over the
                    ! donor range (i, j, and k indices could be swapped and
                    ! therefore stored wrongly in the 1D buffer).
 
                    stepi = 1; if (dinend(ll) < dinbeg(ll)) stepi = -1
                    stepj = 1; if (djnend(ll) < djnbeg(ll)) stepj = -1
                    stepk = 1; if (dknend(ll) < dknbeg(ll)) stepk = -1
 
                    ! Determine the transformation matrix between the donor
                    ! and the current face. As the information stored in l1,
                    ! l2 and l3 is for the transformation matrix to the donor
                    ! face, the transpose must be taken.
 
                    ll1 = l1(ll); ll2 = l2(ll); ll3 = l3(ll)
 
                    trmat(1, 1) = sign(1_inttype, ll1) * delta(ll1, 1_inttype)
                    trmat(1, 2) = sign(1_inttype, ll1) * delta(ll1, 2_inttype)
                    trmat(1, 3) = sign(1_inttype, ll1) * delta(ll1, 3_inttype)
 
                    trmat(2, 1) = sign(1_inttype, ll2) * delta(ll2, 1_inttype)
                    trmat(2, 2) = sign(1_inttype, ll2) * delta(ll2, 2_inttype)
                    trmat(2, 3) = sign(1_inttype, ll2) * delta(ll2, 3_inttype)
 
                    trmat(3, 1) = sign(1_inttype, ll3) * delta(ll3, 1_inttype)
                    trmat(3, 2) = sign(1_inttype, ll3) * delta(ll3, 2_inttype)
                    trmat(3, 3) = sign(1_inttype, ll3) * delta(ll3, 3_inttype)
 
                    ! Store the coordinates in realBuf by looping over the
                    ! points of the subface.
 
                    ii = nccount(proc)
 
                    do k = dknbeg(ll), dknend(ll), stepk
                        do j = djnbeg(ll), djnend(ll), stepj
                            do i = dinbeg(ll), dinend(ll), stepi
 
                                ! Determine the nodal indices in the current block.
 
                                ll1 = i - dinbeg(ll)
                                ll2 = j - djnbeg(ll)
                                ll3 = k - dknbeg(ll)
 
                                ic = inbeg(ll) &
                                     + trmat(1, 1) * ll1 + trmat(1, 2) * ll2 + trmat(1, 3) * ll3
                                jc = jnbeg(ll) &
                                     + trmat(2, 1) * ll1 + trmat(2, 2) * ll2 + trmat(2, 3) * ll3
                                kc = knbeg(ll) &
                                     + trmat(3, 1) * ll1 + trmat(3, 2) * ll2 + trmat(3, 3) * ll3
 
                                ! Store the coordinates in the buffer.
 
                                ii = ii + 1
                                realbuf(1, ii) = x(ic, jc, kc, 1)
                                realbuf(2, ii) = x(ic, jc, kc, 2)
                                realbuf(3, ii) = x(ic, jc, kc, 3)
                            end do
                        end do
                    end do
 
                    ! Set ii to the number of points stored for this subface.
                    ! This number is needed when for the periodic correction.
 
                    ii = ii - nccount(proc)
 
                    ! Check if a periodic correction must be applied and if so
                    ! call the routine to do so for the coordinates of this
                    ! subface. Note that internally created subfaces due to
                    ! block splitting must be excluded from this test.
 
                    k = cgnssubface(ll)
                    if (k > 0) then
                        if (cgnsdoms(nbkglobal)%conn1to1(k)%periodic) then
 
                            j = nccount(proc) + 1
                            call periodictransformsubface(realbuf(1, j), ii, &
                                                          cgnsdoms(nbkglobal)%conn1to1(k)%rotationCenter, &
                                                          cgnsdoms(nbkglobal)%conn1to1(k)%rotationAngles, &
                                                          cgnsdoms(nbkglobal)%conn1to1(k)%translation)
                        end if
                    end if
 
                    ! Update the counter nCCount(proc).
 
                    nccount(proc) = nccount(proc) + ii
 
                end do
            end do
        end do
        !
        !       Determine the number of messages I will send and receive.
        !       The term message in this routine means a the complete subface
        !       info, i.e. a combination of an integer and real message.
        !
        ! Fill the array nCCount with 0's and 1's. A 0 indicates that no
        ! message is sent to the corresponding processor.
 
        do nn = 1, nproc
            nccount(nn) = 0
            if (nfsend(nn) > nfsend(nn - 1)) nccount(nn) = 1
        end do
 
        ! Make sure that no message is sent to myself. An offset of + 1
        ! must be added, because the processor numbers start at 0 and
        ! nCCount at 1.
 
        nccount(myid + 1) = 0
 
        ! Determine the number of message I have to receive.
 
        sizemessage = 1
        call mpi_reduce_scatter(nccount, nmessagesreceive, &
                                sizemessage, adflow_integer, mpi_sum, &
                                adflow_comm_world, ierr)
 
        ! Send the data I have to send. Do not send a message to myself.
        ! That is handled separately. As nonblocking sends must be used
        ! to avoid deadlock and two messages are sent to a processor,
        ! the sendRequests (stored in the module communication) are used
        ! for the integer messages and the recvRequests (same module)
        ! are used for the real messages.
 
        ii = 0
        do nn = 1, nproc
 
            ! Check if something must be sent to this processor.
 
            if (nccount(nn) == 1) then
 
                ! Send the integer buffer.
 
                ii = ii + 1
                procid = nn - 1
                size = 10 * (nfsend(nn) - nfsend(nn - 1))
                mm = nfsend(nn - 1) + 1
 
                call mpi_isend(intbuf(1, mm), size, adflow_integer, procid, &
                               procid, adflow_comm_world, sendrequests(ii), &
                               ierr)
 
                ! Send the real buffer.
 
                size = 3 * (ncoor(nn) - ncoor(nn - 1))
                mm = ncoor(nn - 1) + 1
 
                call mpi_isend(realbuf(1, mm), size, adflow_real, procid, &
                               procid + 1, adflow_comm_world, &
                               recvrequests(ii), ierr)
 
            end if
        end do
 
        ! Store the number of messages sent.
 
        nmessagessend = ii
        !
        !       Check the coordinates of the subfaces which should have been
        !       sent to myself.
        !
        ! Initialize nBad to 0 and allocate the memory to store possible
        ! bad subfaces.
 
        nbad = 0
        allocate (badsubfaces(4, nfcheck), baddist(nfcheck), stat=ierr)
        if (ierr /= 0) &
             call terminate("check1to1Subfaces", &
             "Memory allocation failure for badSubfaces &
             &and badDist")
 
        ! Determine the number of local subfaces that must be checked.
        ! If any present, check them.
 
        nn = nfsend(myid + 1) - nfsend(myid)
        if (nn > 0) then
            ii = nfsend(myid) + 1
            mm = ncoor(myid) + 1
            call checksubfacecoor(intbuf(1, ii), realbuf(1, mm), nn, &
                                  nbad, badsubfaces, baddist, &
                                  ntimeintervalsspectral)
        end if
        !
        !       Check the coordinates of the subfaces which are received from
        !       other processors.
        !
        ! Loop over the number of messages to be received.
 
        do ii = 1, nmessagesreceive
 
            ! Wait until an integer message arrives and determine the
            ! source and size of the message.
 
            call mpi_probe(mpi_any_source, myid, adflow_comm_world, &
                           mpistatus, ierr)
 
            procid = mpistatus(mpi_source)
            call mpi_get_count(mpistatus, adflow_integer, size, ierr)
 
            ! Check in debug mode that the incoming message is of
            ! correct size.
 
            if (debug) then
                if (size == mpi_undefined .or. mod(size, 10) /= 0) &
                    call terminate("check1to1Subfaces", &
                                   "Unexpected size of integer message")
            end if
 
            ! Determine the number of subfaces this message contains and
            ! allocate the memory for the integer receive buffer.
 
            nn = size / 10
            allocate (intrecv(10, nn), stat=ierr)
            if (ierr /= 0) &
                call terminate("check1to1Subfaces", &
                               "Memory allocation failure for intRecv")
 
            ! Receive the integer buffer. Blocking receives can be used,
            ! because the message has already arrived.
 
            call mpi_recv(intrecv, size, adflow_integer, procid, &
                          myid, adflow_comm_world, mpistatus, ierr)
 
            ! Probe for the corresponding real buffer and determine its
            ! size.
 
            call mpi_probe(procid, myid + 1, adflow_comm_world, &
                           mpistatus, ierr)
            call mpi_get_count(mpistatus, adflow_real, size, ierr)
 
            ! Check in debug mode that the incoming message is of
            ! correct size.
 
            if (debug) then
                if (size == mpi_undefined .or. mod(size, 3) /= 0) &
                    call terminate("check1to1Subfaces", &
                                   "Unexpected size of real message")
            end if
 
            ! Determine the total number of coordinates in the message and
            ! allocate the memory for the real receive buffer.
 
            mm = size / 3
            allocate (realrecv(3, mm), stat=ierr)
            if (ierr /= 0) &
                call terminate("check1to1Subfaces", &
                               "Memory allocation failure for realRecv")
 
            ! Receive the real buffer. Blocking receives can be used,
            ! because the message has already arrived.
 
            call mpi_recv(realrecv, size, adflow_real, procid, &
                          myid + 1, adflow_comm_world, mpistatus, ierr)
 
            ! Check the subfaces stored in these messages.
 
            call checksubfacecoor(intrecv, realrecv, nn, nbad, &
                                  badsubfaces, baddist, &
                                  ntimeintervalsspectral)
 
            ! Release the memory of the integer and real receive buffer.
 
            deallocate (intrecv, realrecv, stat=ierr)
            if (ierr /= 0) &
                 call terminate("check1to1Subfaces", &
                 "Deallocation failure for intRecv &
                 &and realRecv")
        end do
 
        ! Complete the nonblocking sends.
 
        size = nmessagessend
        do nn = 1, nmessagessend
            call mpi_waitany(size, sendrequests, procid, mpistatus, ierr)
            call mpi_waitany(size, recvrequests, procid, mpistatus, ierr)
        end do
 
        ! Deallocate the memory for the integer and real buffers.
 
        deallocate (intbuf, realbuf, stat=ierr)
        if (ierr /= 0) &
            call terminate("check1to1Subfaces", &
                           "Deallocation failure for intBuf and realBuf")
        !
        !       Determine the global number of bad subfaces and gather this
        !       information.
        !
        ! Determine the global number of bad subfaces.
 
        call mpi_allgather(nbad, 1, adflow_integer, nccount, 1, &
                           adflow_integer, adflow_comm_world, ierr)
 
        ! Determine the global number of bad subfaces and the arrays
        ! recvcounts and displs needed for the call to allgatherv
 
        nbadglobal = nccount(1)
        recvcounts(1) = nccount(1)
        displs(1) = 0
 
        do nn = 2, nproc
            nbadglobal = nbadglobal + nccount(nn)
            recvcounts(nn) = nccount(nn)
            displs(nn) = displs(nn - 1) + recvcounts(nn - 1)
        end do
 
        ! Allocate the memory to store the global bad surfaces.
 
        allocate (badglobal(4, nbadglobal), baddistglobal(nbadglobal), &
                  stat=ierr)
        if (ierr /= 0) &
             call terminate("check1to1Subfaces", &
             "Memory allocation failure for badGlobal &
             &and badDistGlobal")
 
        ! Gather the data. First the distance info.
 
        size = nbad
        call mpi_allgatherv(baddist, size, adflow_real, baddistglobal, &
                            recvcounts, displs, adflow_real, &
                            adflow_comm_world, ierr)
 
        ! And the integer info. Multiply recvcounts and displs
        ! by 4, because 4 integers are received.
 
        do nn = 1, nproc
            recvcounts(nn) = 4 * recvcounts(nn)
            displs(nn) = 4 * displs(nn)
        end do
 
        size = 4 * nbad
        call mpi_allgatherv(badsubfaces, size, adflow_integer, badglobal, &
                            recvcounts, displs, adflow_integer, &
                            adflow_comm_world, ierr)
 
        ! Sort the bad subfaces and get rid of the multiple entries.
 
        call sortbadentities(nbadglobal, badglobal, baddistglobal, .true.)
 
        ! Check for the presence of any internally created subfaces.
        ! This only occurs when something goes wrong in the block
        ! splitting and therefore the program is terminated.
 
        do nn = 1, nbadglobal
            if (badglobal(2, nn) == 0) then
                if (myid == 0) &
                     call terminate("check1to1Subfaces", &
                     "Non-matching internally created &
                     &face found.")
                call mpi_barrier(adflow_comm_world, ierr)
            end if
        end do
 
        ! Print the bad subfaces, if present. Only processor 0 performs
        ! this task.
 
        devformat = "(7(A), ES12.5)"
        spectraldevformat = "(9(A), ES12.5)"
 
        if (myid == 0 .and. nbadglobal > 0) then
 
            write (intstring, "(i6)") nbadglobal
            intstring = adjustl(intstring)
 
            print "(a)", "#"
            print strings, "#                      Warning"
            print strings, "# Found ", trim(intstring), " one to one subfaces which do not coincide."
            print strings, "# Computation continues, but be aware of it."
            print strings, "# List of nonmatching one to one subfaces."
            print "(a)", "#"
 
            do nn = 1, nbadglobal
                i = badglobal(1, nn)
                j = badglobal(2, nn)
 
                write (intstring, "(i6)") badglobal(3, nn)
                intstring = adjustl(intstring)
 
                ! Write a different error message if more than one grid has
                ! been read.
 
                if (ntimeintervalsspectral > 1) then
 
                    if (badglobal(4, nn) == 1) then
                        print spectraldevformat, "# Spectral grid ", trim(intstring), &
                            ", zone ", trim(cgnsdoms(i)%zoneName), &
                            ", periodic subface ", trim(cgnsdoms(i)%conn1to1(j)%connectName), &
                            " does not match donor ", trim(cgnsdoms(i)%conn1to1(j)%donorName), &
                            ". Maximum deviation: ", baddistglobal(nn)
                    else
                        print spectraldevformat, "# Spectral grid ", trim(intstring), &
                            ", zone ", trim(cgnsdoms(i)%zoneName), &
                            ", subface ", trim(cgnsdoms(i)%conn1to1(j)%connectName), &
                            " does not match donor ", trim(cgnsdoms(i)%conn1to1(j)%donorName), &
                            ". Maximum deviation: ", baddistglobal(nn)
                    end if
 
                else
 
                    if (badglobal(4, nn) == 1) then
                        print devformat, "Zone ", trim(cgnsdoms(i)%zoneName), &
                            ", periodic subface ", trim(cgnsdoms(i)%conn1to1(j)%connectName), &
                            " does not match donor ", trim(cgnsdoms(i)%conn1to1(j)%donorName), &
                            ". Maximum deviation: ", baddistglobal(nn)
                    else
                        print devformat, "Zone ", trim(cgnsdoms(i)%zoneName), &
                            ", subface ", trim(cgnsdoms(i)%conn1to1(j)%connectName), &
                            " does not match donor ", trim(cgnsdoms(i)%conn1to1(j)%donorName), &
                            ". Maximum deviation: ", baddistglobal(nn)
                    end if
 
                end if
            end do
 
            print "(a)", "#"
 
        end if
 
        ! Deallocate the memory of to store the bad subfaces.
 
        deallocate (badsubfaces, badglobal, baddist, baddistglobal, &
                    stat=ierr)
        if (ierr /= 0) &
             call terminate("check1to1Subfaces", &
             "Deallocation failure for badSubfaces, &
             &badGlobal, badDist and badDistGlobal")
 
        ! Synchronize the processors, because wild cards have been used.
 
        call mpi_barrier(adflow_comm_world, ierr)
 
    end subroutine check1to1subfaces
 
    !      ==================================================================
 
    subroutine periodictransformsubface(coor, nn, rotCenter, &
                                        rotAngles, translation)
        !
        !       periodicTransformSubface transforms the given set of
        !       coordinates using the periodic transformation defined by
        !       rotCenter, rotAngles and translation.
        !
        use constants
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: nn
 
        real(kind=realtype), dimension(3), intent(in) :: rotcenter
        real(kind=realtype), dimension(3), intent(in) :: rotangles
        real(kind=realtype), dimension(3), intent(in) :: translation
 
        real(kind=realtype), dimension(3, nn), intent(inout) :: coor
        !
        !      Local variables.
        !
        integer(kind=intType) :: i
 
        real(kind=realtype) :: costheta, cosphi, cospsi
        real(kind=realtype) :: sintheta, sinphi, sinpsi
        real(kind=realtype) :: dx, dy, dz
 
        real(kind=realtype), dimension(3) :: trans
        real(kind=realtype), dimension(3, 3) :: rotmatrix
 
        ! Construct from the given rotation angles the rotation matrix
        ! from the current coordinates to the donor coordinates.
        ! Note that the sequence of rotation is first rotation around the
        ! x-axis, followed by rotation around the y-axis and finally
        ! rotation around the z-axis.
 
        costheta = cos(rotangles(1)); sintheta = sin(rotangles(1))
        cosphi = cos(rotangles(2)); sinphi = sin(rotangles(2))
        cospsi = cos(rotangles(3)); sinpsi = sin(rotangles(3))
 
        rotmatrix(1, 1) = cosphi * cospsi
        rotmatrix(2, 1) = cosphi * sinpsi
        rotmatrix(3, 1) = -sinphi
 
        rotmatrix(1, 2) = sintheta * sinphi * cospsi - costheta * sinpsi
        rotmatrix(2, 2) = sintheta * sinphi * sinpsi + costheta * cospsi
        rotmatrix(3, 2) = sintheta * cosphi
 
        rotmatrix(1, 3) = costheta * sinphi * cospsi + sintheta * sinpsi
        rotmatrix(2, 3) = costheta * sinphi * sinpsi - sintheta * cospsi
        rotmatrix(3, 3) = costheta * cosphi
 
        ! Store the translation plus the rotation center in trans.
 
        trans(1) = rotcenter(1) + translation(1)
        trans(2) = rotcenter(2) + translation(2)
        trans(3) = rotcenter(3) + translation(3)
 
        ! Loop over the number of coordinates to be corrected.
 
        do i = 1, nn
 
            ! Determine the relative position w.R.T. The rotation center.
 
            dx = coor(1, i) - rotcenter(1)
            dy = coor(2, i) - rotcenter(2)
            dz = coor(3, i) - rotcenter(3)
 
            ! Determine the coordinates after the transformation.
 
            coor(1, i) = rotmatrix(1, 1) * dx + rotmatrix(1, 2) * dy &
                         + rotmatrix(1, 3) * dz + trans(1)
            coor(2, i) = rotmatrix(2, 1) * dx + rotmatrix(2, 2) * dy &
                         + rotmatrix(2, 3) * dz + trans(2)
            coor(3, i) = rotmatrix(3, 1) * dx + rotmatrix(3, 2) * dy &
                         + rotmatrix(3, 3) * dz + trans(3)
 
        end do
 
    end subroutine periodictransformsubface
 
    !      ==================================================================
 
    subroutine checksubfacecoor(subfaceInfo, coor, nFace, &
                                nBad, badSubfaces, badDist, &
                                nSpectral)
        !
        !       checkSubfaceCoor checks if the coordinates of the subfaces
        !       defined in subfaceInfo and coor match the coordinates stored
        !       in flowDoms.
        !
        use constants
        use blockpointers
        use cgnsgrid
        use communication
        use utils, only: setpointers, terminate
        use commonformats, only: strings
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: nFace, nSpectral
        integer(kind=intType), intent(inout) :: nBad
 
        integer(kind=intType), dimension(10, *), intent(in) :: subfaceInfo
        integer(kind=intType), dimension(4, *), intent(inout) :: badSubfaces
 
        real(kind=realtype), dimension(*), intent(inout) :: baddist
        real(kind=realtype), dimension(3, *), intent(in) :: coor
        !
        !      Local parameter.
        !
        real(kind=realtype), parameter :: reltol = 0.05_realtype
        !
        !      Local variables.
        !
        integer(kind=intType) :: i, j, k, ii, jj, mm, nn
        integer(kind=intType) :: stepI, stepJ, stepK
        integer(kind=intType) :: im1, ip1, jm1, jp1, km1, kp1
        integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, kBeg, kEnd
        integer(kind=intType) :: blockID, sps, globalDonID, subfaceID
 
        real(kind=realtype) :: dist2, dist2i, dist2j, dist2k, tol
        real(kind=realtype) :: facti, factj, factk, diffmax
 
        character(len=maxStringLen) :: errorMessage
        character(len=7) :: intString
 
        logical :: badFace
 
        ! Initialize the counter ii for the coordinates and loop over the
        ! number of subfaces to be checked.
 
        ii = 0
        subfaceloop: do nn = 1, nface
 
            ! Store the integer info a bit easier.
 
            ibeg = subfaceinfo(1, nn)
            iend = subfaceinfo(2, nn)
            jbeg = subfaceinfo(3, nn)
            jend = subfaceinfo(4, nn)
            kbeg = subfaceinfo(5, nn)
            kend = subfaceinfo(6, nn)
 
            blockid = subfaceinfo(7, nn)
            sps = subfaceinfo(8, nn)
            globaldonid = subfaceinfo(9, nn)
            subfaceid = subfaceinfo(10, nn)
 
            ! Set the pointers to the block to be searched.
 
            call setpointers(blockid, 1_inttype, sps)
 
            ! Find the matching 1 to 1 subface.
 
            do mm = 1, n1to1
 
                ! Add the offset of nBocos to nn such that it contains the
                ! correct index of the arrays.
 
                jj = mm + nbocos
 
                ! If this is the subface exit the loop.
 
                if (min(ibeg, iend) == min(inbeg(jj), inend(jj)) .and. &
                    max(ibeg, iend) == max(inbeg(jj), inend(jj)) .and. &
                    min(jbeg, jend) == min(jnbeg(jj), jnend(jj)) .and. &
                    max(jbeg, jend) == max(jnbeg(jj), jnend(jj)) .and. &
                    min(kbeg, kend) == min(knbeg(jj), knend(jj)) .and. &
                    max(kbeg, kend) == max(knbeg(jj), knend(jj))) exit
 
            end do
 
            ! Check if a subface was found. If not terminate.
 
            if (mm > n1to1) then
 
                if (nspectral > 1) then
                    write (intstring, "(i6)") sps
                    intstring = adjustl(intstring)
                    write (errormessage, strings) "Spectral grid ", trim(intstring), &
                        ", Zone ", trim(cgnsdoms(globaldonid)%zoneName), &
                        "Connectivity ", trim(cgnsdoms(globaldonid)%conn1to1(subfaceid)%connectName), &
                        ": 1 to 1 subface not found. Something is seriously wrong with the zone connectivity."
 
                else
                    write (errormessage, strings) "Zone ", trim(cgnsdoms(globaldonid)%zoneName), &
                        "Connectivity", trim(cgnsdoms(globaldonid)%conn1to1(subfaceid)%connectName), &
                        ": 1 to 1 subface not found. Something is seriously wrong with the zone connectivity."
                end if
 
                call terminate("checkSubfaceCoor", errormessage)
 
            end if
 
            ! Determine whether positive or negative running indices must
            ! be used for the subface. This depends on the sequence stored
            ! in the coordinate buffer and not of the sequence of the
            ! 1 to 1 subface mm.
 
            stepi = 1; if (iend < ibeg) stepi = -1
            stepj = 1; if (jend < jbeg) stepj = -1
            stepk = 1; if (kend < kbeg) stepk = -1
 
            ! Initialize badFace to .false. to indicate that it is
            ! a correct subface. Set the maximum difference to zero.
 
            badface = .false.
            diffmax = zero
 
            ! Loop over the coordinates of the subface to see if they
            ! match of to a certain tolerance.
 
            do k = kbeg, kend, stepk
 
                ! Determine the indices to the left and to the right.
                ! Do not exceed the block boundary. Determine the
                ! scaling factor of the tolerance accordingly.
 
                km1 = max(1_inttype, k - 1_inttype)
                kp1 = min(kl, k + 1_inttype)
 
                factk = one
                if (km1 == k) factk = factk * four
                if (kp1 == k) factk = factk * four
 
                ! Loop in j-direction.
 
                do j = jbeg, jend, stepj
 
                    ! Determine the neighbors to the left and right and factJ.
 
                    jm1 = max(1_inttype, j - 1_inttype)
                    jp1 = min(jl, j + 1_inttype)
 
                    factj = one
                    if (jm1 == j) factj = factj * four
                    if (jp1 == j) factj = factj * four
 
                    ! Loop in i-direction.
 
                    do i = ibeg, iend, stepi
 
                        ! Determine the neighbors to the left and right and factI.
 
                        im1 = max(1_inttype, i - 1_inttype)
                        ip1 = min(il, i + 1_inttype)
 
                        facti = one
                        if (im1 == i) facti = facti * four
                        if (ip1 == i) facti = facti * four
 
                        ! Determine the distances squared in i, j and k direction.
                        ! The reason for squared is that some square roots are
                        ! avoided this way.
 
                        dist2i = (x(ip1, j, k, 1) - x(im1, j, k, 1))**2 &
                                 + (x(ip1, j, k, 2) - x(im1, j, k, 2))**2 &
                                 + (x(ip1, j, k, 3) - x(im1, j, k, 3))**2
 
                        dist2j = (x(i, jp1, k, 1) - x(i, jm1, k, 1))**2 &
                                 + (x(i, jp1, k, 2) - x(i, jm1, k, 2))**2 &
                                 + (x(i, jp1, k, 3) - x(i, jm1, k, 3))**2
 
                        dist2k = (x(i, j, kp1, 1) - x(i, j, km1, 1))**2 &
                                 + (x(i, j, kp1, 2) - x(i, j, km1, 2))**2 &
                                 + (x(i, j, kp1, 3) - x(i, j, km1, 3))**2
 
                        ! Make sure that singular lines are excluded.
 
                        if (dist2i < eps) dist2i = large
                        if (dist2j < eps) dist2j = large
                        if (dist2k < eps) dist2k = large
 
                        ! Determine the tolerance for identical points.
                        ! Note the multiplication with the square of the relative
                        ! tolerance, because the distances are squared as well.
 
                        tol = factk * dist2k
                        tol = min(tol, factj * dist2j)
                        tol = min(tol, facti * dist2i)
 
                        tol = tol * reltol * reltol
 
                        ! Update the counter for the coordinate in the buffer
                        ! and determine the distance squared between the points
                        ! that should be identical.
 
                        ii = ii + 1
                        dist2 = (x(i, j, k, 1) - coor(1, ii))**2 &
                                + (x(i, j, k, 2) - coor(2, ii))**2 &
                                + (x(i, j, k, 3) - coor(3, ii))**2
 
                        ! Flag the subface to bad if the nodes do not coincide
                        ! within the given tolerance. Store the difference if
                        ! it is larger than the currently stored value.
 
                        if (dist2 > tol) then
                            badface = .true.
                            diffmax = max(diffmax, dist2)
                        end if
 
                    end do
                end do
            end do
 
            ! Store this subface if it is not matching. The following data
            ! is stored: global block number of the donor, the subface ID
            ! on this block, the spectral solution, whether or not this
            ! is a periodic face and the maximum difference.
 
            if (badface) then
                nbad = nbad + 1
 
                badsubfaces(1, nbad) = globaldonid
                badsubfaces(2, nbad) = subfaceid
                badsubfaces(3, nbad) = sps
 
                badsubfaces(4, nbad) = 0
                if (subfaceid > 0) then
                    if (cgnsdoms(globaldonid)%conn1to1(subfaceid)%periodic) &
                        badsubfaces(4, nbad) = 1
                end if
 
                baddist(nbad) = sqrt(diffmax)
 
                !  if(subfaceID == 0 .and. myID == 1) then
                !    write(*,*) "myID: ", myID, badDist(nBad)
                !    write(*,*) blockID, mm, globalDonID, nbkGlobal
                !    write(*,"(6I4)") iBegor, iEndor, jBegor, jEndor, kBegor, kEndor
 
                !    jj = mm + nBocos
                !    write(*,"(6I4)") inBeg(jj),inEnd(jj), jnBeg(jj),jnEnd(jj), knBeg(jj),knEnd(jj)
                !  endif
            end if
 
        end do subfaceloop
 
    end subroutine checksubfacecoor
end module gridchecking