partitioning.F90

Go to the documentation of this file.

module partitioning
 
contains
 
    subroutine partitionandreadgrid(partitionOnly)
        !
        !       partitionAndReadGrid determines the partitioning of the
        !       multiblock grid over the processors and reads the grid of the
        !       blocks (or block parts) assigned to this processor. Other
        !       preprocessing activities, such as the proper setup of the halo
        !       communication structure, creation of coarse grids and wall
        !       distance computation, are performed in the preprocessing
        !       library.
        !
        use constants
        use iomodule, only: iovar
        use partitionmod, only: fileids, gridfiles
        use utils, only: terminate
        use loadbalance, only: loadbalancegrid
        use gridchecking, only: checkfaces, check1to1subfaces
        use readcgnsgrid, only: readblocksizes, readgrid
        implicit none
 
        logical, intent(in) :: partitionOnly
        !
        !      Local variables
        !
        integer :: ierr
 
        ! Determine the number of grid files that must be read,
        ! as well as the corresponding file names.
 
        call determinegridfilenames
 
        ! Read the number of blocks and the block sizes of the grid stored
        ! in the cgns grid file. This info is stored on all processors
 
        call readblocksizes
        call determineneighborids
 
        ! Determine the grid sections.
        call determinesections
 
        ! If we are just doing a partition test, return
        if (partitiononly) then
            return
        end if
 
        ! Determine the number of blocks to be stored on this processor
        ! and the relation to the original grid. Remember that blocks
        ! can be split for load balancing reasons.
 
        call loadbalancegrid
 
        ! Initialize the iblank array for the fine grid domains
 
        call initfinegridiblank
 
        ! Allocate the coordinates of the fine grid level and the
        ! derived data type used to read them.
 
        call alloccoorfinegrid
 
        ! Read the grid of the blocks (block parts) to be stored
        ! on this processor.
 
        call readgrid
 
        ! Determine for the time spectral mode the time of one period,
        ! the rotation matrices for the velocity components and
        ! create the fine grid coordinates of all time spectral locations.
 
        call timeperiodspectral
        call timerotmatricesspectral
        call finegridspectralcoor
 
        ! Release the memory of fileIDs, gridFiles and IOVar.
        ! They are not needed anymore.
 
        deallocate (fileids, gridfiles, iovar, stat=ierr)
        if (ierr /= 0) &
             call terminate("partitionAndReadGrid", &
             "Deallocation failure for fileIDs, gridFiles &
             &and IOVar")
 
        ! Check if for all faces on the block boundaries either a
        ! physical boundary condition or a connectivity has been
        ! specified and check if the 1 to 1 subfaces match.
 
        call checkfaces
        call check1to1subfaces
 
    end subroutine partitionandreadgrid
 
    subroutine determinegridfilenames
        !
        !       determineGridFileNames determines the number and names of the
        !       files that contain the grids. For steady computations only one
        !       file must be present no matter if a restart is performed or
        !       not. For unsteady the situation is a little more complicated.
        !       If no restart is performed only one file must be present. If a
        !       restart is performed in unsteady or time spectral mode and a
        !       rigid body motion is prescribed again only one grid file is
        !       required; however for a consistent restart with deforming
        !       meshes the grids in the past must be read as well. If this is
        !       not possible only a first order restart can be made in
        !       unsteady mode and some kind of interpolation is used for the
        !       time spectral method.
        !
        use constants
        use communication, only: myid, adflow_comm_world
        use inputio, only: gridfile
        use inputphysics, only: equationmode
        use inputtimespectral, only: ntimeintervalsspectral
        use inputunsteady, only: noldgridread
        use iteration, only: deforming_grid, noldlevels
        use partitionmod, only: gridfiles, interpolspectral, &
                                ngridsread, fileids
        use utils, only: terminate
        implicit none
        !
        !      Local variables
        !
        integer :: ierr
 
        integer(kind=intType) :: ii, nn, restartID
 
        character(len=7) :: integerString
        character(len=maxStringLen) :: tmpName
 
        ! Initialization of nOldGridRead and interpolSpectral.
 
        noldgridread = 1
        interpolspectral = .true.
 
        ! Determine the desired number of files from which grids at
        ! certain time levels should be read.  This depends on the
        ! equation mode we have to solve for. Also set the corresponding
        ! file names.
 
        select case (equationmode)
 
        case (steady)
 
            ! Steady computation. Only one grid needs to be read.
 
            ngridsread = 1
            allocate (fileids(ngridsread), &
                      gridfiles(ngridsread), stat=ierr)
            if (ierr /= 0) &
                 call terminate("determineGridFileNames", &
                 "Memory allocation failure for fileIDs &
                 &and gridFiles")
 
            gridfiles(1) = gridfile
 
            !===============================================================
 
        case (unsteady)
 
            ! Unsteady computation. A further check is required.
            ! EJ: replaced boolean variable restart with .false. for now
            ! Need to refactor this code as well with ALE restart
 
            testmultipleunsteady: if (deforming_grid .and. .false.) then
 
                ! A restart is made with deforming meshes. For a consistent
                ! restart nOldLevels grids must be read. First determine
                ! the prefix of the grid file and the time step number
                ! from which a restart should be made.
 
                ii = len_trim(gridfile)
                do
                    if (gridfile(ii:ii) < "0" .or. gridfile(ii:ii) > "9") exit
                    ii = ii - 1
                end do
 
                ! If the last characters of the file name do not contain a
                ! number, the grid file does not come from a previous
                ! unsteady deforming mesh computation and therefore only
                ! one grid will be read.
 
                if (ii == len_trim(gridfile)) then
 
                    ngridsread = 1
                    allocate (fileids(ngridsread), &
                              gridfiles(ngridsread), stat=ierr)
                    if (ierr /= 0) &
                         call terminate("determineGridFileNames", &
                         "Memory allocation failure for fileIDs &
                         &and gridFiles")
 
                    gridfiles(1) = gridfile
 
                else
 
                    ! Read the integer number from the last characters
                    ! of the grid file.
 
                    read (gridfile(ii + 1:), *) restartid
 
                    ! Allocate the memory for the file names and set them.
 
                    ngridsread = noldlevels
                    allocate (fileids(ngridsread), &
                              gridfiles(ngridsread), stat=ierr)
                    if (ierr /= 0) &
                         call terminate("determineGridFileNames", &
                         "Memory allocation failure for fileIDs &
                         &and gridFiles")
 
                    do nn = 1, ngridsread
                        write (integerstring, "(i6)") restartid - nn + 1
                        integerstring = adjustl(integerstring)
                        gridfiles(nn) = gridfile(:ii)//trim(integerstring)
                    end do
 
                end if
 
                else testmultipleunsteady
 
                ! The computation either starts from scratch or an unsteady
                ! restart for a rigid grid (possibly moving) is made. In
                ! all cases only one grid file is needed.
 
                ngridsread = 1
                allocate (fileids(ngridsread), &
                          gridfiles(ngridsread), stat=ierr)
                if (ierr /= 0) &
                     call terminate("determineGridFileNames", &
                     "Memory allocation failure for fileIDs &
                     &and gridFiles")
 
                gridfiles(1) = gridfile
 
            end if testmultipleunsteady
 
            ! Check if the files can be opened.
 
            do nn = 1, ngridsread
                open (unit=21, file=gridfiles(nn), status="old", iostat=ierr)
                if (ierr /= 0) exit
                close (unit=21)
            end do
 
            ! Possibly correct nGridsRead and set nOldGridRead.
            ! If nOldGridRead == 0, i.e. not a valid grid is found,
            ! print an error message and terminate.
 
            ngridsread = nn - 1
            noldgridread = ngridsread
 
            if (noldgridread == 0) then
                if (myid == 0) &
                    call terminate("determineGridFileNames", &
                                   "Grid file(s) could not be opened")
                call mpi_barrier(adflow_comm_world, ierr)
            end if
 
            !===============================================================
 
        case (timespectral)
 
            ! Time spectral computation. A further check is required.
            ! EJ: replaced boolean variable restart with .false. for now
            ! Need to refactor this code as well with ALE restart
 
            testmultiplets: if (deforming_grid .and. .false.) then
 
                ! A restart is made with deforming meshes. For a consistent
                ! restart multiple grids must be read. First determine the
                ! the prefix of the grid file from which a restart should
                ! be made.
 
                ii = len_trim(gridfile)
                do
                    if (gridfile(ii:ii) < "0" .or. gridfile(ii:ii) > "9") exit
                    ii = ii - 1
                end do
 
                ! If the last characters of the file name do not contain a
                ! number, the grid file does not come from a previous
                ! time spectral deforming mesh computation and therefore
                ! only one grid will be read.
 
                if (ii == len_trim(gridfile)) then
 
                    ngridsread = 1
                    allocate (fileids(ngridsread), &
                              gridfiles(ngridsread), stat=ierr)
                    if (ierr /= 0) &
                         call terminate("determineGridFileNames", &
                         "Memory allocation failure for fileIDs &
                         &and gridFiles")
 
                    gridfiles(1) = gridfile
 
                else
 
                    ! Loop to find out how many time instances were used in
                    ! the previous computation from which a restart is made.
 
                    nn = 0
                    do
                        nn = nn + 1
                        write (integerstring, "(i6)") nn
                        integerstring = adjustl(integerstring)
                        tmpname = gridfile(:ii)//trim(integerstring)
 
                        open (unit=21, file=tmpname, status="old", iostat=ierr)
                        if (ierr /= 0) exit
                        close (unit=21)
                    end do
 
                    nn = nn - 1
 
                    ! Take care of the exceptional situation that nn == 0.
                    ! This happens when the restart file ends at with an
                    ! integer, but does not correspond to a time spectral
                    ! solution. Allocate the memory.
 
                    ngridsread = max(nn, 1_inttype)
                    allocate (fileids(ngridsread), &
                              gridfiles(ngridsread), stat=ierr)
                    if (ierr /= 0) &
                         call terminate("determineGridFileNames", &
                         "Memory allocation failure for fileIDs &
                         &and gridFiles")
 
                    if (nn == 0) then
                        gridfiles(1) = gridfile
                    else
                        do nn = 1, ngridsread
                            write (integerstring, "(i6)") nn
                            integerstring = adjustl(integerstring)
                            gridfiles(nn) = gridfile(:ii)//trim(integerstring)
                        end do
                    end if
 
                    ! Check whether or not the coordinates must be interpolated,
                    ! i.e. check if nGridsRead == nTimeIntervalsSpectral.
 
                    if (ngridsread == ntimeintervalsspectral) &
                        interpolspectral = .false.
 
                end if
 
                else testmultiplets
 
                ! The computation either starts from scratch or a
                ! restart for a rigid grid (possibly moving) is made. In
                ! all cases only one grid file is needed.
 
                ngridsread = 1
                allocate (fileids(ngridsread), &
                          gridfiles(ngridsread), stat=ierr)
                if (ierr /= 0) &
                     call terminate("determineGridFileNames", &
                     "Memory allocation failure for fileIDs &
                     &and gridFiles")
 
                gridfiles(1) = gridfile
 
            end if testmultiplets
 
        end select
 
    end subroutine determinegridfilenames
 
    subroutine determineneighborids
        !
        !       determineNeighborIDs determines for every internal block
        !       boundary the block ID of the neighbor. In the cgns file only
        !       the zone name is stored, but the ID's are more useful
        !       internally.
        !       Although for this case a quadratic search algorithm is not too
        !       bad (number of blocks are O(1000) maximum), I don't like the
        !       idea of having a quadratic loop in the code. That's why a
        !       O(n log(n)) algorithm is used here.
        !
        use constants
        use cgnsgrid, only: cgnsdoms, cgnsndom
        use utils, only: terminate
        use sorting, only: qsortstrings, bsearchstrings
        implicit none
        !
        !      Local variables
        !
        character(len=maxCGNSNameLen), dimension(cgnsNDom) :: zoneNames
 
        integer(kind=intType), dimension(cgnsNDom) :: zoneNumbers
 
        integer(kind=intType) :: i, j, k, ii
        !
        ! Copy the zone name from the derived data type into zoneNames.
 
        do i = 1, cgnsndom
            zonenames(i) = cgnsdoms(i)%zoneName
        end do
 
        ! Sort zoneNames in increasing order.
 
        call qsortstrings(zonenames, cgnsndom)
 
        ! Initialize zoneNumbers to -1. This serves as a check during
        ! the search.
 
        zonenumbers = -1
 
        ! Find the original zone ids for the sorted zone names.
 
        do i = 1, cgnsndom
            ii = bsearchstrings(cgnsdoms(i)%zoneName, zonenames)
 
            ! Check if the zone number is not already taken. If this is the
            ! case, this means that the grid file contains two identical
            ! zone names.
 
            if (zonenumbers(ii) /= -1) &
                 call terminate("determineNeighborIDs", &
                 "Error occurs only when two identical zone &
                 &names are present")
 
            ! And set the zone number.
 
            zonenumbers(ii) = i
        end do
 
        ! Loop over the blocks and its connectivities to find out the
        ! neighbors.
 
        domains: do i = 1, cgnsndom
 
            ! The 1-to-1 connectivities.
 
            do j = 1, cgnsdoms(i)%n1to1
 
                ! Determine the neighbor ID for this internal block boundary.
 
                ii = bsearchstrings(cgnsdoms(i)%conn1to1(j)%donorName, zonenames)
                if (ii == 0) &
                    call terminate("determineNeighborIDs", &
                                   "donor name not found in sorted zone names")
 
                cgnsdoms(i)%conn1to1(j)%donorBlock = zonenumbers(ii)
 
            end do
 
            ! The non-matching abutting connectivities.
 
            do j = 1, cgnsdoms(i)%nNonMatchAbutting
 
                ! Determine the neighbor ID's for this subface.
 
                do k = 1, cgnsdoms(i)%connNonMatchAbutting(j)%nDonorBlocks
 
                    ii = bsearchstrings( &
                         cgnsdoms(i)%connNonMatchAbutting(j)%donorNames(k), zonenames)
                    if (ii == 0) &
                        call terminate("determineNeighborIDs", &
                                       "donor name not found in sorted zone names")
 
                    cgnsdoms(i)%connNonMatchAbutting(j)%donorBlocks(k) = &
                        zonenumbers(ii)
                end do
            end do
 
        end do domains
 
    end subroutine determineneighborids
 
    subroutine determineinterfaceids
        !
        !       DetermineInterfaceIDs determines more information for both the
        !       sliding mesh and domain interfaces with other codes, which are
        !       both specified as user defined boundary conditions in CGNS. In
        !       particular the number of sliding mesh interfaces and their
        !       pairings, and number of interfaces with other codes are
        !       determined.
        !       There are some parts in the coupler API routines where
        !       the family-specified domain interfaces are implicitly assumed.
        !       Therefore, it is recommended for the time being that all the
        !       domain interfaces should be family-specified.
        !
        use constants
        use cgnsgrid, only: cgnsdoms, cgnsndom, cgnsfamilies, &
                            famidsdomaininterfaces, bcidsdomaininterfaces, famidssliding, &
                            cgnsnfamilies, cgnsnsliding, cgnsndomaininterfaces
        use communication, only: adflow_comm_world, myid
        use iteration, only: standalonemode
        use utils, only: terminate
        use sorting, only: qsortstrings, bsearchstrings
        use commonformats, only: stringspace
 
        implicit none
        !
        !      Local variables.
        !
        integer :: ierr
 
        integer(kind=intType) :: ii, jj, mm, nn
        integer(kind=intType) :: nSlidingFam, nSlidingBC, nSliding
        integer(kind=intType) :: nDomainFam, nDomainBC
 
        integer(kind=intType), dimension(:), allocatable :: famSlidingID
        integer(kind=intType), dimension(:, :), allocatable :: bcSlidingID
        integer(kind=intType), dimension(:), allocatable :: orID
 
        character(len=maxStringLen) :: errorMessage
 
        character(len=maxCGNSNameLen), dimension(:), allocatable :: &
            namesSliding
        character(len=maxCGNSNameLen), dimension(:), allocatable :: &
            namesSorted
 
        logical :: validInterface
        !
        ! Count the total number of each type of user-defined BC that
        ! needs to be treated here. Note that if a BC is specified by the
        ! the family it is not counted, such that BCs making up a user-
        ! defined family are only counted once.
 
        nslidingfam = 0
        nslidingbc = 0
        ndomainfam = 0
        ndomainbc = 0
 
        do nn = 1, cgnsnfamilies
 
            select case (cgnsfamilies(nn)%BCType)
            case (slidinginterface)
                nslidingfam = nslidingfam + 1
            case (domaininterfaceall, domaininterfacerhouvw, &
                  domaininterfacep, domaininterfacerho, &
                  domaininterfacetotal)
                ndomainfam = ndomainfam + 1
            end select
 
        end do
 
        do nn = 1, cgnsndom
            do mm = 1, cgnsdoms(nn)%nBocos
                if (cgnsdoms(nn)%bocoInfo(mm)%actualFace .and. &
                    cgnsdoms(nn)%bocoInfo(mm)%familyID == 0) then
 
                    select case (cgnsdoms(nn)%bocoInfo(mm)%BCType)
                    case (slidinginterface)
                        nslidingbc = nslidingbc + 1
                    case (domaininterfaceall, domaininterfacerhouvw, &
                          domaininterfacep, domaininterfacerho, &
                          domaininterfacetotal)
                        ndomainbc = ndomainbc + 1
                    end select
 
                end if
            end do
        end do
 
        ! Domain interfaces are only allowed when the code is run in a
        ! multi-disciplinary environment. Check this.
 
        cgnsndomaininterfaces = ndomainfam + ndomainbc
 
        if (standalonemode .and. cgnsndomaininterfaces > 0) then
            if (myid == 0) &
                 call terminate("determineInterfaceIDs", &
                 "Domain interfaces are not allowed in &
                 &stand alone mode.")
            call mpi_barrier(adflow_comm_world, ierr)
        end if
 
        ! The number of sliding mesh interfaces must be even, because each
        ! sliding interface should have two sides. Check this.
 
        nsliding = nslidingfam + nslidingbc
        if (mod(nsliding, 2) == 1) then
            if (myid == 0) &
                call terminate("determineInterfaceIDs", &
                               "Odd number of sliding mesh families found")
            call mpi_barrier(adflow_comm_world, ierr)
        end if
 
        ! Allocate memory to store the names and IDs for the interfaces.
 
        allocate (famidsdomaininterfaces(ndomainfam), &
                  bcidsdomaininterfaces(2, ndomainbc), &
                  namessliding(nsliding), namessorted(nsliding), &
                  orid(nsliding), famslidingid(nslidingfam), &
                  bcslidingid(2, nslidingbc), stat=ierr)
        if (ierr /= 0) &
            call terminate("determineInterfaceIDs", &
                           "Memory allocation failure for names, IDs")
 
        ! Loop back over the families again and store the names and
        ! IDs this time around. Note the ID is just the index of the
        ! corresponding family.
 
        ii = 0
        jj = 0
 
        do nn = 1, cgnsnfamilies
 
            select case (cgnsfamilies(nn)%BCType)
            case (slidinginterface)
                ii = ii + 1
                namessliding(ii) = cgnsfamilies(nn)%familyName
                famslidingid(ii) = nn
            case (domaininterfaceall, domaininterfacerhouvw, &
                  domaininterfacep, domaininterfacerho, &
                  domaininterfacetotal)
                jj = jj + 1
                famidsdomaininterfaces(jj) = nn
            end select
 
        end do
 
        ! Loop back over the boundary conditions again and store the
        ! names and IDs this time around. Note the ID has two parts:
        ! the domain and the index of the BC info in that domain.
 
        jj = 0
 
        do nn = 1, cgnsndom
            do mm = 1, cgnsdoms(nn)%nBocos
                if (cgnsdoms(nn)%bocoInfo(mm)%actualFace .and. &
                    cgnsdoms(nn)%bocoInfo(mm)%familyID == 0) then
 
                    select case (cgnsdoms(nn)%bocoInfo(mm)%BCType)
                    case (slidinginterface)
                        ii = ii + 1
                        namessliding(ii) = cgnsdoms(nn)%bocoInfo(mm)%bocoName
                        bcslidingid(:, ii - nslidingfam) = (/nn, mm/)
                    case (domaininterfaceall, domaininterfacerhouvw, &
                          domaininterfacep, domaininterfacerho, &
                          domaininterfacetotal)
                        jj = jj + 1
                        bcidsdomaininterfaces(:, jj) = (/nn, mm/)
                    end select
 
                end if
            end do
        end do
 
        ! Initialize orID to -1, which serves as a check later on, and
        ! copy the names of the sliding mesh families in namesSorted.
 
        do ii = 1, nsliding
            orid(ii) = -1
            namessorted(ii) = namessliding(ii)
        end do
 
        ! Sort the names of the sliding mesh families in increasing
        ! order and find the corresponding entry in namesSliding.
 
        call qsortstrings(namessorted, nsliding)
 
        do ii = 1, nsliding
 
            ! Search the sorted strings.
 
            mm = bsearchstrings(namessliding(ii), namessorted)
 
            if (orid(mm) /= -1) then
 
                ! Family name occurs more than once. This is not allowed.
 
                write (errormessage, stringspace) "Family name", trim(namessliding(ii)), "occurs more than once."
                if (myid == 0) &
                    call terminate("determineInterfaceIDs", errormessage)
                call mpi_barrier(adflow_comm_world, ierr)
 
            end if
 
            ! Store the entry in orID.
 
            orid(mm) = ii
 
        end do
 
        ! Set the number of sliding mesh interfaces and allocate the
        ! memory for famIDsSliding.
 
        cgnsnsliding = nsliding / 2
        allocate (famidssliding(cgnsnsliding, 2), stat=ierr)
        if (ierr /= 0) &
            call terminate("determineInterfaceIDs", &
                           "Memory allocation failure for famIDsSliding")
 
        ! Check if the sorted family names indeed form a
        ! sliding mesh interface.
 
        do ii = 1, cgnsnsliding
 
            ! Store the entries in namesSorted, which should form a sliding
            ! interface, a bit easier.
 
            nn = 2 * ii
            mm = nn - 1
 
            ! Store the original values in famIDsSliding.
 
            jj = orid(nn) / 2
 
            famidssliding(jj, 1) = famslidingid(orid(mm))
            famidssliding(jj, 2) = famslidingid(orid(nn))
 
            ! Check if the names form a valid sliding interface.
 
            validinterface = .true.
            if (len_trim(namessorted(nn)) /= len_trim(namessorted(mm))) &
                validinterface = .false.
 
            jj = len_trim(namessorted(nn)) - 2
            if (jj <= 0) then
                validinterface = .false.
            else if (namessorted(nn) (:jj) /= namessorted(mm) (:jj)) then
                validinterface = .false.
            end if
 
            ! Print an error message and exit if the two families do not
            ! form a valid sliding interface.
 
            if (.not. validinterface) then
                write (errormessage, stringspace) "Family names", trim(namessliding(nn)), "and", &
                    trim(namessliding(mm)), "do not form a valid sliding mesh interface"
                if (myid == 0) &
                    call terminate("determineInterfaceIDs", errormessage)
                call mpi_barrier(adflow_comm_world, ierr)
            end if
 
            ! The two names form a valid sliding interface. Store the
            ! sliding interface ID for the original family or BC. Note that
            ! one gets a positive ID and the other one a negative ID.
            ! In this way a distinction is made between the two sides.
 
            if (orid(mm) > nslidingfam) then
                jj = bcslidingid(2, orid(mm))
                mm = bcslidingid(1, orid(mm))
                cgnsdoms(mm)%bocoInfo(jj)%slidingID = -ii
            else
                mm = famslidingid(orid(mm))
                cgnsfamilies(mm)%slidingID = -ii
            end if
 
            if (orid(nn) > nslidingfam) then
                jj = bcslidingid(2, orid(nn))
                nn = bcslidingid(1, orid(nn))
                cgnsdoms(nn)%bocoInfo(jj)%slidingID = ii
            else
                nn = famslidingid(orid(nn))
                cgnsfamilies(nn)%slidingID = ii
            end if
 
        end do
 
        ! Loop over all the boundary conditions again, this time
        ! copying the sliding interface IDs for those BCs that
        ! were part of a sliding family.
 
        do nn = 1, cgnsndom
            do mm = 1, cgnsdoms(nn)%nBocos
                if (cgnsdoms(nn)%bocoInfo(mm)%actualFace .and. &
                    cgnsdoms(nn)%bocoInfo(mm)%familyID > 0 .and. &
                    cgnsdoms(nn)%bocoInfo(mm)%BCType == slidinginterface) then
 
                    cgnsdoms(nn)%bocoInfo(mm)%slidingID = &
                        cgnsfamilies(cgnsdoms(nn)%bocoInfo(mm)%familyID)%slidingID
 
                end if
            end do
        end do
 
        ! Deallocate the memory used to determine the sliding IDs.
 
        deallocate (namessliding, namessorted, orid, &
                    famslidingid, bcslidingid, stat=ierr)
        if (ierr /= 0) &
            call terminate("determineInterfaceIDs", &
                           "Deallocation failure for names, IDs")
 
    end subroutine determineinterfaceids
 
    subroutine initfinegridiblank
        !
        !       InitFineGridIblank allocates the fine grid iblank array and
        !       initializes the values for the holes, boundary, and halos. The
        !       holes read into the cgns domains are distributed amongst its
        !       sublocks in the form of iblanks. That is, we do not store a
        !       list of indices for the holes of the flow domains as done in
        !       the CGNS. The number of holes in each domain are also counted.
        !
        use constants
        use block, only: ndom, flowdoms
        use utils, only: terminate
        use inputtimespectral, only: ntimeintervalsspectral
        implicit none
        !
        !      Local variables.
        !
        integer :: ierr
 
        integer(kind=intType) :: i, j, k, l, m, n, cgnsId, sps
 
        ! Loop over the local blocks.
        spectralloop: do sps = 1, ntimeintervalsspectral
            domains: do n = 1, ndom
 
                ! Allocate memory for the iblank array of this block.
 
                i = flowdoms(n, 1, sps)%ib
                j = flowdoms(n, 1, sps)%jb
                k = flowdoms(n, 1, sps)%kb
                allocate (flowdoms(n, 1, sps)%iblank(0:i, 0:j, 0:k), &
                          flowdoms(n, 1, sps)%forcedRecv(0:i, 0:j, 0:k), &
                          flowdoms(n, 1, sps)%status(0:i, 0:j, 0:k), &
                          stat=ierr)
                if (ierr /= 0) &
                    call terminate("initFineGridIblank", &
                                   "Memory allocation failure for iblank")
 
                ! Initialize iblank to 1 everywhere, and the number of holes
                ! for this domain to 0.
 
                flowdoms(n, 1, sps)%iblank = 1
                flowdoms(n, 1, sps)%forcedRecv = 0
                flowdoms(n, 1, sps)%status = 0
            end do domains
        end do spectralloop
 
    end subroutine initfinegridiblank
 
    subroutine timeperiodspectral
        !
        !       timePeriodSpectral determines the time of one period for the
        !       time spectral method. It is possible that sections have
        !       different periodic times.
        !
        use constants
        use communication, only: myid, adflow_comm_world
        use inputmotion, only: degreefourxrot, degreefouryrot, degreefourzrot, &
                               omegafouralpha, omegafourbeta, omegafourmach, omegafourxrot, &
                               omegafouryrot, omegafourzrot, degreefourmach, degreefouralpha, &
                               degreefourbeta
        use inputphysics, only: equationmode, flowtype
        use inputtimespectral, only: omegafourier
        use section, only: sections, nsections
        use utils, only: terminate
        implicit none
        !
        !      Local parameter.
        !
        real(kind=realtype), parameter :: tol = 1.e-5_realtype
        !
        !      Local variables.
        !
        integer :: ierr
        integer(kind=intType) :: nn
 
        real(kind=realtype) :: tt, omega
        real(kind=realtype) :: timeperiod
 
        logical :: timeDetermined
 
        ! This routine is only used for the spectral solutions. Return
        ! immediately if a different mode is solved.
 
        if (equationmode /= timespectral) return
 
        ! First check if a rotational frequency has been specified.
        ! Only for external flows.
 
        timedetermined = .false.
 
        externaltest: if (flowtype == externalflow) then
 
            ! X-rotation.
 
            if (degreefourxrot > 0) then
                timeperiod = two * pi / omegafourxrot
                timedetermined = .true.
            end if
 
            ! Y-rotation.
 
            if (degreefouryrot > 0) then
                tt = two * pi / omegafouryrot
 
                ! Check if a time period was already determined. If so, try
                ! to determine a common time. Otherwise just copy the data.
 
                if (timedetermined) then
                    timeperiod = commontimespectral(timeperiod, tt)
                else
                    timeperiod = tt
                    timedetermined = .true.
                end if
            end if
 
            ! Z-rotation.
 
            if (degreefourzrot > 0) then
                tt = two * pi / omegafourzrot
 
                ! Check if a time period was already determined. If so, try
                ! to determine a common time. Otherwise just copy the data.
 
                if (timedetermined) then
                    timeperiod = commontimespectral(timeperiod, tt)
                else
                    timeperiod = tt
                    timedetermined = .true.
                end if
            end if
 
            ! Alpha
            !print *,'degreeFourAlpha',degreefouralpha,omegafouralpha,sincoeffouralpha
            if (degreefouralpha > 0) then
                tt = two * pi / omegafouralpha
                !print *,'timePeriod',tt
                ! Check if a time period was already determined. If so, try
                ! to determine a common time. Otherwise just copy the data.
 
                if (timedetermined) then
                    timeperiod = commontimespectral(timeperiod, tt)
                else
                    timeperiod = tt
                    timedetermined = .true.
                end if
            end if
 
            ! Beta
 
            if (degreefourbeta > 0) then
                tt = two * pi / omegafourbeta
 
                ! Check if a time period was already determined. If so, try
                ! to determine a common time. Otherwise just copy the data.
 
                if (timedetermined) then
                    timeperiod = commontimespectral(timeperiod, tt)
                else
                    timeperiod = tt
                    timedetermined = .true.
                end if
            end if
 
            ! Mach
 
            if (degreefourmach > 0) then
                tt = two * pi / omegafourmach
 
                ! Check if a time period was already determined. If so, try
                ! to determine a common time. Otherwise just copy the data.
 
                if (timedetermined) then
                    timeperiod = commontimespectral(timeperiod, tt)
                else
                    timeperiod = tt
                    timedetermined = .true.
                end if
            end if
 
            ! aeroelastic case
            if (omegafourier > 0) then
                tt = two * pi / omegafourier
 
                ! Check if a time period was already determined. If so, try
                ! to determine a common time. Otherwise just copy the data.
 
                if (timedetermined) then
                    timeperiod = commontimespectral(timeperiod, tt)
                else
                    timeperiod = tt
                    timedetermined = .true.
                end if
            end if
 
!!$         ! Altitude.
!!$
!!$         if(degreeFourAltitude > 0) then
!!$           tt = two*pi/omegaFourAltitude
!!$
!!$           ! Check if a time period was already determined. If so, try
!!$           ! to determine a common time. Otherwise just copy the data.
!!$
!!$           if( timeDetermined ) then
!!$             timePeriod = commonTimeSpectral(timePeriod, tt)
!!$           else
!!$             timePeriod     = tt
!!$             timeDetermined = .true.
!!$           endif
!!$         endif
 
        end if externaltest
 
        ! If it was possible to determine the time, copy it to the
        ! sections and return.
 
        if (timedetermined) then
            do nn = 1, nsections
                sections(nn)%timePeriod = timeperiod / sections(nn)%nSlices
                !print *,'sectionTimePeriod',sections(nn)%timePeriod,nn
            end do
            return
        end if
 
        ! Try to determine the periodic time via the rotation rate of the
        ! sections and its number of slices.
 
        sectionloop: do nn = 1, nsections
 
            ! Test if the section is rotating, because only for rotating
            ! sections the periodic time can be determined.
 
            testrotating: if (sections(nn)%rotating) then
 
                ! Determine the magnitude of the rotation rate and the
                ! corresponding periodic time period.
 
                omega = sqrt(sections(nn)%rotRate(1)**2 &
                             + sections(nn)%rotRate(2)**2 &
                             + sections(nn)%rotRate(3)**2)
 
                tt = two * pi / omega
 
                ! If a time period was already determined, check if this is
                ! identical to tt. If not print an error message and exit.
 
                if (timedetermined) then
 
                    tt = abs(tt - timeperiod) / timeperiod
                    if (tt > tol) then
                        if (myid == 0) &
                             call terminate("timePeriodSpectral", &
                             "Rotational frequencies of the rotating &
                             &sections are not identical.")
                        call mpi_barrier(adflow_comm_world, ierr)
                    end if
 
                else
 
                    ! Just copy the data.
 
                    timeperiod = tt
                    timedetermined = .true.
 
                end if
 
            end if testrotating
        end do sectionloop
 
        ! Divide the periodic time by the number of slices to get the
        ! characteristic time for every section.
 
        do nn = 1, nsections
            sections(nn)%timePeriod = timeperiod / sections(nn)%nSlices
        end do
 
        ! Return if it was possible to determine the time.
 
        if (timedetermined) return
 
        ! Periodic time could not be determined. Print an error
        ! message and exit.
 
        if (myid == 0) &
             call terminate("timePeriodSpectral", &
             "Not possible to determine the periodic time &
             &for the time spectral method")
        call mpi_barrier(adflow_comm_world, ierr)
 
    end subroutine timeperiodspectral
 
    function commontimespectral(t1, t2)
        !
        !       The function commonTimeSpectral determines the smallest
        !       possible common time between t1 and t2, such that
        !       tcommon = n1*t1 = n2*t2 and n1, n2 integers.
        !
        use communication
        use utils, only: terminate
        implicit none
        !
        !      Function definition
        !
        real(kind=realtype) :: commontimespectral
        !
        !      Function arguments.
        !
        real(kind=realtype), intent(in) :: t1, t2
        !
        !      Local parameters.
        !
        integer(kind=intType), parameter :: nmax = 100
        real(kind=realtype), parameter :: tol = 1.e-5_realtype
        !
        !      Local variables.
        !
        integer :: ierr
        integer(kind=intType) :: n1, n2
        real(kind=realtype) :: tt1, tt2, tt, ratio
 
        ! Store the largest time in tt1 and the smallest in tt2 and
        ! compute the ratio tt1/tt2, which is >= 1
 
        tt1 = max(t1, t2)
        tt2 = min(t1, t2)
        ratio = tt1 / tt2
 
        ! Loop to find the smallest integer values of n1 and n2, such
        ! that n1*tt1 = n2*tt2. Note that due to the previous definition
        ! n2 >= n1.
 
        do n1 = 1, nmax
            tt = n1 * ratio
            n2 = nint(tt)
            tt = abs(tt - n2)
            if (tt <= tol) exit
        end do
 
        ! Check if a common time was found
 
        if (n1 > nmax) then
            if (myid == 0) &
                 call terminate("commonTimeSpectral", &
                 "No common time periodic time found. &
                 &Problem may not be periodic")
            call mpi_barrier(adflow_comm_world, ierr)
        end if
 
        ! Set the common time.
 
        commontimespectral = n1 * tt1
 
    end function commontimespectral
    subroutine timerotmatricesspectral
        !
        !       timeRotMatricesSpectral determines the rotation matrices
        !       used in the time derivatives for the velocity components in
        !       the time spectral method. These matrices are the identity
        !       matrices for non-rotating sections and something different for
        !       rotating sections. Therefore the rotation matrices are stored
        !       for every section.
        !
        use constants
        use inputphysics, only: equationmode
        use inputtimespectral, only: rotmatrixspectral
        use section, only: sections, nsections
        use utils, only: terminate
        implicit none
        !
        !      Local variables.
        !
        integer :: ierr
 
        integer(kind=intType) :: nn
        real(kind=realtype) :: tmp, theta, costheta, sintheta
 
        real(kind=realtype), dimension(3) :: xt, yt, zt
 
        ! This routine is only used for the spectral solutions. Return
        ! immediately if a different mode is solved.
 
        if (equationmode /= timespectral) return
 
        ! Allocate the memory for rotMatrixSpectral, which will store
        ! the rotation matrices for all the sections.
 
        if (allocated(rotmatrixspectral)) deallocate (rotmatrixspectral)
        allocate (rotmatrixspectral(nsections, 3, 3), stat=ierr)
 
        if (ierr /= 0) &
             call terminate("timeRotMatricesSpectral", &
             "Memory allocation failure for &
             &rotMatrixSpectral")
 
        ! Loop over the number of sections.
 
        sectionloop: do nn = 1, nsections
 
            ! Test if the rotation matrix is the unity matrix. This is the
            ! case if this section is not rotating or if the number of
            ! slices is only 1, i.e. the true physical model is computed.
 
            testunity: if (.not. sections(nn)%rotating .or. &
                           sections(nn)%nSlices == 1) then
 
                ! Set the rotation matrix to the unity matrix.
 
                rotmatrixspectral(nn, 1, 1) = one
                rotmatrixspectral(nn, 1, 2) = zero
                rotmatrixspectral(nn, 1, 3) = zero
 
                rotmatrixspectral(nn, 2, 1) = zero
                rotmatrixspectral(nn, 2, 2) = one
                rotmatrixspectral(nn, 2, 3) = zero
 
                rotmatrixspectral(nn, 3, 1) = zero
                rotmatrixspectral(nn, 3, 2) = zero
                rotmatrixspectral(nn, 3, 3) = one
 
                else testunity
 
                ! Section is rotating and only a part of the physical problem
                ! is modelled. Consequently a rotation matrix is present for
                ! the velocity components.
 
                ! First transform to a frame where the xt-axis points in the
                ! direction of the rotation vector.
 
                xt(1) = sections(nn)%rotAxis(1)
                xt(2) = sections(nn)%rotAxis(2)
                xt(3) = sections(nn)%rotAxis(3)
 
                ! Construct the yt axis. It does not matter exactly as long
                ! as it is normal to xt.
 
                if (abs(xt(2)) < 0.707107_realtype) then
                    yt(1) = zero
                    yt(2) = one
                    yt(3) = zero
                else
                    yt(1) = zero
                    yt(2) = zero
                    yt(3) = one
                end if
 
                ! Make sure that yt is normal to xt.
 
                tmp = xt(1) * yt(1) + xt(2) * yt(2) + xt(3) * yt(3)
                yt(1) = yt(1) - tmp * xt(1)
                yt(2) = yt(2) - tmp * xt(2)
                yt(3) = yt(3) - tmp * xt(3)
 
                ! And create a unit vector.
 
                tmp = one / sqrt(yt(1)**2 + yt(2)**2 + yt(3)**2)
                yt(1) = tmp * yt(1)
                yt(2) = tmp * yt(2)
                yt(3) = tmp * yt(3)
 
                ! Create the vector zt by taking the cross product xt*yt.
 
                zt(1) = xt(2) * yt(3) - xt(3) * yt(2)
                zt(2) = xt(3) * yt(1) - xt(1) * yt(3)
                zt(3) = xt(1) * yt(2) - xt(2) * yt(1)
 
                ! Compute the periodic angle theta and its sine and cosine.
 
                theta = two * pi / real(sections(nn)%nSlices, realtype)
                costheta = cos(theta)
                sintheta = sin(theta)
 
                ! The rotation matrix in the xt,yt,zt frame is given by
                !
                ! R = | 1      0           0      |
                !     | 0  cos(theta) -sin(theta) |
                !     | 0  sin(theta)  cos(theta) |
                !
                ! The rotation matrix in the standard cartesian frame is then
                ! given by t * r * t^t, where the colums of the transformation
                ! matrix t are the unit vectors xt,yt,zt. One can easily check
                ! this by checking rotation around the y- and z-axis. The
                ! result of this is the expression below.
 
                rotmatrixspectral(nn, 1, 1) = xt(1) * xt(1) &
                                              + costheta * (yt(1) * yt(1) + zt(1) * zt(1))
                rotmatrixspectral(nn, 1, 2) = xt(1) * xt(2) &
                                              + costheta * (yt(1) * yt(2) + zt(1) * zt(2)) &
                                              - sintheta * (yt(1) * zt(2) - yt(2) * zt(1))
                rotmatrixspectral(nn, 1, 3) = xt(1) * xt(3) &
                                              + costheta * (yt(1) * yt(3) + zt(1) * zt(3)) &
                                              - sintheta * (yt(1) * zt(3) - yt(3) * zt(1))
 
                rotmatrixspectral(nn, 2, 1) = xt(1) * xt(2) &
                                              + costheta * (yt(1) * yt(2) + zt(1) * zt(2)) &
                                              + sintheta * (yt(1) * zt(2) - yt(2) * zt(1))
                rotmatrixspectral(nn, 2, 2) = xt(2) * xt(2) &
                                              + costheta * (yt(2) * yt(2) + zt(2) * zt(2))
                rotmatrixspectral(nn, 2, 3) = xt(2) * xt(3) &
                                              + costheta * (yt(2) * yt(3) + zt(2) * zt(3)) &
                                              - sintheta * (yt(2) * zt(3) - yt(3) * zt(2))
 
                rotmatrixspectral(nn, 3, 1) = xt(1) * xt(3) &
                                              + costheta * (yt(1) * yt(3) + zt(1) * zt(3)) &
                                              + sintheta * (yt(1) * zt(3) - yt(3) * zt(1))
                rotmatrixspectral(nn, 3, 2) = xt(2) * xt(3) &
                                              + costheta * (yt(2) * yt(3) + zt(2) * zt(3)) &
                                              + sintheta * (yt(2) * zt(3) - yt(3) * zt(2))
                rotmatrixspectral(nn, 3, 3) = xt(3) * xt(3) &
                                              + costheta * (yt(3) * yt(3) + zt(3) * zt(3))
 
            end if testunity
 
        end do sectionloop
 
    end subroutine timerotmatricesspectral
 
    subroutine finegridspectralcoor
        !
        !       fineGridSpectralCoor computes the coordinates of all but
        !       the first spectral solution from the known coordinates of the
        !       first time instance.
        !
        use constants
        use block, only: flowdoms, ndom
        use inputphysics, only: equationmode
        use inputtimespectral, only: ntimeintervalsspectral
        use iomodule, only: iovar
        use iteration, only: currentlevel
        use monitor, only: timeunsteady
        use section, only: nsections, sections
        use partitionmod, only: interpolspectral, ngridsread
        use utils, only: terminate
        implicit none
        !
        !      Local variables.
        !
        integer :: ierr
 
        integer(kind=intType) :: nn, ll, i, j, k
        integer(kind=intType) :: il, jl, kl
 
        real(kind=realtype), dimension(nSections) :: dt, t
 
        ! This routine is only used for the spectral solutions. Return
        ! immediately if a different mode is solved.
 
        if (equationmode /= timespectral) return
 
        ! Also return immediately if all coordinates were already read
        ! from the grid file.
 
        if (.not. interpolspectral) return
        !
        !       Step 1. Perform a rigid body motion of the coordinates of the
        !               1st time instance to the other instances.
        !
        ! Set currentLevel to 1, such that updateCoorFineMesh
        ! updates the coordinates of the correct level.
 
        currentlevel = 1
 
        ! Determine the delta t for every section. Remember it is possible
        ! that every section has a different periodic time.
 
        do nn = 1, nsections
            dt(nn) = sections(nn)%timePeriod &
                     / real(ntimeintervalsspectral, realtype)
        end do
 
        ! Loop over the number of spectral solutions, starting at 2,
        ! because the first is already known.
 
        timeunsteady = zero
        spectralloop: do ll = 2, ntimeintervalsspectral
 
            ! Set the owned coordinates to the coordinates of the first
            ! solution such that updateCoorFineMesh can modify them.
 
            do nn = 1, ndom
 
                do k = 1, flowdoms(nn, 1, 1)%kl
                    do j = 1, flowdoms(nn, 1, 1)%jl
                        do i = 1, flowdoms(nn, 1, 1)%il
                            flowdoms(nn, 1, ll)%x(i, j, k, 1) = flowdoms(nn, 1, 1)%x(i, j, k, 1)
                            flowdoms(nn, 1, ll)%x(i, j, k, 2) = flowdoms(nn, 1, 1)%x(i, j, k, 2)
                            flowdoms(nn, 1, ll)%x(i, j, k, 3) = flowdoms(nn, 1, 1)%x(i, j, k, 3)
                        end do
                    end do
                end do
 
            end do
 
            ! Compute the corresponding times for this spectral solution
            ! and call updateCoorFineMesh to determine the coordinates.
 
            do nn = 1, nsections
                t(nn) = (ll - 1) * dt(nn)
            end do
 
            call updatecoorfinemesh(t, ll)
 
        end do spectralloop
 
        ! Return if only one grid has been read.
 
        if (ngridsread == 1) return
        !
        !       Step 2. Multiple grids have been read, but the number is not
        !               equal to the number of time instances used in the
        !               computation. As multiple grids have been read this
        !               means that a time spectral computation on a deforming
        !               mesh is performed. Therefore the deformations,
        !               relative to the first grid, must be interpolated.
        !
        ! First allocate the memory of IOVar(..,1)%w.
        ! This will serve as temporary storage for the coordinates of
        ! the 1st spectral instance, because those will be used in the
        ! call to updateCoorFineMesh. In this way this routine can be
        ! used without modification.
 
        do nn = 1, ndom
 
            kl = flowdoms(nn, 1, 1)%kl
            jl = flowdoms(nn, 1, 1)%jl
            il = flowdoms(nn, 1, 1)%il
 
            allocate (iovar(nn, 1)%w(il, jl, kl, 3), stat=ierr)
            if (ierr /= 0) &
                 call terminate("fineGridSpectralCoor", &
                 "Memory allocation failure for &
                 &IOVar(nn,1)%w")
 
            do k = 1, kl
                do j = 1, jl
                    do i = 1, il
                        iovar(nn, 1)%w(i, j, k, 1) = flowdoms(nn, 1, 1)%x(i, j, k, 1)
                        iovar(nn, 1)%w(i, j, k, 2) = flowdoms(nn, 1, 1)%x(i, j, k, 2)
                        iovar(nn, 1)%w(i, j, k, 3) = flowdoms(nn, 1, 1)%x(i, j, k, 3)
                    end do
                end do
            end do
 
        end do
 
        ! Determine the delta t for every section. Remember it is possible
        ! that every section has a different periodic time.
 
        do nn = 1, nsections
            dt(nn) = sections(nn)%timePeriod &
                     / real(ngridsread, realtype)
        end do
 
        ! Loop over the number of spectral solutions read, starting at
        ! 2, and determine the displacements relative to the rigid body
        ! motion of the grid of the 1st time instance.
 
        timeunsteady = zero
        spectralloopread: do ll = 2, ngridsread
 
            ! Compute the corresponding times for this spectral solution
            ! and call updateCoorFineMesh to determine the coordinates.
 
            do nn = 1, nsections
                t(nn) = (ll - 1) * dt(nn)
            end do
 
            call updatecoorfinemesh(t, 1_inttype)
 
            ! Determine the relative displacements for this time instance
            ! and initialize flowDoms(nn,1,1) for the next round.
 
            do nn = 1, ndom
 
                do k = 1, flowdoms(nn, 1, 1)%kl
                    do j = 1, flowdoms(nn, 1, 1)%jl
                        do i = 1, flowdoms(nn, 1, 1)%il
                            iovar(nn, ll)%w(i, j, k, 1) = iovar(nn, ll)%w(i, j, k, 1) &
                                                          - flowdoms(nn, 1, 1)%x(i, j, k, 1)
                            iovar(nn, ll)%w(i, j, k, 2) = iovar(nn, ll)%w(i, j, k, 2) &
                                                          - flowdoms(nn, 1, 1)%x(i, j, k, 2)
                            iovar(nn, ll)%w(i, j, k, 3) = iovar(nn, ll)%w(i, j, k, 3) &
                                                          - flowdoms(nn, 1, 1)%x(i, j, k, 3)
 
                            flowdoms(nn, 1, 1)%x(i, j, k, 1) = iovar(nn, 1)%w(i, j, k, 1)
                            flowdoms(nn, 1, 1)%x(i, j, k, 2) = iovar(nn, 1)%w(i, j, k, 2)
                            flowdoms(nn, 1, 1)%x(i, j, k, 3) = iovar(nn, 1)%w(i, j, k, 3)
                        end do
                    end do
                end do
 
            end do
 
        end do spectralloopread
 
        ! The coordinates of IOVar now contain the relative
        ! displacements compared to the rigid body motion of the first
        ! time instance, except for the first time instance.
        ! Set these to zero.
 
        do nn = 1, ndom
            iovar(nn, 1)%w = zero
        end do
 
        ! Interpolate the displacements and add them to the currently
        ! stored coordinates.
 
        call terminate("fineGridSpectralCoor", &
                       "Arti should do this interpolation stuff")
 
        ! Release the memory of the variable w in IOVar.
 
        do nn = 1, ndom
            do ll = 1, ngridsread
                deallocate (iovar(nn, ll)%w, stat=ierr)
                if (ierr /= 0) &
                    call terminate("fineGridSpectralCoor", &
                                   "Deallocation failure for IOVar%w")
            end do
        end do
 
    end subroutine finegridspectralcoor
    subroutine updatecoorfinemesh(dtAdvance, sps)
        !
        !       updateCoorFineMesh updates the coordinates of the
        !       moving parts of the current finest mesh by the given amount of
        !       time, possibly different per section. In unsteady mode all the
        !       times will be equal, but in time spectral mode they can be
        !       different.
        !       This routine is called in the full mg cycle to put the fine
        !       mesh to the position previously calculated on the coarser
        !       grid levels, in the unsteady time loop to advance the
        !       coordinates only one time step and in the partitioning part
        !       of the spectral mode to compute the coordinates of the given
        !       spectral solution sps. As it is used in the full MG cycle,
        !       currentLevel points to the correct grid level and not
        !       ground level.
        !
        use constants
        use block
        use blockpointers
        use flowvarrefstate
        use cgnsgrid
        use inputmotion
        use iteration
        use monitor
        use utils, only: setpointers, rotmatrixrigidbody
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: sps
 
        real(kind=realtype), dimension(*), intent(in) :: dtadvance
        !
        !      Local variables.
        !
        integer(kind=intType) :: i, j, k, nn
 
        real(kind=realtype) :: displx, disply, displz
        real(kind=realtype) :: tnew, told
        real(kind=realtype) :: anglex, angley, anglez, dx, dy, dz
        real(kind=realtype) :: xix, xiy, xiz, etax, etay, etaz
        real(kind=realtype) :: zetax, zetay, zetaz, xp, yp, zp, t
        real(kind=realtype) :: phi, cosphi, sinphi, eta, zeta
 
        real(kind=realtype), dimension(3) :: rotationpoint
        real(kind=realtype), dimension(3, 3) :: rotationmatrix
 
        ! Compute the displacements due to the rigid motion of the mesh.
 
        displx = zero
        disply = zero
        displz = zero
 
        ! Determine the time values of the old and new time level.
        ! It is assumed that the rigid body rotation of the mesh is only
        ! used when only 1 section is present.
 
        tnew = timeunsteady + timeunsteadyrestart
        told = tnew - dtadvance(1)
 
        ! Compute the rotation matrix of the rigid body rotation as
        ! well as the rotation point; the latter may vary in time due
        ! to rigid body translation.
 
        call rotmatrixrigidbody(tnew, told, rotationmatrix, rotationpoint)
 
        ! Loop over the number of local blocks.
 
        blockloop: do nn = 1, ndom
 
            ! Set the pointers for this block on the current level.
            ! Note that currentLevel must be used and not groundLevel,
            ! because groundLevel is 1 level too coarse when this routine
            ! is called in the full mg cycle.
 
            call setpointers(nn, currentlevel, sps)
            !
            !         The rigid body motion of the entire mesh.
            !
            ! First the rotation.
 
            do k = 1, kl
                do j = 1, jl
                    do i = 1, il
 
                        ! Determine the vector relative to the rotation point.
 
                        xp = x(i, j, k, 1) - rotationpoint(1)
                        yp = x(i, j, k, 2) - rotationpoint(2)
                        zp = x(i, j, k, 3) - rotationpoint(3)
 
                        ! Apply the transformation matrix to the vector (xp,yp,zp)
                        ! and set the new coordinates.
 
                        x(i, j, k, 1) = rotationmatrix(1, 1) * xp &
                                        + rotationmatrix(1, 2) * yp &
                                        + rotationmatrix(1, 3) * zp + rotationpoint(1)
                        x(i, j, k, 2) = rotationmatrix(2, 1) * xp &
                                        + rotationmatrix(2, 2) * yp &
                                        + rotationmatrix(2, 3) * zp + rotationpoint(2)
                        x(i, j, k, 3) = rotationmatrix(3, 1) * xp &
                                        + rotationmatrix(3, 2) * yp &
                                        + rotationmatrix(3, 3) * zp + rotationpoint(3)
                    end do
                end do
            end do
 
            ! Add the translation.
 
            do k = 1, kl
                do j = 1, jl
                    do i = 1, il
                        x(i, j, k, 1) = x(i, j, k, 1) + displx
                        x(i, j, k, 2) = x(i, j, k, 2) + disply
                        x(i, j, k, 3) = x(i, j, k, 3) + displz
                    end do
                end do
            end do
 
            !
            !         Determine whether the corresponding cgns block is a rotating
            !         block. If it is, apply the rotation.
            !         Note that now the section ID of the block is taken into
            !         account to allow for different periodic times per section.
            !
            if (cgnsdoms(nbkglobal)%rotatingFrameSpecified) then
 
                ! Compute the rotation angles.
 
                anglex = dtadvance(sectionid) * cgnsdoms(nbkglobal)%rotRate(1)
                angley = dtadvance(sectionid) * cgnsdoms(nbkglobal)%rotRate(2)
                anglez = dtadvance(sectionid) * cgnsdoms(nbkglobal)%rotRate(3)
 
                ! Compute the unit vector in the direction of the rotation
                ! axis, which will be called the xi-direction.
 
                t = one / max(eps, sqrt(anglex**2 + angley**2 + anglez**2))
                xix = t * anglex
                xiy = t * angley
                xiz = t * anglez
 
                ! Determine the rotation angle in xi-direction and its sine
                ! and cosine. Due to the definition of the xi-direction this
                ! angle will always be positive.
 
                phi = xix * anglex + xiy * angley + xiz * anglez
                cosphi = cos(phi)
                sinphi = sin(phi)
 
                ! Loop over the owned coordinates of this block.
 
                do k = 1, kl
                    do j = 1, jl
                        do i = 1, il
 
                            ! Compute the vector relative to center of rotation.
 
                            dx = x(i, j, k, 1) - cgnsdoms(nbkglobal)%rotCenter(1)
                            dy = x(i, j, k, 2) - cgnsdoms(nbkglobal)%rotCenter(2)
                            dz = x(i, j, k, 3) - cgnsdoms(nbkglobal)%rotCenter(3)
 
                            ! Compute the coordinates of the point p, which is the
                            ! closest point on the rotation axis.
 
                            t = dx * xix + dy * xiy + dz * xiz
                            xp = cgnsdoms(nbkglobal)%rotCenter(1) + t * xix
                            yp = cgnsdoms(nbkglobal)%rotCenter(2) + t * xiy
                            zp = cgnsdoms(nbkglobal)%rotCenter(3) + t * xiz
 
                            ! Determine the unit vector in eta direction, which
                            ! is defined from point p to the current point.
 
                            etax = x(i, j, k, 1) - xp
                            etay = x(i, j, k, 2) - yp
                            etaz = x(i, j, k, 3) - zp
 
                            eta = sqrt(etax**2 + etay**2 + etaz**2)
                            t = one / max(eps, eta)
 
                            etax = t * etax
                            etay = t * etay
                            etaz = t * etaz
 
                            ! Determine the unit vector in zeta-direction. This is
                            ! the cross product of the unit vectors in xi and in
                            ! eta-direction.
 
                            zetax = xiy * etaz - xiz * etay
                            zetay = xiz * etax - xix * etaz
                            zetaz = xix * etay - xiy * etax
 
                            ! Compute the new eta and zeta coordinates.
 
                            zeta = eta * sinphi
                            eta = eta * cosphi
 
                            ! Compute the new cartesian coordinates.
 
                            x(i, j, k, 1) = xp + eta * etax + zeta * zetax
                            x(i, j, k, 2) = yp + eta * etay + zeta * zetay
                            x(i, j, k, 3) = zp + eta * etaz + zeta * zetaz
 
                        end do
                    end do
                end do
 
            end if
 
        end do blockloop
 
    end subroutine updatecoorfinemesh
    subroutine alloccoorfinegrid
        !
        !       allocCoorFineGrid allocates the memory for all the coordinates
        !       of all local blocks. Also the memory for the derived data type
        !       used for the reading is allocated. If an interpolation must be
        !       performed for the time spectral method the variables of this
        !       IO type are allocated as well. For all other cases the pointer
        !       of the variables are set to the appropriate entry in flowDoms.
        !
        use constants
        use block, only: ndom, flowdoms
        use inputphysics, only: equationmode
        use inputtimespectral, only: ntimeintervalsspectral
        use iomodule, only: iovar
        use iteration, only: noldlevels, deforming_grid
        use partitionmod, only: interpolspectral, ngridsread
        use utils, only: terminate
        implicit none
        !
        !      Local variables.
        !
        integer :: ierr
 
        integer(kind=intType) :: nn, mm
        integer(kind=intType) :: il, jl, kl, ie, je, ke
 
        ! Loop over the local blocks and allocate the memory for the
        ! coordinates.
 
        blockloop: do nn = 1, ndom
 
            ! Some abbreviations of the block dimensions.
 
            il = flowdoms(nn, 1, 1)%il
            jl = flowdoms(nn, 1, 1)%jl
            kl = flowdoms(nn, 1, 1)%kl
 
            ie = flowdoms(nn, 1, 1)%ie
            je = flowdoms(nn, 1, 1)%je
            ke = flowdoms(nn, 1, 1)%ke
 
            ! Loop over the number of spectral modes and allocate the
            ! memory of the coordinates, single halo's included. The halo
            ! values will be initialized in the preprocessing part.
 
            do mm = 1, ntimeintervalsspectral
                allocate (flowdoms(nn, 1, mm)%x(0:ie, 0:je, 0:ke, 3), stat=ierr)
                if (ierr /= 0) &
                    call terminate("allocCoorFineGrid", &
                                   "Memory allocation failure for flowDoms%x")
 
                flowdoms(nn, 1, mm)%x = 0.0
 
                !allocate xInit for all time spectral intervals for mesh warping
 
                if (ierr /= 0) &
                    call terminate("allocCoorFineGrid", &
                                   "Memory allocation failure for flowDoms%xInit")
                !flowDoms(nn,1,mm)%xInit=0.0
                !for the first grid also allocate xPlus and xMinus for the
                !mesh warping verification...
            end do
 
            ! For a time accurate computation on deforming meshes, allocate
            ! the memory for the old coordinates. As this is not the time
            ! spectral mode, the third index of flowDoms is always 1.
 
            if (deforming_grid .and. equationmode == unsteady) then
 
                allocate (flowdoms(nn, 1, 1)%xOld(noldlevels, 0:ie, 0:je, 0:ke, 3), &
                          stat=ierr)
                if (ierr /= 0) &
                    call terminate("allocCoorFineGrid", &
                                   "Memory allocation failure for xOld")
            end if
 
            ! Added by HDN
            if (deforming_grid .and. equationmode == unsteady) then
 
                allocate (flowdoms(nn, 1, 1)%xALE(0:ie, 0:je, 0:ke, 3), &
                          stat=ierr)
                if (ierr /= 0) &
                    call terminate("allocCoorFineGrid", &
                                   "Memory allocation failure for xALE")
            end if
 
        end do blockloop
 
        ! Allocate the memory for IOVar.
 
        allocate (iovar(ndom, ngridsread), stat=ierr)
        if (ierr /= 0) &
            call terminate("allocCoorFineGrid", &
                           "Memory allocation failure for IOVar")
 
        ! Determine the equation mode we are solving and set the pointers
        ! of IOVar accordingly, or even allocate the memory, if needed.
 
        select case (equationmode)
 
        case (steady)
 
            ! Steady computation. Only one grid needs to be read.
            ! Loop over the number of blocks and set the pointer.
            ! No pointer offsets are needed for the coordinates.
 
            do nn = 1, ndom
                iovar(nn, 1)%pointerOffset = 0
                iovar(nn, 1)%w => flowdoms(nn, 1, 1)%x(1:, 1:, 1:, :)
            end do
 
            !===============================================================
 
        case (unsteady)
 
            ! Unsteady computation. The first set of coordinates should
            ! be stored in x, other sets (if present) in xOld.
            ! No pointer offsets are needed for the coordinates.
 
            do nn = 1, ndom
                iovar(nn, 1)%pointerOffset = 0
                iovar(nn, 1)%w => flowdoms(nn, 1, 1)%x(1:, 1:, 1:, :)
 
                do mm = 2, ngridsread
                    iovar(nn, mm)%pointerOffset = 0
                    iovar(nn, mm)%w => flowdoms(nn, 1, 1)%xOld(mm - 1, 1:, 1:, 1:, :)
                end do
            end do
 
            !===============================================================
 
        case (timespectral)
 
            ! Time spectral mode. A further check is required.
 
            testallociovar: if (interpolspectral .and. &
                                ngridsread > 1) then
 
                ! A restart is performed for a deforming mesh using a
                ! different number of time instances than the previous
                ! computation. Consequently the coordinates, or better
                ! the deformations, will be interpolated later on. Hence
                ! some additional storage is required for the coordinates
                ! to be read and thus the memory for the variables w of
                ! IOVar is allocated. No halo data is needed here.
                ! Note that for the 1st time instance the pointer is set
                ! to the coordinates of flowDoms, because these are not
                ! interpolated. Only the higher time instances are
                ! interpolated. No pointer offsets are needed for the
                ! coordinates.
 
                do nn = 1, ndom
                    il = flowdoms(nn, 1, 1)%il
                    jl = flowdoms(nn, 1, 1)%jl
                    kl = flowdoms(nn, 1, 1)%kl
 
                    iovar(nn, 1)%pointerOffset = 0
                    iovar(nn, 1)%w => flowdoms(nn, 1, 1)%x(1:, 1:, 1:, :)
 
                    do mm = 2, ngridsread
                        iovar(nn, mm)%pointerOffset = 0
                        allocate (iovar(nn, mm)%w(il, jl, kl, 3), stat=ierr)
                        if (ierr /= 0) &
                             call terminate("allocCoorFineGrid", &
                             "Memory allocation failure for &
                             &IOVar%w")
                    end do
                end do
 
                else testallociovar
 
                ! One of the following options is true.
                ! - The computation starts from scratch.
                ! - A restart is performed using a rigid grid, possibly
                !   moving. The number of time instances does not have
                !   to be the same.
                ! - A restart with a deforming mesh is performed with the
                !   same number of time instances.
                !
                ! In all these situations the pointers of IOVar are
                ! simply set to the coordinates of flowDoms.
 
                do nn = 1, ndom
                    do mm = 1, ngridsread
                        iovar(nn, mm)%pointerOffset = 0
                        iovar(nn, mm)%w => flowdoms(nn, 1, mm)%x(1:, 1:, 1:, :)
                    end do
                end do
 
            end if testallociovar
 
        end select
 
    end subroutine alloccoorfinegrid
    subroutine checkpartitioning(np, load_inbalance, face_inbalance)
 
        ! This subroutine runs the load balancing and partitioning algorithm
        ! to determine what the load balancing will be for a given number of
        ! procs np. The output is load_inbalance and face_inbalance.
 
        use constants
        use communication, only: nproc
        use partitionmod, only: ubvec
        use loadbalance, only: blockdistribution
        implicit none
 
        integer(kind=intType), intent(in) :: np
        integer(kind=intType) :: nproc_save
        real(kind=realtype), intent(out) :: load_inbalance, face_inbalance
 
        ! Note: This file follows mostly partitionAndReadGrid. See
        ! partitionAndReadGrid.f90 for more infromation
 
        ! Trick it into thinking we have np processors:
        nproc_save = nproc
        nproc = np
 
        call blockdistribution
 
        ! Restore the number of procs
        nproc = nproc_save
 
        ! Extract the inbalance info:
        load_inbalance = ubvec(1)
        face_inbalance = ubvec(2)
 
    end subroutine checkpartitioning
 
    subroutine determinesections
        !
        !       determineSections determines the number of sections, i.e.
        !       grid parts between sliding mesh interfaces, present in the
        !       entire grid.
        !
        use constants
        use block
        use cgnsgrid
        use communication
        use inputtimespectral
        use section
        use su_cgns
        use sorting, only: qsortintegers, bsearchintegers
        use utils, only: terminate
        implicit none
        !
        !      Local parameter, threshold for allowed angle difference between
        !                       the theoretical and true value, 0.1 degrees.
        !
        real(kind=realtype), parameter :: threshold = 0.1_realtype
        !
        !      Local variables.
        !
        integer :: ierr
 
        integer(kind=intType) :: nn, mm, ii, jj
        integer(kind=intType) :: nLevel, nSlices, slideID
        integer(kind=intType), dimension(cgnsNDom) :: sectionID, sorted
        integer(kind=intType), dimension(cgnsNsliding, 2) :: secSliding
 
        real(kind=realtype) :: costheta, cosphi, cospsi
        real(kind=realtype) :: sintheta, sinphi, sinpsi
        real(kind=realtype) :: r11, r12, r13, r21, r22, r23
        real(kind=realtype) :: r31, r32, r33
        real(kind=realtype) :: d1, d2, a1, a0, lamr, lami, angle, dangle
 
        logical :: situationChanged
 
        ! Initialize sectionID to the cgns block number.
 
        do nn = 1, cgnsndom
            sectionid(nn) = nn
        end do
 
        ! Loop to determine the highest cgns block id in every section.
 
        loophighestblock: do
 
            ! Initialize situationChanged to .false.
 
            situationchanged = .false.
 
            ! Loop over the internal block faces of the blocks.
 
            do nn = 1, cgnsndom
 
                ! First the 1 to 1 block connectivities.
 
                do mm = 1, cgnsdoms(nn)%n1to1
 
                    ! Store the neighboring block a bit easier and change the
                    ! value of sectionID if the neighboring block has a
                    ! higher value. In that case situationChanged is .true.
 
                    ii = cgnsdoms(nn)%conn1to1(mm)%donorBlock
                    if (sectionid(ii) > sectionid(nn)) then
                        sectionid(nn) = sectionid(ii)
                        situationchanged = .true.
                    end if
 
                end do
 
                ! No general connectivities yet.
 
            end do
 
            ! Criterion to exit the loop.
 
            if (.not. situationchanged) exit
 
        end do loophighestblock
 
        ! Copy sectionID in sorted and sort sorted in increasing order.
 
        do nn = 1, cgnsndom
            sorted(nn) = sectionid(nn)
        end do
 
        call qsortintegers(sorted, cgnsndom)
 
        ! Determine the number of sections. Note there is at least one
        ! section, because there is at least one block in the grid.
 
        nsections = 1
        do nn = 2, cgnsndom
            if (sorted(nn) > sorted(nsections)) then
                nsections = nsections + 1
                sorted(nsections) = sorted(nn)
            end if
        end do
 
        ! Determine the sections to which the owned blocks belong.
 
        nlevel = ubound(flowdoms, 2)
        do nn = 1, ndom
 
            ! Determine the corresponding cgns block id and search its
            ! section id in sorted.
 
            mm = flowdoms(nn, 1, 1)%cgnsBlockID
            mm = bsearchintegers(sectionid(mm), sorted)
 
            if (debug) then
                if (mm == 0) call terminate("determineSections", &
                                            "Entry not found in sorted.")
            end if
 
            ! Set the section id for all grid levels for all spectral
            ! time intervals.
 
            do jj = 1, ntimeintervalsspectral
                do ii = 1, nlevel
                    flowdoms(nn, ii, jj)%sectionID = mm
                end do
            end do
 
        end do
 
        ! Allocate the memory for sections.
 
        allocate (sections(nsections), stat=ierr)
        if (ierr /= 0) &
            call terminate("determineSections", &
                           "Memory allocation failure for sections.")
 
        ! Initialize the number of slices for each of the sections to 1,
        ! periodic and rotating to .false. and rotCenter, rotAxis
        ! and rotRate to zero.
 
        do nn = 1, nsections
            sections(nn)%nSlices = 1
            sections(nn)%periodic = .false.
            sections(nn)%rotating = .false.
            sections(nn)%rotCenter = zero
            sections(nn)%rotAxis = zero
            sections(nn)%rotRate = zero
        end do
 
        ! Determine the number of slices and the periodic transformation
        ! for the sections.
 
        loopcgnsdom: do nn = 1, cgnsndom
 
            ! Search for the corresponding section.
 
            ii = bsearchintegers(sectionid(nn), sorted(1:nsections))
            if (debug) then
                if (ii == 0) call terminate("determineSections", &
                                            "Entry not found in sorted.")
            end if
 
            ! It is assumed that periodic info is correct. So if this
            ! section has already been treated, there is no need to do
            ! it again.
 
            if (sections(ii)%periodic) cycle
 
            ! Loop over the 1 to 1 subfaces of the cgns block and try to
            ! find a periodic one.
 
            do mm = 1, cgnsdoms(nn)%n1to1
                if (cgnsdoms(nn)%conn1to1(mm)%periodic) exit
            end do
 
            ! Continue with the next block if this block does not have
            ! periodic subfaces.
 
            if (mm > cgnsdoms(nn)%n1to1) cycle
 
            ! Subface mm is a periodic one. Set periodic to .true.
 
            sections(ii)%periodic = .true.
 
            ! Set the rotation axis of the section to the rotation
            ! angles of the periodic transformation. This may be
            ! overwritten later on using the rotation rate, but for
            ! some cases this is the only rotation information present.
 
            sections(ii)%rotAxis = cgnsdoms(nn)%conn1to1(mm)%rotationAngles
 
            ! Construct the rotation matrix, where it is assumed 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(cgnsdoms(nn)%conn1to1(mm)%rotationAngles(1))
            sintheta = sin(cgnsdoms(nn)%conn1to1(mm)%rotationAngles(1))
 
            cosphi = cos(cgnsdoms(nn)%conn1to1(mm)%rotationAngles(2))
            sinphi = sin(cgnsdoms(nn)%conn1to1(mm)%rotationAngles(2))
 
            cospsi = cos(cgnsdoms(nn)%conn1to1(mm)%rotationAngles(3))
            sinpsi = sin(cgnsdoms(nn)%conn1to1(mm)%rotationAngles(3))
 
            r11 = cosphi * cospsi
            r21 = cosphi * sinpsi
            r31 = -sinphi
 
            r12 = sintheta * sinphi * cospsi - costheta * sinpsi
            r22 = sintheta * sinphi * sinpsi + costheta * cospsi
            r32 = sintheta * cosphi
 
            r13 = costheta * sinphi * cospsi + sintheta * sinpsi
            r23 = costheta * sinphi * sinpsi - sintheta * cospsi
            r33 = costheta * cosphi
 
            ! Store the rotation matrix, rotation center and translation
            ! vector for this section.
 
            sections(ii)%rotCenter = &
                cgnsdoms(nn)%conn1to1(mm)%rotationCenter
            sections(ii)%translation = &
                cgnsdoms(nn)%conn1to1(mm)%translation
 
            sections(ii)%rotMatrix(1, 1) = r11
            sections(ii)%rotMatrix(2, 1) = r21
            sections(ii)%rotMatrix(3, 1) = r31
 
            sections(ii)%rotMatrix(1, 2) = r12
            sections(ii)%rotMatrix(2, 2) = r22
            sections(ii)%rotMatrix(3, 2) = r32
 
            sections(ii)%rotMatrix(1, 3) = r13
            sections(ii)%rotMatrix(2, 3) = r23
            sections(ii)%rotMatrix(3, 3) = r33
 
            ! Determine the coefficients of lambda and lambda^2 of the
            ! characteristic polynomial.
 
            d2 = -r11 - r22 - r33
            d1 = r11 * r22 + r11 * r33 + r22 * r33 - r12 * r21 - r13 * r31 - r23 * r32
 
            ! Make use of the fact that one eigenvalue of the transformation
            ! matrix is 1 and determine the coefficients of the quadratic
            ! equation for the other two eigenvalues.
 
            a1 = d2 + one
            a0 = a1 + d1
 
            ! Determine the real and imaginary part of the two eigenvalues.
            ! Neglect the factor 1/2 here.
 
            lamr = -a1
            lami = sqrt(abs(a1 * a1 - four * a0))
 
            ! Determine the angle in the imaginary plane. Due to the
            ! positive definition of lami, this angle will be between
            ! 0 and pi. Take care of the exceptional case that the angle
            ! is zero.
 
            angle = atan2(lami, lamr)
            if (angle == zero) angle = two * pi
 
            ! Determine the number of slices.
 
            nslices = nint(two * pi / angle)
 
            ! Determine the angle difference in degrees between
            ! nSlices*angle and a complete rotation. If this is larger than
            ! the threshold processor 0 will print an error message and exit.
 
            dangle = abs(180.0_realtype * (two * pi - nslices * angle) / pi)
 
            if (dangle >= threshold) then
                if (myid == 0) &
                     call terminate("determineSections", &
                     "Periodic angle not a integer divide of &
                     &360 degrees")
                call mpi_barrier(adflow_comm_world, ierr)
            end if
 
            ! Set the number of slices for this section.
 
            sections(ii)%nSlices = nslices
 
        end do loopcgnsdom
 
        ! Again loop over the number of block of the original mesh,
        ! but now determine whether or not the section is rotating.
 
        do nn = 1, cgnsndom
 
            ! If the block is rotating, copy that information to
            ! the corresponding section. If the section is not
            ! periodic, also set the rotation center.
 
            if (cgnsdoms(nn)%rotatingFrameSpecified) then
 
                ii = bsearchintegers(sectionid(nn), sorted(1:nsections))
                sections(ii)%rotating = .true.
                sections(ii)%rotAxis = cgnsdoms(nn)%rotRate
                sections(ii)%rotRate = cgnsdoms(nn)%rotRate
 
                if (.not. sections(ii)%periodic) &
                    sections(ii)%rotCenter = cgnsdoms(nn)%rotCenter
 
            end if
        end do
 
        ! Determine the two sections for every sliding mesh
        ! interface.
 
        secsliding = 0
        do nn = 1, cgnsndom
            do mm = 1, cgnsdoms(nn)%nBocos
                if (cgnsdoms(nn)%bocoInfo(mm)%actualFace .and. &
                    cgnsdoms(nn)%bocoInfo(mm)%BCType == slidinginterface) then
 
                    ! Boundary face is part of a sliding mesh interface.
                    ! Determine the ID of the interface.
 
                    slideid = abs(cgnsdoms(nn)%bocoInfo(mm)%slidingID)
 
                    ! Determine the section to which this block belongs and
                    ! store its id for this sliding interface.
 
                    ii = bsearchintegers(sectionid(nn), sorted(1:nsections))
                    if (secsliding(slideid, 1) == 0) then
                        secsliding(slideid, 1) = ii
                    else if (secsliding(slideid, 1) /= ii) then
                        secsliding(slideid, 2) = ii
                    end if
 
                end if
            end do
        end do
 
        ! Loop over the sliding mesh interfaces to set the rotation axis
        ! for non-rotating sections.
 
        do ii = 1, cgnsnsliding
 
            ! Store the two sections id's a bit easier.
 
            mm = secsliding(ii, 1)
            nn = secsliding(ii, 2)
 
            ! Print an error message if both sections are not rotating.
 
            if ((.not. sections(mm)%rotating) .and. &
                (.not. sections(nn)%rotating)) then
                if (myid == 0) &
                     call terminate("determineSections", &
                     "Encountered sliding interface between &
                     &two non-rotating sections")
                call mpi_barrier(adflow_comm_world, ierr)
            end if
 
            ! Set the rotation axis if section mm is not rotating.
            ! If it is not periodic also set the rotation point.
 
            if (.not. sections(mm)%rotating) then
                sections(mm)%rotAxis = sections(nn)%rotAxis
                if (.not. sections(mm)%periodic) &
                    sections(mm)%rotCenter = sections(nn)%rotCenter
            end if
 
            ! Idem for section nn.
 
            if (.not. sections(nn)%rotating) then
                sections(nn)%rotAxis = sections(mm)%rotAxis
                if (.not. sections(nn)%periodic) &
                    sections(nn)%rotCenter = sections(mm)%rotCenter
            end if
 
        end do
 
        ! Determine the unit rotation axis for the sections.
 
        do nn = 1, nsections
            d1 = one / max(eps, sqrt(sections(nn)%rotAxis(1)**2 &
                                     + sections(nn)%rotAxis(2)**2 &
                                     + sections(nn)%rotAxis(3)**2))
 
            sections(nn)%rotAxis(1) = d1 * sections(nn)%rotAxis(1)
            sections(nn)%rotAxis(2) = d1 * sections(nn)%rotAxis(2)
            sections(nn)%rotAxis(3) = d1 * sections(nn)%rotAxis(3)
        end do
 
    end subroutine determinesections
 
end module partitioning