partitionMod.F90

Go to the documentation of this file.

module partitionmod
    !
    !       This local module contains definitions of derived datatypes
    !       as well as variables used in the partitioning directory.
    !
    use constants
    implicit none
    save
    !
    !       The definition of the derived datatype distributionBlockType
    !
    type distributionblocktype
        !
        !         Block dimensions and local block ID.
        !
        !  nx, ny, nz - block integer dimensions for no halo cell based
        !               quantities.
        !  il, jl, kl - block integer dimensions for no halo node based
        !               quantities.
        !  blockID    - local block ID on the processor this block is
        !               stored.
 
        integer(kind=intType) :: nx, ny, nz, &
                                 il, jl, kl
        integer(kind=intType) :: blockid
        !
        !         Total number cells and faces inside the block. In the number
        !         faces the work for nonmatching block boundaries is included,
        !         such that the load balance is still guaranteed.
        !
        ! Ncell     : total number of cells in this block.
        ! Nface     : total number of faces in this block.
 
        integer(kind=intType) :: ncell, nface
        !
        !         Block boundary conditions.
        !
        !  nSubface             - Number of subfaces on this block.
        !  n1to1                - Number of 1 to 1 block boundaries.
        !  nBocos               - Number of physical boundary subfaces.
        !  BCType(:)            - Boundary condition type for each
        !                         subface. See the module BCTypes for
        !                         the possibilities.
        !  BCFaceID(:)          - Block face location of each subface.
        !                         Possible values are: iMin, iMax, jMin,
        !                         jMax, kMin, kMax.
        !  cgnsSubface(:)       - The subface in the corresponding cgns
        !                         block. As cgns distinguishes between
        !                         boundary and internal boundaries, the
        !                         BCType of the subface is needed to
        !                         know which one to take. A zero
        !                         indicates that this face was obtained
        !                         by splitting a cgns block.
        !  inBeg(:), inEnd(:)   - Lower and upper limits for the nodes
        !  jnBeg(:), jnEnd(:)     in each of the index directions on a
        !  knBeg(:), knEnd(:)     given subface. Note that one of these
        !                         indices will not change since we will
        !                         be moving on a face.
        !  dinBeg(:), dinEnd(:) - Lower and upper limits for the nodes
        !  djnBeg(:), djnEnd(:)   in the each of the index directions
        !  dknBeg(:), dknEnd(:)   of the donor subface for this
        !                         particular subface. Note that one of
        !                         these indices will not change since we
        !                         will be moving on a face.
        !  neighBlock(:)        - Block number to which this subface
        !                         connects. This value is set to zero if
        !                         this subface is not connected to
        !                         another block.
        !  l1(:), l2(:),        - Short hand for the transformation
        !  l3(:)                  matrix between this subface and the
        !                         neighbor block. These value are set to
        !                         zero if this subface is not connected
        !                         to another block.
        !  groupNum(:)          - Group number to which this subface
        !                         belongs. If this subface does not
        !                         belong to any group, the corresponding
        !                         entry in this array is zeroed out.
 
        integer(kind=intType) :: nsubface, n1to1, nbocos
 
        integer(kind=intType), dimension(:), pointer :: bctype
        integer(kind=intType), dimension(:), pointer :: bcfaceid
        integer(kind=intType), dimension(:), pointer :: cgnssubface
 
        integer(kind=intType), dimension(:), pointer :: inbeg, inend
        integer(kind=intType), dimension(:), pointer :: jnbeg, jnend
        integer(kind=intType), dimension(:), pointer :: knbeg, knend
 
        integer(kind=intType), dimension(:), pointer :: dinbeg, dinend
        integer(kind=intType), dimension(:), pointer :: djnbeg, djnend
        integer(kind=intType), dimension(:), pointer :: dknbeg, dknend
 
        integer(kind=intType), dimension(:), pointer :: neighblock
        integer(kind=intType), dimension(:), pointer :: l1, l2, l3
        integer(kind=intType), dimension(:), pointer :: groupnum
 
        !
        !         Relation to the original cgns grid.
        !
        ! cgnsBlockID    - block/zone number of the cgns grid to which
        !                  this block is related.
        ! iBegOr, iEndOr - range of points of this block in the
        ! jBegOr, jEndOr   corresponding cgns block, i.e. for this block
        ! kBegOr, kEndOr   iBegOr <= i <= iEndOr, jBegOr <= j <= jEndOr,
        !                  kBegOr <= k <= kEndOr.
        !                  It is of course possible that the entire
        !                  block is stored.
 
        integer(kind=intType) :: cgnsblockid
        integer(kind=intType) :: ibegor, iendor, jbegor, jendor
        integer(kind=intType) :: kbegor, kendor
 
    end type distributionblocktype
 
    ! nBlocks        : Number of blocks to be distributed.
    ! blocks(nBlocks): The array with the block info.
 
    integer(kind=intType) :: nblocks
    type(distributionblocktype), dimension(:), allocatable :: blocks
 
    !
    !       Type definition to store the way the original cgns blocks are
    !       split for load balancing reasons.
    !
    type splitcgnstype
 
        ! nSubBlocks:             # of subblocks into which the original
        !                         block is split.
        ! ranges(nsubblocks,3,2): The nodal range in the original
        !                         block of these subblocks.
 
        integer(kind=intType) :: nsubblocks
        integer(kind=intType), dimension(:, :, :), pointer :: ranges
 
    end type splitcgnstype
 
    !
    !       Type definition needed to determine the processor ID's and
    !       nodal ranges of the subblocks for every CGNS block.
    !
    type subblocksofcgnstype
 
        ! cgnsBlockID    - block/zone number of the cgns grid to which
        !                  this subblock is related.
        ! procID         - Processor ID on which this subblock is
        !                  stored.
        ! blockID        - local block ID on the processor this block is
        !                  stored.
        ! iBegOr, iEndOr - range of points of this block in the
        ! jBegOr, jEndOr   corresponding cgns block, i.e. for this block
        ! kBegOr, kEndOr   iBegOr <= i <= iEndOr, jBegOr <= j <= jEndOr,
        !                  kBegOr <= k <= kEndOr.
        !                  It is of course possible that the entire
        !                  block is stored.
 
        integer :: cgnsblockid
        integer :: procid
        integer :: blockid
        integer :: ibegor, iendor, jbegor, jendor, kbegor, kendor
 
    end type subblocksofcgnstype
 
    ! Interface for the extension of the operators <= and <.
    ! These are needed for the sorting of subblocksOfCGNSType. Note
    ! that the = operator does not need to be defined, because
    ! subblocksOfCGNSType only contains primitive types.
 
    interface operator(<=)
        module procedure lessequalsubblocksofcgnstype
    end interface operator(<=)
 
    interface operator(<)
        module procedure lesssubblocksofcgnstype
    end interface operator(<)
    !
    !       Type definition needed to determine the number of distinct
    !       non-matching abutting subfaces in the CGNS file.
    !
    type subfacenonmatchtype
 
        ! iBeg, iEnd - Nodal subface range om the CGNS block, i-direction
        ! jBeg, jEnd - Idem in the j-direction
        ! kBeg, kEnd - Idem in the k-direction
        ! connID     - The cgns connectivity ID.
 
        integer :: ibeg, jbeg, kbeg, iend, jend, kend
        integer :: connid
 
    end type subfacenonmatchtype
 
    ! Interface for the extension of the operators <= and <.
    ! These are needed for the sorting of subfaceNonMatchType. Note
    ! that the = operator does not need to be defined, because
    ! subfaceNonMatchType only contains primitive types.
 
    interface operator(<=)
        module procedure lessequalsubfacenonmatchtype
    end interface operator(<=)
 
    interface operator(<)
        module procedure lesssubfacenonmatchtype
    end interface operator(<)
 
    type sortsubrangetype
 
        ! iMin: minimum i-index in the subrange.
        ! jMin: minimum j-index in the subrange.
        ! kMin: minimum k-index in the subrange.
        ! iMax: maximum i-index in the subrange.
        ! jMax: maximum j-index in the subrange.
        ! kMax: maximum k-index in the subrange.
 
        integer(kind=intType) :: imin, jmin, kmin
        integer(kind=intType) :: imax, jmax, kmax
 
    end type sortsubrangetype
 
    ! Interfaces for the definitions of the operators <=, < and /=.
    ! These are needed for the sorting of this derived data type.
    ! Note that the = operator does not need to be defined, because
    ! sortSubRangeType only contains primitive types.
 
    interface operator(<=)
        module procedure lessequalsortsubrangetype
    end interface operator(<=)
 
    interface operator(<)
        module procedure lesssortsubrangetype
    end interface operator(<)
 
    ! Definition of the derived data type.
 
    type fourintplusrealtype
        integer(kind=intType) :: n1, n2, n3, n4
        real(kind=realtype) :: dist
    end type fourintplusrealtype
 
    ! Interfaces for the definitions of the operators <=, < and /=.
    ! These are needed for the sorting of this derived data type.
    ! Note that the = operator does not need to be defined, because
    ! fourIntPlusRealType only contains primitive types.
 
    interface operator(<=)
        module procedure lessequalfourintplusrealtype
    end interface operator(<=)
 
    interface operator(<)
        module procedure lessfourintplusrealtype
    end interface operator(<)
 
    interface operator(/=)
        module procedure notequalfourintplusrealtype
    end interface operator(/=)
 
    ! ==========================================================================
    !
    !       Variable to store the partition number (processor ID) of the
    !       computational blocks.
    !
    ! ubvec(2):      Tolerance for the constraints.
    ! part(nBlocks): The processor ID for each block, starting at 0.
 
    real, dimension(2) :: ubvec
 
    integer(kind=intType), dimension(:), allocatable :: part
    !
    !       Variables needed for the reading of the grid files.
    !
    ! nGridsRead:            Number of grids to read.
    ! fileIDs(nGridsRead):   The file ID's.
    ! gridFiles(nGridsRead): Names of the grid files to read.
    ! interpolSpectral:      Whether or not to interpolate the
    !                        coordinates for the time spectral mode.
 
    integer(kind=intType) :: ngridsread
    logical :: interpolspectral
 
    integer, dimension(:), allocatable :: fileids
 
    character(len=maxStringLen), dimension(:), allocatable :: gridfiles
    ! ==========================================================================
 
contains
    !
    !         Functions to simulate the operators <= and < for the derived
    !         datatypes subblocksOfCGNSType and subfaceNonMatchType.
    !
    logical function lessequalsubblocksofcgnstype(g1, g2)
        !
        !         This function returns .true. if g1 <= g2 and .false.
        !         otherwise. The comparison is firstly based on the CGNS block
        !         ID, then the processor ID and finally the local block ID.
        !
        implicit none
        !
        !        Function arguments.
        !
        type(subblocksofcgnstype), intent(in) :: g1, g2
 
        ! Comparison of the CGNS block ID. If not equal, set
        ! lessEqualSubblocksOfCGNSType appropriately and return.
 
        if (g1%cgnsBlockID < g2%cgnsBlockID) then
            lessequalsubblocksofcgnstype = .true.
            return
        else if (g1%cgnsBlockID > g2%cgnsBlockID) then
            lessequalsubblocksofcgnstype = .false.
            return
        end if
 
        ! Compare the processor ID's.
 
        if (g1%procID < g2%procID) then
            lessequalsubblocksofcgnstype = .true.
            return
        else if (g1%procID > g2%procID) then
            lessequalsubblocksofcgnstype = .false.
            return
        end if
 
        ! Compare the local block ID's.
 
        if (g1%blockID < g2%blockID) then
            lessequalsubblocksofcgnstype = .true.
            return
        else if (g1%blockID > g2%blockID) then
            lessequalsubblocksofcgnstype = .false.
            return
        end if
 
        ! It does not make sense to compare the subranges, because
        ! these are identical if we came so far. Both entities are
        ! identical. So set lessEqualSubblocksOfCGNSType to .true.
 
        lessequalsubblocksofcgnstype = .true.
 
    end function lessequalsubblocksofcgnstype
 
    !===============================================================
 
    logical function lesssubblocksofcgnstype(g1, g2)
        !
        !         This function returns .true. if g1 < g2 and .false.
        !         otherwise. The comparison is firstly based on the CGNS block
        !         ID, then the processor ID and finally the local blockID.
        !
        implicit none
        !
        !        Function arguments.
        !
        type(subblocksofcgnstype), intent(in) :: g1, g2
 
        ! Comparison of the CGNS block ID. If not equal, set
        ! lessSubblocksOfCGNSType appropriately and return.
 
        if (g1%cgnsBlockID < g2%cgnsBlockID) then
            lesssubblocksofcgnstype = .true.
            return
        else if (g1%cgnsBlockID > g2%cgnsBlockID) then
            lesssubblocksofcgnstype = .false.
            return
        end if
 
        ! Compare the processor ID's.
 
        if (g1%procID < g2%procID) then
            lesssubblocksofcgnstype = .true.
            return
        else if (g1%procID > g2%procID) then
            lesssubblocksofcgnstype = .false.
            return
        end if
 
        ! Compare the local block ID's.
 
        if (g1%blockID < g2%blockID) then
            lesssubblocksofcgnstype = .true.
            return
        else if (g1%blockID > g2%blockID) then
            lesssubblocksofcgnstype = .false.
            return
        end if
 
        ! It does not make sense to compare the subranges, because
        ! these are identical if we came so far. Both entities are
        ! identical. So set lessSubblocksOfCGNSType to .false.
 
        lesssubblocksofcgnstype = .false.
 
    end function lesssubblocksofcgnstype
 
    !===============================================================
 
    logical function lessequalsubfacenonmatchtype(g1, g2)
        !
        !         This function returns .true. if g1 <= g2 and .false.
        !         otherwise. The comparison is firstly based on the i-range,
        !         followed by the j-range and k-range. If these are all the
        !         same the connectivity ID is compared.
        !
        implicit none
        !
        !        Function arguments.
        !
        type(subfacenonmatchtype), intent(in) :: g1, g2
 
        ! Comparison of the iBeg value. If different set
        ! lessEqualSubfaceNonMatchType appropriately and return.
 
        if (g1%iBeg < g2%iBeg) then
            lessequalsubfacenonmatchtype = .true.
            return
        else if (g1%iBeg > g2%iBeg) then
            lessequalsubfacenonmatchtype = .false.
            return
        end if
 
        ! The iEnd value.
 
        if (g1%iEnd < g2%iEnd) then
            lessequalsubfacenonmatchtype = .true.
            return
        else if (g1%iEnd > g2%iEnd) then
            lessequalsubfacenonmatchtype = .false.
            return
        end if
 
        ! The jBeg value.
 
        if (g1%jBeg < g2%jBeg) then
            lessequalsubfacenonmatchtype = .true.
            return
        else if (g1%jBeg > g2%jBeg) then
            lessequalsubfacenonmatchtype = .false.
            return
        end if
 
        ! The jEnd value.
 
        if (g1%jEnd < g2%jEnd) then
            lessequalsubfacenonmatchtype = .true.
            return
        else if (g1%jEnd > g2%jEnd) then
            lessequalsubfacenonmatchtype = .false.
            return
        end if
 
        ! The kBeg value.
 
        if (g1%kBeg < g2%kBeg) then
            lessequalsubfacenonmatchtype = .true.
            return
        else if (g1%kBeg > g2%kBeg) then
            lessequalsubfacenonmatchtype = .false.
            return
        end if
 
        ! The kEnd value.
 
        if (g1%kEnd < g2%kEnd) then
            lessequalsubfacenonmatchtype = .true.
            return
        else if (g1%kEnd > g2%kEnd) then
            lessequalsubfacenonmatchtype = .false.
            return
        end if
 
        ! The subrange is identical. Compare the connectivity ID.
 
        if (g1%connID < g2%connID) then
            lessequalsubfacenonmatchtype = .true.
            return
        else if (g1%connID > g2%connID) then
            lessequalsubfacenonmatchtype = .false.
            return
        end if
 
        ! g1 and g2 are identical. Return .true.
 
        lessequalsubfacenonmatchtype = .true.
 
    end function lessequalsubfacenonmatchtype
 
    !===============================================================
 
    logical function lesssubfacenonmatchtype(g1, g2)
        !
        !         This function returns .true. if g1 < g2 and .false.
        !         otherwise. The comparison is firstly based on the i-range,
        !         followed by the j-range and k-range. If these are all the
        !         same the connectivity ID is compared.
        !
        implicit none
        !
        !        Function arguments.
        !
        type(subfacenonmatchtype), intent(in) :: g1, g2
 
        ! Comparison of the iBeg value. If different set
        ! lessEqualSubfaceNonMatchType appropriately and return.
 
        if (g1%iBeg < g2%iBeg) then
            lesssubfacenonmatchtype = .true.
            return
        else if (g1%iBeg > g2%iBeg) then
            lesssubfacenonmatchtype = .false.
            return
        end if
 
        ! The iEnd value.
 
        if (g1%iEnd < g2%iEnd) then
            lesssubfacenonmatchtype = .true.
            return
        else if (g1%iEnd > g2%iEnd) then
            lesssubfacenonmatchtype = .false.
            return
        end if
 
        ! The jBeg value.
 
        if (g1%jBeg < g2%jBeg) then
            lesssubfacenonmatchtype = .true.
            return
        else if (g1%jBeg > g2%jBeg) then
            lesssubfacenonmatchtype = .false.
            return
        end if
 
        ! The jEnd value.
 
        if (g1%jEnd < g2%jEnd) then
            lesssubfacenonmatchtype = .true.
            return
        else if (g1%jEnd > g2%jEnd) then
            lesssubfacenonmatchtype = .false.
            return
        end if
 
        ! The kBeg value.
 
        if (g1%kBeg < g2%kBeg) then
            lesssubfacenonmatchtype = .true.
            return
        else if (g1%kBeg > g2%kBeg) then
            lesssubfacenonmatchtype = .false.
            return
        end if
 
        ! The kEnd value.
 
        if (g1%kEnd < g2%kEnd) then
            lesssubfacenonmatchtype = .true.
            return
        else if (g1%kEnd > g2%kEnd) then
            lesssubfacenonmatchtype = .false.
            return
        end if
 
        ! The subrange is identical. Compare the connectivity ID.
 
        if (g1%connID < g2%connID) then
            lesssubfacenonmatchtype = .true.
            return
        else if (g1%connID > g2%connID) then
            lesssubfacenonmatchtype = .false.
            return
        end if
 
        ! g1 and g2 are identical. Return .false.
 
        lesssubfacenonmatchtype = .false.
 
    end function lesssubfacenonmatchtype
 
    subroutine qsortsubblocksofcgnstype(arr, nn)
        !
        !       qsortSubblocksOfCGNSType sorts the array of the derived
        !       datatype subblocksOfCGNSType in increasing order based on the
        !       <= operator for this derived data type.
        !
        use constants
        use utils, only: terminate
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: nn
 
        type(subblocksofcgnstype), 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(subblocksofcgnstype) :: 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("qsortSubblocksOfCGNSType", &
                           "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("qsortSubblocksOfCGNSType", &
                                       "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("qsortSubblocksOfCGNSType", &
                                       "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("qsortSubblocksOfCGNSType", &
                                       "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("qsortSubblocksOfCGNSType", &
                                       "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("qsortSubblocksOfCGNSType", &
                           "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("qsortSubblocksOfCGNSType", &
                                   "Array is not sorted correctly")
            end do
        end if
 
    end subroutine qsortsubblocksofcgnstype
 
    subroutine qsortsubfacenonmatchtype(arr, nn)
        !
        !       qsortSubfaceNonMatchType sorts the array of the derived
        !       datatype subfaceNonMatchType in increasing order based on the
        !       <= operator for this derived data type.
        !
        use constants
        use utils, only: terminate
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: nn
 
        type(subfacenonmatchtype), 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(subfacenonmatchtype) :: 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("qsortSubfaceNonMatchType", &
                           "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("qsortSubfaceNonMatchType", &
                                       "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("qsortSubfaceNonMatchType", &
                                       "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("qsortSubfaceNonMatchType", &
                                       "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("qsortSubfaceNonMatchType", &
                                       "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("qsortSubfaceNonMatchType", &
                           "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("qsortSubfaceNonMatchType", &
                                   "Array is not sorted correctly")
            end do
        end if
 
    end subroutine qsortsubfacenonmatchtype
 
    logical function lessequalsortsubrangetype(g1, g2)
        !
        !         lessEqualSortSubRangeType defines the operator <= for the
        !         derived datatype sortSubRangeType. The comparison is first
        !         based on kMin, followed by jMin and finally iMin.
        !         The comparison is therefore not based on the max values.
        !
        implicit none
        !
        !        Function arguments.
        !
        type(sortsubrangetype), intent(in) :: g1, g2
        !
        !         Begin executation.
        !
        ! Compare the kMin index and return .true. or .false. if they
        ! differ.
 
        if (g1%kMin < g2%kMin) then
            lessequalsortsubrangetype = .true.
            return
        else if (g1%kMin > g2%kMin) then
            lessequalsortsubrangetype = .false.
            return
        end if
 
        ! kMin indices are equal. Compare the jMin's.
 
        if (g1%jMin < g2%jMin) then
            lessequalsortsubrangetype = .true.
            return
        else if (g1%jMin > g2%jMin) then
            lessequalsortsubrangetype = .false.
            return
        end if
 
        ! Also the jMin's are equal. Compare iMin's.
 
        if (g1%iMin < g2%iMin) then
            lessequalsortsubrangetype = .true.
            return
        else if (g1%iMin > g2%iMin) then
            lessequalsortsubrangetype = .false.
            return
        end if
 
        ! g1 equals g2. Return .true.
 
        lessequalsortsubrangetype = .true.
 
    end function lessequalsortsubrangetype
 
    !===============================================================
 
    logical function lesssortsubrangetype(g1, g2)
        !
        !         lessSortSubRangeType defines the operator < for the derived
        !         datatype sortSubRangeType. The comparison is first based on
        !         kMin, followed by jMin and finally iMin.
        !         The comparison is therefore not based on the max values.
        !
        implicit none
        !
        !        Function arguments.
        !
        type(sortsubrangetype), intent(in) :: g1, g2
        !
        !         Begin executation.
        !
        ! Compare the kMin index and return .true. or .false. if they
        ! differ.
 
        if (g1%kMin < g2%kMin) then
            lesssortsubrangetype = .true.
            return
        else if (g1%kMin > g2%kMin) then
            lesssortsubrangetype = .false.
            return
        end if
 
        ! kMin indices are equal. Compare the jMin's.
 
        if (g1%jMin < g2%jMin) then
            lesssortsubrangetype = .true.
            return
        else if (g1%jMin > g2%jMin) then
            lesssortsubrangetype = .false.
            return
        end if
 
        ! Also the jMin's are equal. Compare iMin's.
 
        if (g1%iMin < g2%iMin) then
            lesssortsubrangetype = .true.
            return
        else if (g1%iMin > g2%iMin) then
            lesssortsubrangetype = .false.
            return
        end if
 
        ! g1 equals g2. Return .false.
 
        lesssortsubrangetype = .false.
 
    end function lesssortsubrangetype
 
    !      ==================================================================
 
    subroutine sortrangessplitinfo(splitInfo)
        !
        !       sortRangesSplitInfo sort the ranges of the given subblocks in
        !       increasing order such that a unique ordering is obtained,
        !       independent of the history of the splitting.
        !
        use constants
        implicit none
        !
        !      Subroutine arguments.
        !
        type(splitcgnstype), intent(inout) :: splitInfo
        !
        !      Local variables.
        !
        integer(kind=intType) :: i, nSubBlocks
 
        type(sortsubrangetype), dimension(splitInfo%nSubBlocks) :: subRanges
 
        ! Copy the subface range from splitInfo into subRanges.
 
        nsubblocks = splitinfo%nSubBlocks
 
        do i = 1, nsubblocks
            subranges(i)%iMin = splitinfo%ranges(i, 1, 1)
            subranges(i)%jMin = splitinfo%ranges(i, 2, 1)
            subranges(i)%kMin = splitinfo%ranges(i, 3, 1)
 
            subranges(i)%iMax = splitinfo%ranges(i, 1, 2)
            subranges(i)%jMax = splitinfo%ranges(i, 2, 2)
            subranges(i)%kMax = splitinfo%ranges(i, 3, 2)
        end do
 
        ! Sort subRanges in increasing order.
 
        call qsortsortsubrangetype(subranges, nsubblocks)
 
        ! Copy the data back into splitInfo.
 
        do i = 1, nsubblocks
            splitinfo%ranges(i, 1, 1) = subranges(i)%iMin
            splitinfo%ranges(i, 2, 1) = subranges(i)%jMin
            splitinfo%ranges(i, 3, 1) = subranges(i)%kMin
 
            splitinfo%ranges(i, 1, 2) = subranges(i)%iMax
            splitinfo%ranges(i, 2, 2) = subranges(i)%jMax
            splitinfo%ranges(i, 3, 2) = subranges(i)%kMax
        end do
 
    end subroutine sortrangessplitinfo
 
    !      ==================================================================
 
    subroutine qsortsortsubrangetype(arr, nn)
        !
        !       qsortSortSubRangeType sorts the given number of halo's in
        !       increasing order based on the <= operator for this derived
        !       data type.
        !
        use utils, only: terminate
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: nn
 
        type(sortsubrangetype), 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(sortsubrangetype) :: 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("qsortSortSubRangeType", &
                           "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("qsortSortSubRangeType", &
                                       "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("qsortSortSubRangeType", &
                                       "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("qsortSortSubRangeType", &
                                       "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("qsortSortSubRangeType", &
                                       "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("qsortSortSubRangeType", &
                           "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("qsortSortSubRangeType", &
                                   "Array is not sorted correctly")
            end do
        end if
 
    end subroutine qsortsortsubrangetype
 
    !
    !         Functions to define the operators <, <= and /=.
    !         Note that the comparison is only based on the integers.
    !         The real contains additional info, the maximum deviation,
    !         which is normally different even if the subfaces are
    !         identical.
    !
    logical function lessequalfourintplusrealtype(g1, g2)
        implicit none
        type(fourintplusrealtype), intent(in) :: g1, g2
 
        ! Compare the first element.
 
        if (g1%n1 < g2%n1) then
            lessequalfourintplusrealtype = .true.
            return
        else if (g1%n1 > g2%n1) then
            lessequalfourintplusrealtype = .false.
            return
        end if
 
        ! Compare the second element.
 
        if (g1%n2 < g2%n2) then
            lessequalfourintplusrealtype = .true.
            return
        else if (g1%n2 > g2%n2) then
            lessequalfourintplusrealtype = .false.
            return
        end if
 
        ! Compare the third element.
 
        if (g1%n3 < g2%n3) then
            lessequalfourintplusrealtype = .true.
            return
        else if (g1%n3 > g2%n3) then
            lessequalfourintplusrealtype = .false.
            return
        end if
 
        ! Compare the fourth element.
 
        if (g1%n4 < g2%n4) then
            lessequalfourintplusrealtype = .true.
            return
        else if (g1%n4 > g2%n4) then
            lessequalfourintplusrealtype = .false.
            return
        end if
 
        ! g1 equals g2. Return .true.
 
        lessequalfourintplusrealtype = .true.
 
    end function lessequalfourintplusrealtype
 
    !===============================================================
 
    logical function lessfourintplusrealtype(g1, g2)
        implicit none
        type(fourintplusrealtype), intent(in) :: g1, g2
 
        ! Compare the first element.
 
        if (g1%n1 < g2%n1) then
            lessfourintplusrealtype = .true.
            return
        else if (g1%n1 > g2%n1) then
            lessfourintplusrealtype = .false.
            return
        end if
 
        ! Compare the second element.
 
        if (g1%n2 < g2%n2) then
            lessfourintplusrealtype = .true.
            return
        else if (g1%n2 > g2%n2) then
            lessfourintplusrealtype = .false.
            return
        end if
 
        ! Compare the third element.
 
        if (g1%n3 < g2%n3) then
            lessfourintplusrealtype = .true.
            return
        else if (g1%n3 > g2%n3) then
            lessfourintplusrealtype = .false.
            return
        end if
 
        ! Compare the fourth element.
 
        if (g1%n4 < g2%n4) then
            lessfourintplusrealtype = .true.
            return
        else if (g1%n4 > g2%n4) then
            lessfourintplusrealtype = .false.
            return
        end if
 
        ! g1 equals g2. Return .false.
 
        lessfourintplusrealtype = .false.
 
    end function lessfourintplusrealtype
 
    !===============================================================
 
    logical function notequalfourintplusrealtype(g1, g2)
        implicit none
        type(fourintplusrealtype), intent(in) :: g1, g2
 
        notequalfourintplusrealtype = .true.
        if (g1%n1 == g2%n1 .and. g1%n2 == g2%n2 .and. &
            g1%n3 == g2%n3 .and. g1%n4 == g2%n4) &
            notequalfourintplusrealtype = .false.
 
    end function notequalfourintplusrealtype
 
    !      ==================================================================
 
    subroutine sortbadentities(nEntities, entities, dist, sortDist)
        !
        !       sortBadEntities sorts the given number of entities in
        !       increasing order and gets rid of the multiple entries.
        !
        use constants
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(inout) :: nEntities
        integer(kind=intType), dimension(4, *), intent(inout) :: entities
        real(kind=realtype), dimension(*), intent(inout) :: dist
 
        logical, intent(in) :: sortDist
        !
        !      Local variables.
        !
        integer(kind=intType) :: nn, mm
 
        type(fourintplusrealtype), dimension(nEntities) :: tmp
 
        ! Return immediately if there are no entities to be sorted.
 
        if (nentities == 0) return
 
        ! Copy the info into tmp. If the distances must be sorted as
        ! well, copy the info. Otherwise simply put zero.
 
        do nn = 1, nentities
            tmp(nn)%n1 = entities(1, nn)
            tmp(nn)%n2 = entities(2, nn)
            tmp(nn)%n3 = entities(3, nn)
            tmp(nn)%n4 = entities(4, nn)
            if (sortdist) then
                tmp(nn)%dist = dist(nn)
            else
                tmp(nn)%dist = zero
            end if
        end do
 
        ! Sort tmp in increasing order.
 
        call qsortfourintplusrealtype(tmp, nentities)
 
        ! Get rid of the multiple entries. Note that the exceptional
        ! case of zero entities does not to be considered, because in
        ! that case this part of the subroutine is not executed.
        ! If multiple entries are present the distance is taken as
        ! the maximum of the two.
 
        mm = 1
        do nn = 2, nentities
            if (tmp(nn) /= tmp(mm)) then
                mm = mm + 1
                tmp(mm) = tmp(nn)
            else
                tmp(mm)%dist = max(tmp(mm)%dist, tmp(nn)%dist)
            end if
        end do
 
        ! Copy the data back info entities and dist. The latter
        ! only if the distances should be sorted as well.
 
        nentities = mm
 
        do nn = 1, nentities
            entities(1, nn) = tmp(nn)%n1
            entities(2, nn) = tmp(nn)%n2
            entities(3, nn) = tmp(nn)%n3
            entities(4, nn) = tmp(nn)%n4
            if (sortdist) dist(nn) = tmp(nn)%dist
        end do
 
    end subroutine sortbadentities
 
    !      ==================================================================
 
    subroutine qsortfourintplusrealtype(arr, nn)
        !
        !       qsortFourIntPlusRealType sorts the given number of halo's in
        !       increasing order based on the <= operator for this derived
        !       data type.
        !
        use constants
        use utils, only: terminate
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: nn
 
        type(fourintplusrealtype), 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(fourintplusrealtype) :: 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("qsortFourIntPlusRealType", &
                           "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("qsortFourIntPlusRealType", &
                                       "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("qsortFourIntPlusRealType", &
                                       "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("qsortFourIntPlusRealType", &
                                       "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("qsortFourIntPlusRealType", &
                                       "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("qsortFourIntPlusRealType", &
                           "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("qsortFourIntPlusRealType", &
                                   "Array is not sorted correctly")
            end do
        end if
 
    end subroutine qsortfourintplusrealtype
 
end module partitionmod