initializeFlow.F90

Go to the documentation of this file.

module initializeflow
 
    use constants, only: inttype, realtype, maxstringlen
 
    implicit none
    save
 
contains
 
    subroutine referencestate
        !
        !       The original version has been nuked since the computations are
        !       no longer necessary when calling from python
        !       This is the most compliclated routine in all of ADflow. It is
        !       stupidly complicated. This is most likely the reason your
        !       derivatives are wrong. You don't understand this routine
        !       and its effects.
        !       This routine *requries* the following as input:
        !       Mach, pInfDim, TInfDim, rhoInfDim, rGasDim (machCoef non-SA
        !        turbulence only)
        !       Optionally, pRef, rhoRef and Tref are used if they are
        !       are non-negative. This only happens when you want the equations
        !       normalized by values other than the freestream
        !      * This routine computes as output:
        !      *   muInfDim, (unused anywhere in code)
        !         pRef, rhoRef, Tref, muRef, timeRef ('dimensional' reference)
        !         pInf, pInfCorr, rhoInf, uInf, rGas, muInf, gammaInf and wInf
        !         (Non-dimensionalized values used in actual computations)
        !
        use constants
        use paramturb
        use inputphysics, only: equations, mach, machcoef, &
                                musuthdim, tsuthdim, veldirfreestream, &
                                rgasdim, ssuthdim, eddyvisinfratio, turbmodel, turbintensityinf
        use flowvarrefstate, only: pinfdim, tinfdim, rhoinfdim, &
                                   muinfdim, &
                                   pref, rhoref, tref, muref, timeref, uref, href, &
                                   pinf, pinfcorr, rhoinf, uinf, rgas, muinf, gammainf, winf, &
                                   nw, nwf, kpresent, winf
        use flowutils, only: computegamma, etot
        use turbutils, only: sanuknowneddyratio
        implicit none
 
        integer(kind=intType) :: sps, nn, mm, ierr
        real(kind=realtype) :: gm1, ratio
        real(kind=realtype) :: nuinf, ktmp, uinf2
        real(kind=realtype) :: vinf, zinf, tmp1(1), tmp2(1)
 
        ! Compute the dimensional viscosity from Sutherland's law
        muinfdim = musuthdim &
                   * ((tsuthdim + ssuthdim) / (tinfdim + ssuthdim)) &
                   * ((tinfdim / tsuthdim)**1.5_realtype)
 
        ! Set the reference values. They *COULD* be different from the
        ! free-stream values for an internal flow simulation. For now,
        ! we just use the actual free stream values.
        pref = pinfdim
        tref = tinfdim
        rhoref = rhoinfdim
 
        ! Compute the value of muRef, such that the nonDimensional
        ! equations are identical to the dimensional ones.
        ! Note that in the non-dimensionalization of muRef there is
        ! a reference length. However this reference length is 1.0
        ! in this code, because the coordinates are converted to
        ! meters.
 
        muref = sqrt(pref * rhoref)
 
        ! Compute timeRef for a correct nonDimensionalization of the
        ! unsteady equations. Some story as for the reference viscosity
        ! concerning the reference length.
 
        timeref = sqrt(rhoref / pref)
        href = pref / rhoref
        uref = sqrt(href)
 
        ! Compute the nonDimensional pressure, density, velocity,
        ! viscosity and gas constant.
 
        pinf = pinfdim / pref
        rhoinf = rhoinfdim / rhoref
        uinf = mach * sqrt(gammainf * pinf / rhoinf)
        rgas = rgasdim * rhoref * tref / pref
        muinf = muinfdim / muref
        tmp1(1) = tinfdim
        call computegamma(tmp1, tmp2, 1)
        gammainf = tmp2(1)
 
        ! ----------------------------------------
        !      Compute the final wInf
        ! ----------------------------------------
 
        ! Allocate the memory for wInf if necessary
#ifndef USE_TAPENADE
        if (allocated(winf)) deallocate (winf)
        allocate (winf(nw), stat=ierr)
#endif
 
        ! zero out the winf first
        winf(:) = zero
 
        ! Set the reference value of the flow variables, except the total
        ! energy. This will be computed at the end of this routine.
 
        winf(irho) = rhoinf
        winf(ivx) = uinf * veldirfreestream(1)
        winf(ivy) = uinf * veldirfreestream(2)
        winf(ivz) = uinf * veldirfreestream(3)
 
        ! Compute the velocity squared based on MachCoef. This gives a
        ! better indication of the 'speed' of the flow so the turubulence
        ! intensity ration is more meaningful especially for moving
        ! geometries. (Not used in SA model)
 
        uinf2 = machcoef * machcoef * gammainf * pinf / rhoinf
 
        ! Set the turbulent variables if transport variables are to be
        ! solved. We should be checking for RANS equations here,
        ! however, this code is included in block res. The issue is
        ! that for frozen turbulence (or ANK jacobian) we call the
        ! block_res with equationType set to Laminar even though we are
        ! actually solving the rans equations. The issue is that, the
        ! freestream turb variables will be changed to zero, thus
        ! changing the solution. Insteady we check if nw > nwf which
        ! will accomplish the same thing.
 
        if (nw > nwf) then
 
            nuinf = muinf / rhoinf
 
            select case (turbmodel)
 
            case (spalartallmaras, spalartallmarasedwards)
 
                winf(itu1) = sanuknowneddyratio(eddyvisinfratio, nuinf)
 
                !=============================================================
 
            case (komegawilcox, komegamodified, mentersst)
 
                winf(itu1) = 1.5_realtype * uinf2 * turbintensityinf**2
                winf(itu2) = winf(itu1) / (eddyvisinfratio * nuinf)
 
                !=============================================================
 
            case (ktau)
 
                winf(itu1) = 1.5_realtype * uinf2 * turbintensityinf**2
                winf(itu2) = eddyvisinfratio * nuinf / winf(itu1)
 
                !=============================================================
 
            case (v2f)
 
                winf(itu1) = 1.5_realtype * uinf2 * turbintensityinf**2
                winf(itu2) = 0.09_realtype * winf(itu1)**2 &
                             / (eddyvisinfratio * nuinf)
                winf(itu3) = 0.666666_realtype * winf(itu1)
                winf(itu4) = 0.0_realtype
 
            end select
 
        end if
 
        ! Set the value of pInfCorr. In case a k-equation is present
        ! add 2/3 times rho*k.
 
        pinfcorr = pinf
        if (kpresent) pinfcorr = pinf + two * third * rhoinf * winf(itu1)
 
        ! Compute the free stream total energy.
 
        ktmp = zero
        if (kpresent) ktmp = winf(itu1)
        vinf = zero
        zinf = zero
        call etot(rhoinf, uinf, vinf, zinf, pinfcorr, ktmp, &
                  winf(irhoe), kpresent)
 
    end subroutine referencestate
 
    ! ----------------------------------------------------------------------
    !                                                                      |
    !                    No Tapenade Routine below this line               |
    !                                                                      |
    ! ----------------------------------------------------------------------
 
#ifndef  USE_TAPENADE
    subroutine infchangecorrection(oldWinf, correctionTol, correctionType)
        ! Adjust the flow states to a change in wInf
        use constants
        use blockpointers, only: il, jl, kl, w, ndom, d2wall
        use flowvarrefstate, only: winf, nwf, nw
        use inputphysics, only: equations
        use inputtimespectral, only: ntimeintervalsspectral
        use haloexchange, only: whalo2
        use flowutils, only: adjustinflowangle
        use oversetdata, only: oversetpresent
        use iteration, only: currentlevel
        use turbbcroutines, only: applyallturbbcthisblock, bcturbtreatment
        use bcroutines, only: applyallbc_block
        use utils, only: setpointers, mynorm2, getrotmatrix
        use communication, only: myid
        implicit none
 
        real(kind=realtype), intent(in), dimension(nwf) :: oldwinf
        real(kind=realtype), intent(in) :: correctiontol
        character*(*), intent(in) :: correctionType
        integer(kind=intType) :: sps, nn, i, j, k, l
        real(kind=realtype) :: deltawinf(nwf)
 
        ! variables for rotation update
        real(kind=realtype), dimension(3) :: vec1, vec2, vcell
        real(kind=realtype), dimension(3, 3) :: rotmat
        real(kind=realtype) :: mag1, mag2
 
        ! Make sure we have the updated wInf
        call adjustinflowangle()
        call referencestate
 
        ! first check if the state changed enough to warrant a correction. This tolerance
        ! can be adjusted from python. Default is 1e-12.
        deltawinf = winf(1:nwf) - oldwinf(1:nwf)
 
        if (mynorm2(deltawinf) < correctiontol) then
            ! The change deltaWinf is so small, (or zero) don't do the
            ! update and just return. This will save some time when the
            ! solver is called with the same AP conditions multiple times,
            ! such as during a GS AS solution
            return
        end if
 
        ! The state changed enough so we will run a correction update.
        ! Check which correction type we are using
        if (correctiontype .eq. "offset") then
 
            ! Loop over all the blocks, adding the offset between the old and new Winf
            do sps = 1, ntimeintervalsspectral
                do nn = 1, ndom
                    call setpointers(nn, 1_inttype, sps)
                    do k = 2, kl
                        do j = 2, jl
                            do i = 2, il
                                do l = 1, nwf
                                    w(i, j, k, l) = w(i, j, k, l) + deltawinf(l)
                                end do
                            end do
                        end do
                    end do
                    call applyallbc_block(.true.)
                end do
            end do
 
        else if (correctiontype .eq. "rotate") then
 
            ! Rotate the flow velocities by the rotation change in Winf.
            ! We also scale the velocity magnitudes based on the magnitude change.
            ! vec1 is the old free stream velocity vector, vec2 is the new one.
            vec1 = oldwinf(ivx:ivz)
            vec2 = winf(ivx:ivz)
 
            ! get the general rotation matrix
            call getrotmatrix(vec1, vec2, rotmat)
 
            ! Compute the magnitude of vec1 and vec2 so that we can adjust the magnitude of the velocity
            mag1 = mynorm2(vec1)
            mag2 = mynorm2(vec2)
 
            ! Loop over all the blocks, rotate each cell's velocity
            do sps = 1, ntimeintervalsspectral
                do nn = 1, ndom
                    call setpointers(nn, 1_inttype, sps)
                    do k = 2, kl
                        do j = 2, jl
                            do i = 2, il
                                ! pick out the cell velocity
                                vcell = w(i, j, k, ivx:ivz)
                                ! rotate and put back
                                w(i, j, k, ivx:ivz) = matmul(rotmat, vcell)
                                ! scale the magnitudes
                                w(i, j, k, ivx:ivz) = w(i, j, k, ivx:ivz) * mag2 / mag1
                                ! we also add the offset to the density and energy
                                w(i, j, k, irho) = w(i, j, k, irho) + deltawinf(irho)
                                w(i, j, k, irhoe) = w(i, j, k, irhoe) + deltawinf(irhoe)
                            end do
                        end do
                    end do
                    call applyallbc_block(.true.)
                end do
            end do
 
        end if
 
        ! Exchange values
        call whalo2(currentlevel, 1_inttype, nw, .true., .true., .true.)
 
        ! Need to re-apply the BCs. The reason is that BC halos behind
        ! interpolated cells need to be recomputed with their new
        ! interpolated values from actual compute cells. Only needed for
        ! overset.
        if (oversetpresent) then
            do sps = 1, ntimeintervalsspectral
                do nn = 1, ndom
                    call setpointers(nn, 1, sps)
                    if (equations == ransequations) then
                        call bcturbtreatment
                        call applyallturbbcthisblock(.true.)
                    end if
                    call applyallbc_block(.true.)
                end do
            end do
        end if
 
    end subroutine infchangecorrection
 
    ! Section out the BCdata setup so that it can by called from python when needed
    subroutine updatebcdataalllevels
        ! sets the prescribed boundary data from the CGNS arrays
 
        use constants
        use iteration, only: groundlevel
        use bcdata, only: setbcdatafinegrid, setbcdatacoarsegrid
        implicit none
 
        ! Allocate the memory for the prescribed boundary data at the
        ! boundary faces and determine the data for both the fine grid.
 
        groundlevel = 1
 
        ! Determine the reference state.
        call referencestate
 
        call setbcdatafinegrid(.true.)
 
        ! Determine the prescribed data on the coarse grid levels
        ! by interpolation.
#ifndef USE_TAPENADE
        call setbcdatacoarsegrid
#endif
 
    end subroutine updatebcdataalllevels
 
    subroutine initflow
        !
        !       initFlow allocates the
        !       memory for and initializes the flow variables. In case a
        !       restart is performed the owned variables are read from the
        !       previous solution file(s).
        !
        use constants
        use block, only: flowdoms
        use inputtimespectral, only: ntimeintervalsspectral
        use variablereading, only: halosread
 
        implicit none
        !
        !      Local variables.
        !
        integer :: ierr
 
        integer(kind=intType) :: sps, level, nLevels
 
        ! Determine the number of multigrid levels.
 
        nlevels = ubound(flowdoms, 2)
 
        ! As some boundary conditions can be treated in multiple ways,
        ! some memory allocated must be released again.
 
        call releaseextramembcs
 
        ! Determine for the time spectral mode the matrices for the
        ! time derivatives.
        call timespectralmatrices
 
        ! Loop over the number of spectral solutions to allocate
        ! the memory for the w-variables and p on the fine grid.
        do sps = 1, ntimeintervalsspectral
            call allocmemflovarpart1(sps, 1_inttype)
        end do
 
        ! Allocate the memory for the solution variables on the coarse
        ! grid levels and the memory for the dependent flow variables,
        ! residuals, etc, on all multigrid levels.
        do sps = 1, ntimeintervalsspectral
            call allocmemflovarpart2(sps, 1_inttype)
 
            do level = 2, nlevels
                call allocmemflovarpart1(sps, level)
                call allocmemflovarpart2(sps, level)
            end do
        end do
 
        ! Initialize free stream field
        call initflowfield
 
        ! Initialize the dependent flow variables and the halo values.
 
        call initdepvarandhalos(halosread)
 
    end subroutine initflow
 
    subroutine allocmemflovarpart1(sps, level)
        !
        !       allocMemFlovarPart1 allocates the memory for the flow
        !       variables w and p for all the blocks on the given multigrid
        !       level and spectral solution sps.
        !
        use constants
        use block, only: flowdoms, ndom
        use flowvarrefstate, only: nw, nwf, nt1, nt2
        use inputphysics, only: equationmode, gammaconstant
        use inputunsteady, only: timeintegrationscheme
        use inputiteration, only: mgstartlevel, turbtreatment
        use iteration, only: noldlevels
        use utils, only: terminate
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: sps, level
        !
        !      Local variables.
        !
        integer :: ierr
 
        integer(kind=intType) :: nn
        integer(kind=intType) :: il, jl, kl, ie, je, ke, ib, jb, kb
 
        ! Loop over the domains.
 
        domains: do nn = 1, ndom
 
            ! Store some dimensions a bit easier.
 
            il = flowdoms(nn, level, sps)%il
            jl = flowdoms(nn, level, sps)%jl
            kl = flowdoms(nn, level, sps)%kl
 
            ie = flowdoms(nn, level, sps)%ie
            je = flowdoms(nn, level, sps)%je
            ke = flowdoms(nn, level, sps)%ke
 
            ib = flowdoms(nn, level, sps)%ib
            jb = flowdoms(nn, level, sps)%jb
            kb = flowdoms(nn, level, sps)%kb
 
            ! Allocate the memory for the independent variables.
            ! Memory is allocated for the turbulent variables (if any) if
            ! the current level is smaller or equal to the multigrid start
            ! level or if the turbulent transport equations are solved in
            ! a coupled manner.
 
            if (level <= mgstartlevel .or. turbtreatment == coupled) then
                allocate (flowdoms(nn, level, sps)%w(0:ib, 0:jb, 0:kb, 1:nw), &
                          stat=ierr)
            else
                allocate (flowdoms(nn, level, sps)%w(0:ib, 0:jb, 0:kb, 1:nwf), &
                          stat=ierr)
            end if
            if (ierr /= 0) &
                call terminate("allocMemFlovarPart1", &
                               "Memory allocation failure for w")
 
            ! Alloc mem for nodal gradients
            allocate (flowdoms(nn, level, sps)%ux(il, jl, kl), stat=ierr)
            allocate (flowdoms(nn, level, sps)%uy(il, jl, kl), stat=ierr)
            allocate (flowdoms(nn, level, sps)%uz(il, jl, kl), stat=ierr)
 
            allocate (flowdoms(nn, level, sps)%vx(il, jl, kl), stat=ierr)
            allocate (flowdoms(nn, level, sps)%vy(il, jl, kl), stat=ierr)
            allocate (flowdoms(nn, level, sps)%vz(il, jl, kl), stat=ierr)
 
            allocate (flowdoms(nn, level, sps)%wx(il, jl, kl), stat=ierr)
            allocate (flowdoms(nn, level, sps)%wy(il, jl, kl), stat=ierr)
            allocate (flowdoms(nn, level, sps)%wz(il, jl, kl), stat=ierr)
 
            allocate (flowdoms(nn, level, sps)%qx(il, jl, kl), stat=ierr)
            allocate (flowdoms(nn, level, sps)%qy(il, jl, kl), stat=ierr)
            allocate (flowdoms(nn, level, sps)%qz(il, jl, kl), stat=ierr)
 
            ! Allocate memory for the pressure.
            allocate (flowdoms(nn, level, sps)%p(0:ib, 0:jb, 0:kb), stat=ierr)
            if (ierr /= 0) &
                call terminate("allocMemFlovarPart1", &
                               "Memory allocation failure for p")
 
            ! Allocate memory for the speed of sound squared
            allocate (flowdoms(nn, level, sps)%aa(0:ib, 0:jb, 0:kb), stat=ierr)
            if (ierr /= 0) &
                call terminate("allocMemFlovarPart1", &
                               "Memory allocation failure for p")
 
            ! The eddy viscosity for eddy viscosity models.
            ! Although a dependent variable, it is allocated on all grid
            ! levels, because the eddy viscosity might be frozen in the
            ! multigrid.
 
            ! Always allocate rev due to reverse mode AD- Peter Lyu. Also
            ! zero so that it doesn't affect laminar cases.
            allocate (flowdoms(nn, level, sps)%rev(0:ib, 0:jb, 0:kb), &
                      stat=ierr)
            flowdoms(nn, level, sps)%rev = zero
            if (ierr /= 0) &
                call terminate("allocMemFlovarPart1", &
                               "Memory allocation failure for rev")
            !endif
 
            ! If this is the finest grid some more memory must be allocated.
 
            fineleveltest: if (level == 1) then
 
                ! Allocate the memory for gamma and initialize it to
                ! the constant gamma value.
 
                allocate (flowdoms(nn, level, sps)%gamma(0:ib, 0:jb, 0:kb), &
                          stat=ierr)
                if (ierr /= 0) &
                    call terminate("allocMemFlovarPart1", &
                                   "Memory allocation failure for gamma.")
 
                flowdoms(nn, level, sps)%gamma = gammaconstant
 
                ! The laminar viscosity for viscous computations.
                ! Always allocate rlv due to reverse mode - Peter Lyu
                !if( viscous ) then
                allocate (flowdoms(nn, level, sps)%rlv(0:ib, 0:jb, 0:kb), &
                          stat=ierr)
                if (ierr /= 0) &
                    call terminate("allocMemFlovarPart1", &
                                   "Memory allocation failure for rlv")
                !endif
 
                ! The state vectors in the past for unsteady computations.
 
                if (equationmode == unsteady .and. &
                    timeintegrationscheme == bdf) then
                    allocate ( &
                        flowdoms(nn, level, sps)%wOld(noldlevels, 2:il, 2:jl, 2:kl, nw), &
                        stat=ierr)
                    if (ierr /= 0) &
                        call terminate("allocMemFlovarPart1", &
                                       "Memory allocation failure for wOld")
 
                    ! Initialize wOld to zero, such that it is initialized.
                    ! The actual values do not matter.
 
                    flowdoms(nn, level, sps)%wOld = zero
 
                    ! Added by HDN
                else if (equationmode == unsteady .and. &
                         timeintegrationscheme == md) then
                    allocate ( &
                        flowdoms(nn, level, sps)%wOld(noldlevels, 2:il, 2:jl, 2:kl, nw), &
                        stat=ierr)
                    if (ierr /= 0) &
                        call terminate("allocMemFlovarPart1", &
                                       "Memory allocation failure for wOld")
 
                    flowdoms(nn, level, sps)%wOld = zero
                end if
 
                ! If this is the 1st spectral solution (note that we are
                ! already on the finest grid) and the rans equations are
                ! solved, allocate the memory for the arrays used for the
                ! implicit boundary condition treatment. Normally this should
                ! only be allocated for RANS but the derivative calcs require
                ! these be allocated.
 
                sps1ranstest: if (sps == 1) then
                    allocate (flowdoms(nn, level, sps)%bmti1(je, ke, nt1:nt2, nt1:nt2), &
                              flowdoms(nn, level, sps)%bmti2(je, ke, nt1:nt2, nt1:nt2), &
                              flowdoms(nn, level, sps)%bmtj1(ie, ke, nt1:nt2, nt1:nt2), &
                              flowdoms(nn, level, sps)%bmtj2(ie, ke, nt1:nt2, nt1:nt2), &
                              flowdoms(nn, level, sps)%bmtk1(ie, je, nt1:nt2, nt1:nt2), &
                              flowdoms(nn, level, sps)%bmtk2(ie, je, nt1:nt2, nt1:nt2), &
                              flowdoms(nn, level, sps)%bvti1(je, ke, nt1:nt2), &
                              flowdoms(nn, level, sps)%bvti2(je, ke, nt1:nt2), &
                              flowdoms(nn, level, sps)%bvtj1(ie, ke, nt1:nt2), &
                              flowdoms(nn, level, sps)%bvtj2(ie, ke, nt1:nt2), &
                              flowdoms(nn, level, sps)%bvtk1(ie, je, nt1:nt2), &
                              flowdoms(nn, level, sps)%bvtk2(ie, je, nt1:nt2), &
                              stat=ierr)
                    if (ierr /= 0) &
                        call terminate("allocMemFlovarPart1", &
                                       "Memory allocation failure for bmti1, etc")
 
                end if sps1ranstest
            end if fineleveltest
 
        end do domains
 
    end subroutine allocmemflovarpart1
 
    !      ==================================================================
 
    subroutine allocmemflovarpart2(sps, level)
        !
        !       AllocMemFlovarPart2 allocates the memory for the dependent
        !       flow variables and iteration variables for all the blocks on
        !       the given multigrid level and spectral solution sps. Some
        !       variables are only allocated on the coarser grids, e.g. the
        !       multigrid forcing terms and the state vector upon entrance on
        !       the mg level. Other variables are only allocated on the finest
        !       mesh. These are typically dependent variables like laminar
        !       viscosity, or residuals, time step, etc. Exceptions are
        !       pressure and eddy viscosity. Although these are dependent
        !       variables, they are allocated on all grid levels.
        !
        use block
        use constants
        use flowvarrefstate
        use inputphysics
        use inputdiscretization
        use inputiteration
        use inputunsteady
        use iteration
        use utils, only: terminate
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: sps, level
        !
        !      Local variables.
        !
        integer :: ierr
 
        integer(kind=intType) :: nn, mm
        integer(kind=intType) :: il, jl, kl, ie, je, ke, ib, jb, kb
 
        ! Loop over the domains.
 
        domains: do nn = 1, ndom
 
            ! Store some dimensions a bit easier.
 
            il = flowdoms(nn, level, sps)%il
            jl = flowdoms(nn, level, sps)%jl
            kl = flowdoms(nn, level, sps)%kl
 
            ie = flowdoms(nn, level, sps)%ie
            je = flowdoms(nn, level, sps)%je
            ke = flowdoms(nn, level, sps)%ke
 
            ib = flowdoms(nn, level, sps)%ib
            jb = flowdoms(nn, level, sps)%jb
            kb = flowdoms(nn, level, sps)%kb
 
            ! Block is moving. Allocate the memory for s, sFaceI,
            ! sFaceJ and sFaceK.
 
            allocate (flowdoms(nn, level, sps)%s(ie, je, ke, 3), &
                      flowdoms(nn, level, sps)%sFaceI(0:ie, je, ke), &
                      flowdoms(nn, level, sps)%sFaceJ(ie, 0:je, ke), &
                      flowdoms(nn, level, sps)%sFaceK(ie, je, 0:ke), stat=ierr)
            if (ierr /= 0) &
                 call terminate("allocMemFlovarPart2", &
                 "Memory allocation failure for s, &
                 &sFaceI, sFaceJ and sFaceK.")
 
            ! Extra face velocities for ALE
            if (equationmode == unsteady .and. useale) then
                allocate ( &
                    flowdoms(nn, level, sps)%sVeloIALE(0:ie, je, ke, 3), &
                    flowdoms(nn, level, sps)%sVeloJALE(ie, 0:je, ke, 3), &
                    flowdoms(nn, level, sps)%sVeloKALE(ie, je, 0:ke, 3), &
                    flowdoms(nn, level, sps)%sFaceIALE(0:nalesteps, 0:ie, je, ke), &
                    flowdoms(nn, level, sps)%sFaceJALE(0:nalesteps, ie, 0:je, ke), &
                    flowdoms(nn, level, sps)%sFaceKALE(0:nalesteps, ie, je, 0:ke), stat=ierr)
                if (ierr /= 0) &
                     call terminate("allocMemFlovarPart2", &
                     "Memory allocation failure for &
                     &sVeloIALE, sVeloJALE and sVeloKALE; &
                     &sFaceIALE, sFaceJALE and sFaceKALE.")
            end if
 
            ! Test if we are on the finest mesh.
 
            fineleveltest: if (level == 1) then
 
                ! Allocate the memory that must always be allocated.
 
                allocate ( &
                    flowdoms(nn, level, sps)%dw(0:ib, 0:jb, 0:kb, 1:nw), &
                    flowdoms(nn, level, sps)%fw(0:ib, 0:jb, 0:kb, 1:nwf), &
                    flowdoms(nn, level, sps)%dtl(1:ie, 1:je, 1:ke), &
                    flowdoms(nn, level, sps)%radI(1:ie, 1:je, 1:ke), &
                    flowdoms(nn, level, sps)%radJ(1:ie, 1:je, 1:ke), &
                    flowdoms(nn, level, sps)%radK(1:ie, 1:je, 1:ke), &
                    flowdoms(nn, level, sps)%scratch(0:ib, 0:jb, 0:kb, 10), &
                    flowdoms(nn, level, sps)%shockSensor(0:ib, 0:jb, 0:kb), &
                    stat=ierr)
                if (ierr /= 0) &
                     call terminate("allocMemFlovarPart2", &
                     "Memory allocation failure for dw, fw, dwOld, fwOld, &
                     &gamma, dtl and the spectral radii.")
 
                ! Initialize dw and fw to zero.
 
                flowdoms(nn, level, sps)%dw = zero
                flowdoms(nn, level, sps)%fw = zero
 
                ! Extra variables for ALE
                if (equationmode == unsteady .and. useale) then
                    allocate ( &
                        flowdoms(nn, level, sps)%dwALE(0:nalesteps, 0:ib, 0:jb, 0:kb, 1:nw), &
                        flowdoms(nn, level, sps)%fwALE(0:nalesteps, 0:ib, 0:jb, 0:kb, 1:nwf), &
                        stat=ierr)
                    if (ierr /= 0) &
                        call terminate("allocMemFlovarPart2", &
                                       "Memory allocation failure for dwALE, fwALE.")
 
                    flowdoms(nn, level, sps)%dwALE = zero
                    flowdoms(nn, level, sps)%fwALE = zero
                end if
 
                ! Allocate the memory for the zeroth runge kutta stage
                allocate (flowdoms(nn, level, sps)%wn(2:il, 2:jl, 2:kl, 1:nwf), &
                          flowdoms(nn, level, sps)%pn(2:il, 2:jl, 2:kl), stat=ierr)
                if (ierr /= 0) &
                    call terminate("allocMemFlovarPart2", &
                                   "Memory allocation failure for wn and pn")
 
                ! For unsteady mode using Runge-Kutta schemes allocate the
                ! memory for dwOldRK.
 
                if (equationmode == unsteady .and. &
                    timeintegrationscheme == explicitrk) then
 
                    mm = nrkstagesunsteady - 1
                    allocate (flowdoms(nn, level, sps)%dwOldRK(mm, il, jl, kl, nw), &
                              stat=ierr)
                    if (ierr /= 0) &
                        call terminate("allocMemFlovarPart2", &
                                       "Memory allocation failure for dwOldRK.")
                end if
 
                else fineleveltest
 
                ! Coarser level. Allocate the memory for the multigrid
                ! forcing term and the state variables upon entry.
 
                allocate (flowdoms(nn, level, sps)%p1(1:ie, 1:je, 1:ke), &
                          flowdoms(nn, level, sps)%w1(1:ie, 1:je, 1:ke, 1:nwf), &
                          flowdoms(nn, level, sps)%wr(2:il, 2:jl, 2:kl, 1:nwf), &
                          stat=ierr)
                if (ierr /= 0) &
                     call terminate("allocMemFlovarPart2", &
                     "Memory allocation failure for p1, w1 &
                     &and wr")
 
                ! Initialize w1 and p1 to zero, just that they
                ! are initialized.
 
                flowdoms(nn, level, sps)%p1 = zero
                flowdoms(nn, level, sps)%w1 = zero
 
            end if fineleveltest
 
        end do domains
 
    end subroutine allocmemflovarpart2
 
    subroutine allocrestartfiles(nFiles)
        !
        !          Allocate memory for the restartfles                         *
        !          The array is populated from Python using setRestartFiles
        !          If memeory has been allocated for the array there exist at
        !          least one element in the array.
        !
        use constants
        use inputio, only: restartfiles
        use utils, only: terminate
        implicit none
        !
        !      Subroutine argument.
        !
        integer(kind=intType) :: nFiles
        !
        !      Local variables.
        !
        integer :: ierr
 
        if (allocated(restartfiles)) then
            deallocate (restartfiles)
        end if
 
        allocate (restartfiles(nfiles), stat=ierr)
        if (ierr /= 0) &
            call terminate("allocRestartFiles", &
                           "Memory allocation failure for restartFiles")
 
        ! Zero Array with empty strings
        restartfiles = ""
 
    end subroutine allocrestartfiles
 
    subroutine copyspectralsolution
        !
        !       copySpectralSolution copies the solution of the 1st spectral
        !       solution to all spectral solutions. This typically occurs when
        !       a for the spectral mode a restart is made from a steady or an
        !       unsteady solution. Possible rotation effects are taken into
        !       account for the velocity components.
        !
        use constants
        use block, only: flowdoms, ndom
        use flowvarrefstate, only: nw
        use inputtimespectral, only: ntimeintervalsspectral
        use iomodule, only: iovar
        use monitor, only: timeunsteadyrestart
        use section, only: sections, nsections
        use utils, only: rotmatrixrigidbody
        implicit none
        !
        !      Local variables.
        !
        integer(kind=intType) :: sps, spsm1, mm, nn, i, j, k, l
 
        real(kind=realtype) :: dt, tnew, told, tmp
        real(kind=realtype) :: theta, costheta, sintheta
 
        real(kind=realtype), dimension(3) :: rotpoint, uu
        real(kind=realtype), dimension(3) :: xt, yt, zt
        real(kind=realtype), dimension(3, 3) :: rotmat
 
        real(kind=realtype), dimension(nSections, 3, 3) :: rotmatsec
 
        ! Determine the rotation matrix from one spectral solution to the
        ! other for every section.
 
        sectionloop: do nn = 1, nsections
 
            ! Test if the section is rotating.
 
            testrotating: if (sections(nn)%rotating) then
 
                ! Section is rotating. Determine the angle between the
                ! spectral solutions and its sine and cosine.
 
                i = sections(nn)%nSlices * ntimeintervalsspectral
                theta = two * pi / real(i, realtype)
                costheta = cos(theta)
                sintheta = sin(theta)
 
                ! 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)
 
                ! 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.
 
                rotmatsec(nn, 1, 1) = xt(1) * xt(1) &
                                      + costheta * (yt(1) * yt(1) + zt(1) * zt(1))
                rotmatsec(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))
                rotmatsec(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))
 
                rotmatsec(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))
                rotmatsec(nn, 2, 2) = xt(2) * xt(2) &
                                      + costheta * (yt(2) * yt(2) + zt(2) * zt(2))
                rotmatsec(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))
 
                rotmatsec(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))
                rotmatsec(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))
                rotmatsec(nn, 3, 3) = xt(3) * xt(3) &
                                      + costheta * (yt(3) * yt(3) + zt(3) * zt(3))
 
                else testrotating
 
                ! Section is not rotating. Set the rotation matrix to the
                ! identity matrix.
 
                rotmatsec(nn, 1, 1) = one
                rotmatsec(nn, 1, 2) = zero
                rotmatsec(nn, 1, 3) = zero
 
                rotmatsec(nn, 2, 1) = zero
                rotmatsec(nn, 2, 2) = one
                rotmatsec(nn, 2, 3) = zero
 
                rotmatsec(nn, 3, 1) = zero
                rotmatsec(nn, 3, 2) = zero
                rotmatsec(nn, 3, 3) = one
 
            end if testrotating
 
        end do sectionloop
 
        ! Initialize told to timeUnsteadyRestart. This takes the
        ! possibility into account that the spectral mode is restarted
        ! from an unsteady computation. Although not likely this
        ! possibility is allowed and should therefore be taken into
        ! account. Anyway told corresponds to the time of the 1st
        ! spectral solution. Also determine the time step between the
        ! spectral solutions. Both told and dt are only used to determine
        ! the rigid body motion and if these are specified it is assumed
        ! that there is only one section present in the grid.
 
        told = timeunsteadyrestart
        dt = sections(1)%timePeriod &
             / real(ntimeintervalsspectral, realtype)
 
        ! Loop over the number of spectral modes, starting at 2.
 
        spectralloop: do sps = 2, ntimeintervalsspectral
 
            ! Determine the corresponding time for this spectral solution
            ! and store sps - 1 a bit easier.
 
            tnew = told + dt
            spsm1 = sps - 1
 
            ! Determine the rotation matrix and rotation point between the
            ! told and tnew for the rigid body rotation of the entire mesh.
            ! The rotation point is not needed for the transformation of the
            ! velocities, but rotMatrixRigidBody happens to compute it.
 
            call rotmatrixrigidbody(tnew, told, rotmat, rotpoint)
 
            ! Loop over the local number of blocks.
 
            domains: do nn = 1, ndom
 
                ! Store the section ID of this block a bit easier in mm.
 
                mm = flowdoms(nn, 1, 1)%sectionId
 
                ! Loop over the owned cells of this block. As the number of
                ! cells is identical for all spectral solutions, it does not
                ! matter which mode is taken for the upper dimensions.
 
                do k = 2, flowdoms(nn, 1, 1)%kl
                    do j = 2, flowdoms(nn, 1, 1)%jl
                        do i = 2, flowdoms(nn, 1, 1)%il
 
                            ! Step 1. Copy the solution variables w from
                            ! the previous spectral solution.
 
                            do l = 1, nw
                                iovar(nn, sps)%w(i, j, k, l) = iovar(nn, spsm1)%w(i, j, k, l)
                            end do
 
                            ! Step 2. Apply the rigid body motion rotation matrix
                            ! to the velocity. Use uu as a temporary storage.
 
                            uu(1) = rotmat(1, 1) * iovar(nn, sps)%w(i, j, k, ivx) &
                                    + rotmat(1, 2) * iovar(nn, sps)%w(i, j, k, ivy) &
                                    + rotmat(1, 3) * iovar(nn, sps)%w(i, j, k, ivz)
 
                            uu(2) = rotmat(2, 1) * iovar(nn, sps)%w(i, j, k, ivx) &
                                    + rotmat(2, 2) * iovar(nn, sps)%w(i, j, k, ivy) &
                                    + rotmat(2, 3) * iovar(nn, sps)%w(i, j, k, ivz)
 
                            uu(3) = rotmat(3, 1) * iovar(nn, sps)%w(i, j, k, ivx) &
                                    + rotmat(3, 2) * iovar(nn, sps)%w(i, j, k, ivy) &
                                    + rotmat(3, 3) * iovar(nn, sps)%w(i, j, k, ivz)
 
                            ! Step 3. Apply the rotation matrix of the section to
                            ! the velocity.
 
                            iovar(nn, sps)%w(i, j, k, ivx) = rotmatsec(mm, 1, 1) * uu(1) &
                                                             + rotmatsec(mm, 1, 2) * uu(2) &
                                                             + rotmatsec(mm, 1, 3) * uu(3)
 
                            iovar(nn, sps)%w(i, j, k, ivy) = rotmatsec(mm, 2, 1) * uu(1) &
                                                             + rotmatsec(mm, 2, 2) * uu(2) &
                                                             + rotmatsec(mm, 2, 3) * uu(3)
 
                            iovar(nn, sps)%w(i, j, k, ivz) = rotmatsec(mm, 3, 1) * uu(1) &
                                                             + rotmatsec(mm, 3, 2) * uu(2) &
                                                             + rotmatsec(mm, 3, 3) * uu(3)
                        end do
                    end do
                end do
 
            end do domains
 
            ! Set told to tnew for the next spectral solution.
 
            told = tnew
 
        end do spectralloop
 
    end subroutine copyspectralsolution
 
    subroutine determinesolfilenames
        !
        !       determineSolFileNames determines the number and names of the
        !       files that contain the solutions. For steady computations only
        !       one file must be present. For unsteady the situation is a
        !       little more complicated. It is attempted to read as many
        !       solutions as needed for a consistent restart. If not possible
        !       as many as possible solutions are read. For an unsteady
        !       computation the order will be reduced; for time spectral mode
        !       the solution will be interpolated.
        !
        use constants
        use communication, only: myid
        use inputio, only: restartfiles
        use inputphysics, only: equationmode
        use inputtimespectral, only: ntimeintervalsspectral
        use iteration, only: noldsolavail, oldsolwritten, noldlevels
        use variablereading, only: solfiles, nsolsread, interpolspectral, copyspectral
        use utils, only: terminate
        implicit none
        !
        !      Local variables
        !
        integer :: ierr
 
        integer(kind=intType) :: ii, nn
 
        character(len=7) :: integerString
        character(len=maxStringLen) :: tmpName
 
        ! Initialize copySpectral and interpolSpectral to .false.
 
        copyspectral = .false.
        interpolspectral = .false.
 
        ! Determine the desired number of files to 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 solution needs to be read.
            ! In case a list of restart files were provided in the python
            ! script, we force a read from only the first solution file.
 
            nsolsread = 1
            allocate (solfiles(nsolsread), stat=ierr)
            if (ierr /= 0) &
                call terminate("determineSolFileNames", &
                               "Memory allocation failure for solFiles")
 
            solfiles(1) = restartfiles(1)
 
            ! Check if the files can be opened, exit if that fails.
            call checksolfilenames()
 
            !===============================================================
 
        case (unsteady)
 
            ! Unsteady computation. For a consistent restart nOldLevels
            ! solutions must be read. All restart files are provided explicitly
            ! from python script.
            call setsolfilenames()
 
            ! Check if the files can be opened, exit if that fails.
            call checksolfilenames()
 
            ! Set nOldSolAvail to nSolsRead and check if a consistent
            ! restart can be made. If not, processor 0 prints a warning.
 
            noldsolavail = nsolsread
            if (myid == 0 .and. noldsolavail < noldlevels) then
 
                print "(a)", "#"
                print "(a)", "#                 Warning"
                print "(a)", "# Not enough data found for a consistent &
                     &time accurate restart."
                print "(a)", "# Order is reduced in the first time steps &
                     &until enough data is available again."
                print "(a)", "#"
 
            end if
 
            ! Set the logicals oldSolWritten, such that nothing is
            ! overwritten when solution files are dumped for this
            ! computation.
 
            ii = min(nsolsread, noldlevels - 1)
            do nn = 1, ii
                oldsolwritten(nn) = .true.
            end do
 
            !===============================================================
 
        case (timespectral)
 
            ! Time spectral computation. For a consistent restart
            ! nTimeIntervalsSpectral solutions must be read. First
            ! determine the the restart files.
            call setsolfilenames()
 
            ! Check if the files can be opened, exit if that fails.
            call checksolfilenames()
 
            ! Check whether or not the spectral solution must be copied
            ! or interpolated.
 
            if (nsolsread == 1) then
                copyspectral = .true.
            else if (nsolsread /= ntimeintervalsspectral) then
                interpolspectral = .true.
            end if
 
        end select
 
    end subroutine determinesolfilenames
 
    subroutine setsolfilenames
        !
        !       setSolFileNames allocates and set the solution files that
        !       will be read and loaded in the restart
        !
        use communication
        use inputio
        use utils, only: terminate
        use variablereading, only: solfiles, nsolsread
        implicit none
        !
        !      Local variables
        !
        integer :: ierr
        integer(kind=intType) :: nn
 
        ! The length of the array provided gives the number of nSolsRead.
        nsolsread = SIZE(restartfiles, 1)
 
        ! Allocate the memory for the file names and set them.
        allocate (solfiles(nsolsread), stat=ierr)
        if (ierr /= 0) &
            call terminate("determineSolFileNames", &
                           "Memory allocation failure for solFiles")
 
        do nn = 1, nsolsread
            solfiles(nn) = restartfiles(nn)
        end do
 
    end subroutine setsolfilenames
 
    subroutine checksolfilenames
        !
        !       checkSolFileNames will check if the provided restart files
        !       are readable on disk. If not readable return fail, if readable
        !       message will be printed to let the user know that restart file
        !       will be tried to read.
        !
        use communication
        use utils, only: terminate
        use variablereading, only: solfiles, nsolsread
        use commonformats, only: strings
        implicit none
        !
        !      Local variables
        !
        integer :: ierr
        character(len=maxStringLen) :: errorMessage
        integer(kind=intType) :: nn
 
        do nn = 1, nsolsread
            open (unit=21, file=solfiles(nn), status="old", iostat=ierr)
            if (ierr /= 0) then
                write (errormessage, *) "Restart file ", trim(solfiles(nn)), &
                    " could not be opened for reading"
                call terminate("checkSolFileNames", errormessage)
                exit
            end if
            close (unit=21)
 
            if (myid == 0) then
                print "(a)", "#"
                write (*, strings) "# Found restart file: ", trim(solfiles(nn))
                print "(a)", "#"
            end if
        end do
 
    end subroutine checksolfilenames
 
    subroutine initdepvarandhalos(halosRead)
        !
        !       InitDepvarAndHalos computes the dependent flow variables,
        !       like viscosities, and initializes the halo cells by applying
        !       the boundary conditions and exchanging the internal halo's.
        !       This is all done on the start level grid.
        !
        use blockpointers
        use flowvarrefstate
        use inputio
        use inputiteration
        use inputphysics
        use inputtimespectral
        use iteration
        use monitor
        use section
        use utils, only: setpointers
        use haloexchange, only: whalo2
        use turbutils, only: computeeddyviscosity
        use turbbcroutines, only: applyallturbbc
        use solverutils, only: gridvelocitiesfinelevel, gridvelocitiescoarselevels, &
                               normalvelocitiesalllevels, slipvelocitiesfinelevel, slipvelocitiescoarselevels
        use flowutils, only: computelamviscosity
        use bcroutines, only: applyallbc
        use residuals, only: residual
        implicit none
        !
        !      Subroutine arguments.
        !
        logical, intent(in) :: halosRead
        !
        !      Local variables.
        !
        integer(kind=intType) :: nn, mm
        real(kind=realtype) :: relaxbleedsor
 
        real(kind=realtype), dimension(nSections) :: t
 
        logical :: initBleeds
 
        ! Set the logical whether or not to initialize the prescribed
        ! data for the bleed regions. If the halos were read the bleeds
        ! have been initialized already and nothing needs to be done.
 
        if (halosread) then
            initbleeds = .false.
        else
            initbleeds = .true.
        end if
 
        ! Compute the face velocities and for viscous walls the slip
        ! velocities. This is done for all the mesh levels.
 
        currentlevel = 1
        groundlevel = 1
 
        do mm = 1, ntimeintervalsspectral
 
            ! Compute the time, which corresponds to this spectral solution.
            ! For steady and unsteady mode this is simply the restart time;
            ! for the spectral mode the periodic time must be taken into
            ! account, which can be different for every section.
 
            t = timeunsteadyrestart
 
            if (equationmode == timespectral) then
                do nn = 1, nsections
                    t(nn) = t(nn) + (mm - 1) * sections(nn)%timePeriod &
                            / real(ntimeintervalsspectral, realtype)
                end do
            end if
            call gridvelocitiesfinelevel(.false., t, mm)
            call gridvelocitiescoarselevels(mm)
            call normalvelocitiesalllevels(mm)
 
            call slipvelocitiesfinelevel(.false., t, mm)
            call slipvelocitiescoarselevels(mm)
 
        end do
 
        ! Loop over the number of spectral solutions and blocks
        ! to compute the laminar viscosity.
 
        do mm = 1, ntimeintervalsspectral
            do nn = 1, ndom
                call setpointers(nn, mgstartlevel, mm)
                call computelamviscosity(.false.)
            end do
        end do
 
        ! Exchange the solution on the multigrid start level.
        ! It is possible that the halo values are needed for the boundary
        ! conditions. Viscosities are not exchanged.
 
        call whalo2(mgstartlevel, 1_inttype, nw, .true., .true., &
                    .false.)
 
        ! Apply all flow boundary conditions to be sure that the halo's
        ! contain the correct values. These might be needed to compute
        ! the eddy-viscosity. Also the data for the outflow bleeds
        ! is determined.
 
        currentlevel = mgstartlevel
        groundlevel = mgstartlevel
 
        call applyallbc(.true.)
 
        ! Loop over the number of spectral solutions and blocks
        ! to compute the eddy viscosities.
 
        do mm = 1, ntimeintervalsspectral
            do nn = 1, ndom
 
                ! Set the pointers for this block.
 
                call setpointers(nn, mgstartlevel, mm)
 
                ! Compute the eddy viscosity for rans computations using
                ! an eddy viscosity model.
 
                call computeeddyviscosity(.false.)
 
                ! In case of a rans computation and no restart, initialize
                ! the turbulent variables a bit better for some turbulence
                ! models.
 
                ! if(equations == RANSEquations .and. (.not. restart)) then
                !
                !   select case (turbModel)
                !
                !     case (komegaWilcox, komegaModified, menterSST)
                !       call initKOmega(0_intType)
                !       call computeEddyViscosity
                !
                !   end select
                !
                ! endif
 
            end do
        end do
 
        ! Exchange the laminar and eddy viscosities.
 
        call whalo2(mgstartlevel, 1_inttype, 0_inttype, .false., &
                    .false., .true.)
 
        if (equations == ransequations) then
            call applyallturbbc(.true.)
        end if
        call applyallbc(.true.)
 
        ! Exchange the solution for the second time to be sure that all
        ! halo's are initialized correctly. As this is the initialization
        ! phase, this is not critical.
 
        call whalo2(mgstartlevel, 1_inttype, nw, .true., .true., &
                    .true.)
 
    end subroutine initdepvarandhalos
 
    subroutine initflowrestart
        !
        !       initFlowRestart loads restart information from the restart
        !       file into the state variables.
        !
        use constants
        use iomodule, only: iovar
        use variablereading, only: solfiles, interpolspectral, copyspectral, &
                                   halosread
        use utils, only: terminate
        implicit none
        !
        !      Local variables.
        !
        integer(kind=intType) :: ierr
 
        ! Initialize halosRead to .false. This will be overwritten
        ! if halo values are read during the restart.
 
        halosread = .false.
 
        ! Determine the number and names of the solution files and
        ! determine the contents of IOVar, the data structure to
        ! generalize the IO.
 
        call determinesolfilenames
        call setiovar
 
        ! Determine the format of the files and read them.
        ! Note that halosRead is possibly overwritten in the
        ! folloing select case statement below
 
        call readrestartfile()
 
        ! Copy or interpolate the spectral solution, if needed.
 
        if (copyspectral) call copyspectralsolution
        if (interpolspectral) call interpolatespectralsolution
 
        ! Release the memory of the file names and IOVar.
 
        deallocate (solfiles, iovar, stat=ierr)
        if (ierr /= 0) &
            call terminate("initFlow", &
                           "Deallocation failure for solFiles and IOVar")
 
        ! At the moment the pressure is stored at the location of the
        ! total energy. Copy the pressure to its own arrays and
        ! compute the total energy.
 
        call setpressureandcomputeenergy(halosread)
 
        ! Initialize the halo cells if a restart is performed and the
        ! entire flow field if this is not the case.
 
        call initializehalos(halosread)
 
        ! Initialize the dependent flow variables and the halo values.
        call initdepvarandhalos(halosread)
 
    end subroutine initflowrestart
 
    subroutine initflowfield
        !
        !       initFlowfield initializes the flow field to a uniform flow on
        !       the start level grid. Exception may be some turbulence
        !       variables, which are initialized a bit smarter.
        !
        use constants
        use communication, only: myid
        use inputiteration, only: ncycles, nsgstartup
        use inputphysics, only: equationmode
        use inputunsteady, only: ntimestepsfine
        use iteration, only: noldsolavail
        use monitor, only: ntimestepsrestart, timeunsteadyrestart
        use utils, only: allocconvarrays, alloctimearrays
        implicit none
        !
        !      Local variables.
        !
        integer(kind=intType) :: nn
 
        ! Initialize nTimeStepsRestart to 0 (no restart
        ! is performed) and allocate the memory for the arrays to store
        ! convergence history. This allocation is only to be done by
        ! processor 0. For an unsteady computation the entire convergence
        ! history is stored and the values after every time step is
        ! stored in a separate array.
 
        ntimestepsrestart = 0
        timeunsteadyrestart = zero
 
        nn = nsgstartup + ncycles
        if (equationmode == unsteady) nn = ntimestepsfine * nn
        if (myid == 0) call allocconvarrays(nn)
        if (equationmode == unsteady .and. myid == 0) &
            call alloctimearrays(ntimestepsfine)
 
        ! Set nOldSolAvail to 1, to indicate that an unsteady
        ! computation should be started with a lower order
        ! time integration scheme, because not enough states
        ! in the past are available.
 
        noldsolavail = 1
 
        ! Initialize the flow field to uniform flow.
 
        call setuniformflow
 
    end subroutine initflowfield
 
    subroutine initializehalos(halosRead)
        !
        !       initializeHalos sets the flow variables in the halo cells
        !       using a constant extrapolation. If the halos are read only the
        !       second halos are initialized, otherwise both.
        !
        use constants
        use blockpointers, only: ndom, w, p, rlv, ib, jb, kb, &
                                 il, ie, jl, je, kl, ke, rev
        use flowvarrefstate, only: viscous, nw, eddymodel, muinf
        use inputiteration, only: mgstartlevel
        use inputphysics, only: eddyvisinfratio
        use inputtimespectral, only: ntimeintervalsspectral
        use utils, only: setpointers
        implicit none
        !
        !      Subroutine arguments.
        !
        logical, intent(in) :: halosRead
        !
        !      Local variables.
        !
        integer(kind=intType) :: nn, mm, i, j, k, l
        integer(kind=intType) :: jj, kk
 
        ! Loop over the number of spectral solutions and blocks.
 
        spectralloop: do mm = 1, ntimeintervalsspectral
            domains: do nn = 1, ndom
 
                ! Set the pointers for this block.
 
                call setpointers(nn, mgstartlevel, mm)
 
                ! Determine the situation we are dealing with.
 
                testhalosread: if (halosread) then
 
                    ! The first layer of halo cells have been read. Initialize
                    ! the second layer from the first layer.
 
                    ! Halo cells in the i-direction.
 
                    do k = 0, kb
                        kk = max(1_inttype, min(k, ke))
                        do j = 0, jb
                            jj = max(1_inttype, min(j, je))
 
                            do l = 1, nw
                                w(0, j, k, l) = w(1, jj, kk, l)
                                w(ib, j, k, l) = w(ie, jj, kk, l)
                            end do
 
                            p(0, j, k) = p(1, jj, kk)
                            p(ib, j, k) = p(ie, jj, kk)
 
                        end do
                    end do
 
                    ! Halo cells in j-direction. Note that the i-halo's have
                    ! already been set.
 
                    do k = 0, kb
                        kk = max(1_inttype, min(k, ke))
                        do i = 1, ie
 
                            do l = 1, nw
                                w(i, 0, k, l) = w(i, 1, kk, l)
                                w(i, jb, k, l) = w(i, je, kk, l)
                            end do
 
                            p(i, 0, k) = p(i, 1, kk)
                            p(i, jb, k) = p(i, je, kk)
 
                        end do
                    end do
 
                    ! Halo cells in k-direction. Note that the halo's in both
                    ! i and j direction have already been set.
 
                    do j = 1, je
                        do i = 1, ie
 
                            do l = 1, nw
                                w(i, j, 0, l) = w(i, j, 1, l)
                                w(i, j, kb, l) = w(i, j, ke, l)
                            end do
 
                            p(i, j, 0) = p(i, j, 1)
                            p(i, j, kb) = p(i, j, ke)
 
                        end do
                    end do
 
                    else testhalosread
 
                    ! No halo cells have been read. Initialize both layers
                    ! using the internal value.
 
                    ! Halo cells in the i-direction.
 
                    do k = 0, kb
                        kk = max(2_inttype, min(k, kl))
                        do j = 0, jb
                            jj = max(2_inttype, min(j, jl))
 
                            do l = 1, nw
                                w(0, j, k, l) = w(2, jj, kk, l)
                                w(1, j, k, l) = w(2, jj, kk, l)
                                w(ie, j, k, l) = w(il, jj, kk, l)
                                w(ib, j, k, l) = w(il, jj, kk, l)
                            end do
 
                            p(0, j, k) = p(2, jj, kk)
                            p(1, j, k) = p(2, jj, kk)
                            p(ie, j, k) = p(il, jj, kk)
                            p(ib, j, k) = p(il, jj, kk)
 
                        end do
                    end do
 
                    ! Halo cells in j-direction. Note that the i-halo's have
                    ! already been set.
 
                    do k = 0, kb
                        kk = max(2_inttype, min(k, kl))
                        do i = 2, il
 
                            do l = 1, nw
                                w(i, 0, k, l) = w(i, 2, kk, l)
                                w(i, 1, k, l) = w(i, 2, kk, l)
                                w(i, je, k, l) = w(i, jl, kk, l)
                                w(i, jb, k, l) = w(i, jl, kk, l)
                            end do
 
                            p(i, 0, k) = p(i, 2, kk)
                            p(i, 1, k) = p(i, 2, kk)
                            p(i, je, k) = p(i, jl, kk)
                            p(i, jb, k) = p(i, jl, kk)
 
                        end do
                    end do
 
                    ! Halo cells in k-direction. Note that the halo's in both
                    ! i and j direction have already been set.
 
                    do j = 2, jl
                        do i = 2, il
 
                            do l = 1, nw
                                w(i, j, 0, l) = w(i, j, 2, l)
                                w(i, j, 1, l) = w(i, j, 2, l)
                                w(i, j, ke, l) = w(i, j, kl, l)
                                w(i, j, kb, l) = w(i, j, kl, l)
                            end do
 
                            p(i, j, 0) = p(i, j, 2)
                            p(i, j, 1) = p(i, j, 2)
                            p(i, j, ke) = p(i, j, kl)
                            p(i, j, kb) = p(i, j, kl)
 
                        end do
                    end do
 
                end if testhalosread
 
                ! Initialize the laminar and eddy viscosity, if appropriate,
                ! such that no uninitialized memory is present.
                ! As the viscosities are dependent variables their values
                ! are not read during a restart.
 
                if (viscous) rlv = muinf
                if (eddymodel) rev = eddyvisinfratio * muinf
 
            end do domains
        end do spectralloop
 
    end subroutine initializehalos
    subroutine interpolatespectralsolution
        !
        !       interpolateSpectralSolution uses a spectral interpolation to
        !       determine the initialization of the flow solution.
        !       The solution is interpolated from the solution read, which
        !       contains a different number of time instances and is stored in
        !       IOVar()%w. This variable can be found in IOModule.
        !
        use constants
        use blockpointers, only: w, il, jl, kl, ndom, sectionid
        use flowvarrefstate, only: nw
        use inputtimespectral, only: ntimeintervalsspectral
        use iomodule, only: iovar
        use section, only: sections, nsections
        use utils, only: setpointers, terminate, spectralinterpolcoef
        use variablereading, only: nsolsread
        implicit none
        !
        !      Local variables.
        !
        integer :: ierr
 
        integer(kind=intType) :: jj, nn, ll, sps, i, j, k, l
 
        real(kind=realtype) :: t
 
        real(kind=realtype), dimension(nSolsRead) :: alpscal
        real(kind=realtype), dimension(nSections, nSolsRead, 3, 3) :: alpmat
 
        ! Loop over the number of spectral solutions to be interpolated.
 
        spectralloop: do sps = 1, ntimeintervalsspectral
 
            ! Determine the ratio of the time of this solution and the
            ! periodic time.
 
            nn = sps - 1
            t = real(nn, realtype) / real(ntimeintervalsspectral, realtype)
 
            ! Determine the interpolation coefficients for both the scalar
            ! and the vector quantities.
 
            call spectralinterpolcoef(nsolsread, t, alpscal, alpmat)
 
            ! Loop over the local number of blocks.
 
            domains: do nn = 1, ndom
 
                ! Set the pointers for this block to the finest grid level.
 
                call setpointers(nn, 1_inttype, sps)
 
                ! Loop over the number of variables to be interpolated.
 
                varloop: do l = 1, nw
 
                    ! Check if this is a velocity variable.
 
                    veltest: if (l == ivx .or. l == ivy .or. l == ivz) then
 
                        ! Velocity variable. Set ll, which is the row in the
                        ! matrix coefficients.
 
                        select case (l)
                        case (ivx)
                            ll = 1
                        case (ivy)
                            ll = 2
                        case (ivz)
                            ll = 3
                        end select
 
                        ! Loop over the owned cells to interpolate the variable.
 
                        do k = 2, kl
                            do j = 2, jl
                                do i = 2, il
 
                                    ! Initialization to zero and loop over the number of
                                    ! spectral solutions used in the interpolation and
                                    ! update the variable accordingly. Note that for the
                                    ! vector variables the matrix coefficients must be
                                    ! used; these matrices can be different for the
                                    ! sections present in the grid.
 
                                    w(i, j, k, l) = zero
                                    do jj = 1, nsolsread
                                        w(i, j, k, l) = w(i, j, k, l) &
                                                        + alpmat(sectionid, jj, ll, 1) &
                                                        * iovar(nn, jj)%w(i, j, k, ivx) &
                                                        + alpmat(sectionid, jj, ll, 2) &
                                                        * iovar(nn, jj)%w(i, j, k, ivy) &
                                                        + alpmat(sectionid, jj, ll, 3) &
                                                        * iovar(nn, jj)%w(i, j, k, ivz)
                                    end do
 
                                end do
                            end do
                        end do
 
                        else veltest
 
                        ! Scalar variable.
                        ! Loop over the owned cells to interpolate the variable.
 
                        do k = 2, kl
                            do j = 2, jl
                                do i = 2, il
 
                                    ! Initialization to zero and loop over the number of
                                    ! spectral solutions used in the interpolation and
                                    ! update the variable accordingly.
 
                                    w(i, j, k, l) = zero
                                    do jj = 1, nsolsread
                                        w(i, j, k, l) = w(i, j, k, l) &
                                                        + alpscal(jj) * iovar(nn, jj)%w(i, j, k, l)
                                    end do
 
                                end do
                            end do
                        end do
 
                    end if veltest
 
                end do varloop
 
            end do domains
 
        end do spectralloop
 
        ! Release the memory of w of IOVar.
 
        do sps = 1, nsolsread
            do nn = 1, ndom
                deallocate (iovar(nn, sps)%w, stat=ierr)
                if (ierr /= 0) &
                    call terminate("interpolateSpectralSolution", &
                                   "Deallocation failure for w.")
            end do
        end do
 
    end subroutine interpolatespectralsolution
 
    subroutine releaseextramembcs
        !
        !       releaseExtraMemBCs releases the extra memory allocated in
        !       allocMemBcdata. This additional memory was allocated, such
        !       that alternative boundary condition treatments can be handled
        !       in setBCDataFineGrid.
        !
        use constants
        use blockpointers, only: flowdoms, ndom, bctype, bcdata, nbocos
        use inputtimespectral, only: ntimeintervalsspectral
        use utils, only: setpointers, terminate
        implicit none
        !
        !      Local variables.
        !
        integer :: ierr
 
        integer(kind=intType) :: mm, nn, sps, level, nLevels, ii
 
        ! Determine the number of multigrid levels.
 
        nlevels = ubound(flowdoms, 2)
 
        ! Loop over the number of multigrid level, spectral solutions
        ! and local blocks.
 
        levelloop: do level = 1, nlevels
            spectralloop: do sps = 1, ntimeintervalsspectral
                domainsloop: do nn = 1, ndom
 
                    ! Have the pointers in blockPointers point to the
                    ! current block to make everything more readable.
 
                    call setpointers(nn, level, sps)
 
                    ! Loop over the number of boundary subfaces for this block.
 
                    bocoloop: do mm = 1, nbocos
 
                        ! Determine the boundary condition we are having here.
 
                        inflowtype:select case(bctype(mm))
 
                        case (subsonicinflow)
 
                        ! Subsonic inflow. Determine the boundary condition
                        ! treatment and release the accordingly. Note that
                        ! the boundary condition treatment of the finest mesh
                        ! must be used, because this data is not available yet
                        ! on the coarse grid.
 
                        ii = &
                            flowdoms(nn, 1, sps)%BCData(mm)%subsonicInletTreatment
 
                        select case (ii)
 
                        case (totalconditions)
 
                            ! Total conditions used. Release the memory for
                            ! the density and velocities and nullify their
                            ! pointers.
 
                            deallocate (bcdata(mm)%rho, bcdata(mm)%velx, &
                                        bcdata(mm)%vely, bcdata(mm)%velz, &
                                        stat=ierr)
                            if (ierr /= 0) &
                                 call terminate("releaseExtraMemBCs", &
                                 "Deallocation failure for rho, &
                                 &velx, vely and velz")
 
                            nullify (bcdata(mm)%rho)
                            nullify (bcdata(mm)%velx)
                            nullify (bcdata(mm)%vely)
                            nullify (bcdata(mm)%velz)
                            !===================================================
 
                        case (massflow)
 
                            ! Full velocity vector and density prescribed at
                            ! inlet boundaries. Release the memory for the
                            ! total conditions and flow directions and nullify
                            ! their pointers.
 
                            deallocate (bcdata(mm)%ptInlet, &
                                        bcdata(mm)%ttInlet, &
                                        bcdata(mm)%htInlet, &
                                        bcdata(mm)%flowXdirInlet, &
                                        bcdata(mm)%flowYdirInlet, &
                                        bcdata(mm)%flowZdirInlet, stat=ierr)
                            if (ierr /= 0) &
                                 call terminate("releaseExtraMemBCs", &
                                 "Deallocation failure for the &
                                 &total conditions.")
 
                            nullify (bcdata(mm)%ptInlet)
                            nullify (bcdata(mm)%ttInlet)
                            nullify (bcdata(mm)%htInlet)
                            nullify (bcdata(mm)%flowXdirInlet)
                            nullify (bcdata(mm)%flowYdirInlet)
                            nullify (bcdata(mm)%flowZdirInlet)
 
                        end select
 
                        end select inflowtype
 
                    end do bocoloop
                end do domainsloop
            end do spectralloop
        end do levelloop
 
    end subroutine releaseextramembcs
 
    subroutine setiovar
        !
        !       setIOVar allocates the memory for the derived data type IOVar,
        !       which is simplifies the reading. If an interpolation must be
        !       performed for the time spectral method also the solution of
        !       this IO type is allocated. For all other cases the pointers of
        !       IOVar are set the the appropriate entries of flowDoms, with
        !       possible offset due to the usage of pointers.
        !
        use constants
        use block, only: flowdoms, ndom
        use flowvarrefstate, only: nw
        use inputphysics, only: equationmode
        use iomodule, only: iovar
        use utils, only: terminate
        use variablereading, only: interpolspectral, halosread, nsolsread
        implicit none
        !
        !      Local variables.
        !
        integer :: ierr
 
        integer(kind=intType) :: nn, mm, il, jl, kl
 
        ! Allocate the memory for IOVar.
 
        allocate (iovar(ndom, nsolsread), stat=ierr)
        if (ierr /= 0) &
            call terminate("setIOVar", &
                           "Memory allocation failure for solRead")
 
        ! Determine the equation mode we are solving and set the pointers
        ! of coorRead accordingly, or even allocate the memory, if needed.
 
        select case (equationmode)
 
        case (steady)
 
            ! Steady computation. Only one solution needs to be read.
            ! Loop over the number of blocks and set the pointers.
            ! No pointer offset is needed.
 
            do nn = 1, ndom
                iovar(nn, 1)%pointerOffset = 0
                iovar(nn, 1)%w => flowdoms(nn, 1, 1)%w(1:, 1:, 1:, :)
            end do
 
            !===============================================================
 
        case (unsteady)
 
            ! Unsteady computation. The first solution should be stored in
            ! w. For the others the pointers point to wOld. As the
            ! starting indices of wOld are 2, a pointer shift takes place
            ! here. I know this is a pain in the butt, but that's what
            ! we have to live with.
 
            do nn = 1, ndom
                iovar(nn, 1)%pointerOffset = 0
                iovar(nn, 1)%w => flowdoms(nn, 1, 1)%w(1:, 1:, 1:, :)
 
                do mm = 2, nsolsread
                    iovar(nn, mm)%pointerOffset = -1
                    iovar(nn, mm)%w => flowdoms(nn, 1, 1)%wOld(mm - 1, 2:, 2:, 2:, :)
                end do
            end do
 
            !===============================================================
 
        case (timespectral)
 
            ! Time spectral mode. A further check is required.
 
            testallocsolread: if (interpolspectral) then
 
                ! A restart is performed using a different number of time
                ! instances than the previous computation. Consequently the
                ! solutions will be interpolated later on. Hence some
                ! additional storage is required for the solutions and thus
                ! the w variables of IOVar are allocated. No halo data is
                ! needed here. No pointer offset either due to the explicit
                ! allocation.
 
                do nn = 1, ndom
                    il = flowdoms(nn, 1, 1)%il
                    jl = flowdoms(nn, 1, 1)%jl
                    kl = flowdoms(nn, 1, 1)%kl
 
                    do mm = 1, nsolsread
                        iovar(nn, mm)%pointerOffset = 0
 
                        allocate (iovar(nn, mm)%w(2:il, 2:jl, 2:kl, nw), stat=ierr)
                        if (ierr /= 0) &
                            call terminate("setIOVar", &
                                           "Memory allocation failure for w")
                    end do
 
                end do
 
                else testallocsolread
 
                ! A restart is made using either 1 solution or the correct
                ! number of solution instances. In both cases simply set
                ! the pointers of IOVar. No pointer offset is needed.
 
                do nn = 1, ndom
                    do mm = 1, nsolsread
                        iovar(nn, mm)%pointerOffset = 0
                        iovar(nn, mm)%w => flowdoms(nn, 1, mm)%w(1:, 1:, 1:, :)
                    end do
                end do
 
            end if testallocsolread
 
        end select
 
    end subroutine setiovar
 
    subroutine setpressureandcomputeenergy(halosRead)
        !
        !       Due to the usage of the variable IOVar, which generalizes the
        !       IO and leads to reuse of code, currently the pressure is
        !       stored at the position of rhoE. In this routine that data is
        !       copied to the pressure array and the total energy is computed.
        !       Note that this routine is only called when a restart is done.
        !
        use constants
        use blockpointers, only: ndom, p, w, il, jl, kl
        use flowvarrefstate, only: kpresent
        use inputtimespectral, only: ntimeintervalsspectral
        use utils, only: setpointers
        use flowutils, only: computeetotblock
        implicit none
        !
        !      Subroutine arguments.
        !
        logical, intent(in) :: halosRead
        !
        !      Local variables.
        !
        integer(kind=intType) :: sps, nn, nHalo
        integer(kind=intType) :: i, j, k
        integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, kBeg, kEnd
 
        ! Set the value of nHalo, depending whether or not the halo cells
        ! have been read from the restart file.
 
        nhalo = 0
        if (halosread) nhalo = 1
 
        ! Loop over the number of time instances and the local blocks.
        ! As this routine is only called when a restart is performed,
        ! the MG start level is the finest level.
 
        do sps = 1, ntimeintervalsspectral
            do nn = 1, ndom
 
                ! Set the pointers to the correct block. As this routine is
                ! only called when a restart is performed, the MG start level
                ! is the finest level.
 
                call setpointers(nn, 1_inttype, sps)
 
                ! Determine the range for which the pressure must be computed.
 
                ibeg = 2 - nhalo; jbeg = 2 - nhalo; kbeg = 2 - nhalo
                iend = il + nhalo; jend = jl + nhalo; kend = kl + nhalo
 
                ! Copy the pressure for the required cells.
 
                do k = kbeg, kend
                    do j = jbeg, jend
                        do i = ibeg, iend
                            p(i, j, k) = w(i, j, k, irhoe)
                        end do
                    end do
                end do
 
                ! Compute the total energy as well.
 
                call computeetotblock(ibeg, iend, jbeg, jend, kbeg, kend, kpresent)
            end do
        end do
 
    end subroutine setpressureandcomputeenergy
    subroutine setrestartfiles(fileName, i)
        !
        !          Populates the restartfiles
        !          The array is populated from Python using setRestartFiles
        !
        use constants
        use inputio, only: restartfiles
        implicit none
        !
        !      Subroutine argument.
        !
        character(len=*), intent(inout) :: fileName
        integer(kind=intType) :: i
 
        restartfiles(i) = filename
 
    end subroutine setrestartfiles
 
    subroutine setuniformflow
        !
        !       setUniformFlow set the flow variables of all local blocks on
        !       the start level to the uniform flow field.
        !
        use constants
        use blockpointers, only: w, dw, fw, flowdoms, ib, jb, kb, &
                                 rev, rlv, ndom, bcdata, p
        use communication
        use flowvarrefstate, only: eddymodel, viscous, muinf, nw, nwf, &
                                   pinfcorr, winf
        use inputiteration, only: mgstartlevel
        use inputphysics, only: equationmode, flowtype, eddyvisinfratio
        use inputtimespectral, only: ntimeintervalsspectral
        use utils, only: setpointers
        implicit none
        !
        !      Local variables.
        !
        integer :: ierr
 
        integer(kind=intType) :: nn, mm, i, j, k, l
 
        real(kind=realtype) :: tmp
 
        real(kind=realtype), dimension(3) :: dirloc, dirglob
 
        ! Loop over the number of spectral solutions and blocks.
        spectralloop: do mm = 1, ntimeintervalsspectral
            domains: do nn = 1, ndom
 
                ! Set the pointers for this block.
 
                call setpointers(nn, mgstartlevel, mm)
 
                ! Set the w-variables to the ones of the uniform flow field.
                do l = 1, nw
                    do k = 0, kb
                        do j = 0, jb
                            do i = 0, ib
                                w(i, j, k, l) = winf(l)
                                dw(i, j, k, l) = zero
                            end do
                        end do
                    end do
                end do
                !set this here for a reinitialize flow to eliminate possible NAN's
                do l = 1, nwf
                    do k = 0, kb
                        do j = 0, jb
                            do i = 0, ib
                                fw(i, j, k, l) = zero
                            end do
                        end do
                    end do
                end do
 
                ! Set the pressure.
 
                p = pinfcorr
 
                ! Initialize the laminar and eddy viscosity, if appropriate,
                ! such that no uninitialized memory is present.
 
                if (viscous) rlv = muinf
                if (eddymodel) rev = eddyvisinfratio * muinf
 
            end do domains
        end do spectralloop
 
        ! Correct for the time spectral method in combination with an
        ! internal flow computation the velocity direction.
        ! It is possible that the prescribed direction is different
        ! for every time instance.
 
        testcorrection: if (equationmode == timespectral .and. &
                            flowtype == internalflow) then
 
            ! Loop over the number of spectral solutions.
 
            spectralloopcorr: do mm = 1, ntimeintervalsspectral
 
                ! Initialize the local direction to zero. In the loop
                ! below this direction will be accumulated.
 
                dirloc = zero
 
                ! Loop over the number of blocks and determine the average
                ! velocity direction prescribed at the inlets.
 
                domainloop1: do nn = 1, ndom
 
                    ! Set the pointer for the BCData to make the code more
                    ! readable.
 
                    bcdata => flowdoms(nn, mgstartlevel, mm)%BCData
 
                    ! Loop over the number of boundary subfaces and update
                    ! the local prescribed velocity direction.
 
                    do l = 1, flowdoms(nn, mgstartlevel, mm)%nBocos
                        call velmagnanddirectionsubface(tmp, dirloc, bcdata, l)
                    end do
                end do domainloop1
 
                ! Determine the sum of dirLoc and create a unit vector
                ! for the global direction.
 
                call mpi_allreduce(dirloc, dirglob, 3, adflow_real, mpi_sum, &
                                   adflow_comm_world, ierr)
 
                tmp = one / max(eps, sqrt(dirglob(1)**2 &
                                          + dirglob(2)**2 &
                                          + dirglob(3)**2))
                dirglob(1) = tmp * dirglob(1)
                dirglob(2) = tmp * dirglob(2)
                dirglob(3) = tmp * dirglob(3)
 
                ! Loop again over the local domains and correct the
                ! velocity direction.
 
                domainsloop2: do nn = 1, ndom
 
                    ! Set the pointers for this block.
 
                    call setpointers(nn, mgstartlevel, mm)
 
                    do k = 0, kb
                        do j = 0, jb
                            do i = 0, ib
                                tmp = sqrt(w(i, j, k, ivx)**2 &
                                           + w(i, j, k, ivy)**2 &
                                           + w(i, j, k, ivz)**2)
                                w(i, j, k, ivx) = tmp * dirglob(1)
                                w(i, j, k, ivy) = tmp * dirglob(2)
                                w(i, j, k, ivz) = tmp * dirglob(3)
                            end do
                        end do
                    end do
 
                end do domainsloop2
            end do spectralloopcorr
 
        end if testcorrection
 
    end subroutine setuniformflow
 
    !=================================================================
 
    subroutine velmagnanddirectionsubface(vmag, dir, BCData, mm)
        !
        !       VelMagnAndDirectionSubface determines the maximum value
        !       of the magnitude of the velocity as well as the sum of the
        !       flow directions for the currently active subface.
        !
        use constants
        use block
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: mm
 
        real(kind=realtype), intent(out) :: vmag
        real(kind=realtype), dimension(3), intent(inout) :: dir
 
        type(bcdatatype), dimension(:), pointer :: bcdata
        !
        !      Local variables.
        !
        integer(kind=intType) :: i, j
        real(kind=realtype) :: vel
 
        ! Initialize vmag to -1.0.
 
        vmag = -one
 
        ! Check if the velocity is prescribed.
 
        if (associated(bcdata(mm)%velx) .and. &
            associated(bcdata(mm)%vely) .and. &
            associated(bcdata(mm)%velz)) then
 
            ! Loop over the owned faces of the subface. As the cell range
            ! may contain halo values, the nodal range is used.
 
            do j = (bcdata(mm)%jnBeg + 1), bcdata(mm)%jnEnd
                do i = (bcdata(mm)%inBeg + 1), bcdata(mm)%inEnd
 
                    ! Compute the magnitude of the velocity and compare it
                    ! with the current maximum. Store the maximum of the two.
 
                    vel = sqrt(bcdata(mm)%velx(i, j)**2 &
                               + bcdata(mm)%vely(i, j)**2 &
                               + bcdata(mm)%velz(i, j)**2)
                    vmag = max(vmag, vel)
 
                    ! Compute the unit vector of the velocity and add it to dir.
 
                    vel = one / max(eps, vel)
                    dir(1) = dir(1) + vel * bcdata(mm)%velx(i, j)
                    dir(2) = dir(2) + vel * bcdata(mm)%vely(i, j)
                    dir(3) = dir(3) + vel * bcdata(mm)%velz(i, j)
 
                end do
            end do
        end if
 
        ! Check if the velocity direction is prescribed.
 
        if (associated(bcdata(mm)%flowXdirInlet) .and. &
            associated(bcdata(mm)%flowYdirInlet) .and. &
            associated(bcdata(mm)%flowZdirInlet)) then
 
            ! Add the unit vectors to dir by looping over the owned
            ! faces of the subfaces. Again the nodal range must be
            ! used for this.
 
            do j = (bcdata(mm)%jnBeg + 1), bcdata(mm)%jnEnd
                do i = (bcdata(mm)%inBeg + 1), bcdata(mm)%inEnd
 
                    dir(1) = dir(1) + bcdata(mm)%flowXdirInlet(i, j)
                    dir(2) = dir(2) + bcdata(mm)%flowYdirInlet(i, j)
                    dir(3) = dir(3) + bcdata(mm)%flowZdirInlet(i, j)
 
                end do
            end do
 
        end if
 
    end subroutine velmagnanddirectionsubface
    subroutine timespectralcoef(coefSpectral, matrixCoefSpectral, &
                                diagMatCoefSpectral)
        !
        !       timeSpectralCoef computes the time integration coefficients
        !       for the time spectral method. As it is possible that sections
        !       have different periodic times these coefficients are
        !       determined for all the sections. For vector quantities, such
        !       as momentum, these coefficients can also be different due to
        !       rotation and the fact that only a part of the wheel is
        !       simulated.
        !
        use constants
        use flowvarrefstate, only: timeref
        use inputtimespectral, only: ntimeintervalsspectral, rotmatrixspectral
        use section, only: nsections, sections
        implicit none
        !
        !      Subroutine arguments.
        !
        real(kind=realtype), &
            dimension(nSections, nTimeIntervalsSpectral - 1), &
            intent(out) :: coefspectral
        real(kind=realtype), &
            dimension(nSections, nTimeIntervalsSpectral - 1, 3, 3), &
            intent(out) :: matrixcoefspectral
        real(kind=realtype), dimension(nSections, 3, 3), &
            intent(out) :: diagmatcoefspectral
        !
        !      Local variables.
        !
        integer(kind=intType) :: pp, nn, mm, ii, i, j, ntot
        real(kind=realtype) :: coef, dangle, angle, fact, slicesfact
 
        real(kind=realtype), dimension(3, 3) :: rotmat, tmp
 
        ! Loop over the number of sections.
 
        sectionloop: do mm = 1, nsections
 
            ! Initialize dAngle (smallest angle in the cotangent function)
            ! and coef, which is the multiplication factor in front of the
            ! cotangent/cosecant function. Coef is a combination of the 1/2
            ! and the 2*pi/timePeriod
 
            dangle = pi / real(ntimeintervalsspectral, realtype)
            coef = pi * timeref / sections(mm)%timePeriod
 
            ! Computation of the scalar coefficients.
 
            scalarloop: do nn = 1, (ntimeintervalsspectral - 1)
 
                angle = nn * dangle
 
                ! The coefficient for an odd and even number of time
                ! instances are different; the former is 1/sin, the
                ! latter cos/sin or 1/tan.
 
                coefspectral(mm, nn) = coef / sin(angle)
 
                if (mod(ntimeintervalsspectral, 2_inttype) == 0) &
                    coefspectral(mm, nn) = coefspectral(mm, nn) * cos(angle)
 
                ! Negate coef for the next spectral coefficient.
 
                coef = -coef
 
            end do scalarloop
 
            ! Initialize dAngle to the smallest angle in the cotangent
            ! or cosecant function. Now this angle is for the entire wheel,
            ! i.e. the number of slices must be taken into account.
 
            ntot = ntimeintervalsspectral * sections(mm)%nSlices
            dangle = pi / real(ntot, realtype)
 
            ! Initialize the rotation matrix to the unity matrix.
 
            rotmat(1, 1) = one; rotmat(1, 2) = zero; rotmat(1, 3) = zero
            rotmat(2, 1) = zero; rotmat(2, 2) = one; rotmat(2, 3) = zero
            rotmat(3, 1) = zero; rotmat(3, 2) = zero; rotmat(3, 3) = one
 
            ! Loop over the number of spectral coefficient to initialize the
            ! matrix coefficients; this is basically pp == 0 in the loop
            ! over the number slices. Use is made of the fact that the
            ! rotation matrix is the identity for pp == 0.
            ! coef changes sign at every time instance
 
            slicesfact = one / real(sections(mm)%nSlices, realtype)
            fact = one
 
            do nn = 1, (ntimeintervalsspectral - 1)
 
                ! Determine the scalar coefficient. This value depends now
                ! whether the total number of time instances in the wheel is
                ! odd or even.
 
                angle = nn * dangle
                coef = one / sin(angle)
 
                if (mod(ntot, 2_inttype) == 0) &
                    coef = coef * cos(angle)
 
                coef = coef * fact * slicesfact
 
                ! The first part of matrixCoefSpectral is a diagonal matrix,
                ! because this indicates the contribution of the current
                ! slice to the time derivative.
 
                matrixcoefspectral(mm, nn, 1, 1) = coef
                matrixcoefspectral(mm, nn, 1, 2) = zero
                matrixcoefspectral(mm, nn, 1, 3) = zero
 
                matrixcoefspectral(mm, nn, 2, 1) = zero
                matrixcoefspectral(mm, nn, 2, 2) = coef
                matrixcoefspectral(mm, nn, 2, 3) = zero
 
                matrixcoefspectral(mm, nn, 3, 1) = zero
                matrixcoefspectral(mm, nn, 3, 2) = zero
                matrixcoefspectral(mm, nn, 3, 3) = coef
 
                fact = -fact
 
            end do
 
            ! Initialize diagMatCoefSpectral to zero, because the
            ! starting index in the loop over the number of slices -1 is
            ! 1, i.e. the slice where the actual computation takes places
            ! does not contribute to diagMatCoefSpectral.
 
            do j = 1, 3
                do i = 1, 3
                    diagmatcoefspectral(mm, i, j) = zero
                end do
            end do
 
            ! Loop over the additional slices which complete an entire
            ! revolution. To be able to compute the coefficients a bit
            ! easier the loop runs from 1 to nSlices-1 and not from
            ! 2 to nSlices.
 
            slicesloop: do pp = 1, (sections(mm)%nSlices - 1)
 
                ! Compute the rotation matrix for this slice. This is the
                ! old one multiplied by the transformation matrix going from
                ! one slices to the next. Use tmp as temporary storage.
 
                do j = 1, 3
                    do i = 1, 3
                        tmp(i, j) = rotmatrixspectral(mm, i, 1) * rotmat(1, j) &
                                    + rotmatrixspectral(mm, i, 2) * rotmat(2, j) &
                                    + rotmatrixspectral(mm, i, 3) * rotmat(3, j)
                    end do
                end do
 
                rotmat = tmp
 
                slicesfact = one / real(sections(mm)%nSlices, realtype)
 
                ! Loop over the number of spectral coefficients and update
                ! matrixCoefSpectral. The multiplication with (-1)**nn
                ! takes place here too.
 
                ! Multiply also by the term (-1)**(pN+1)
 
                fact = one
                if (mod(pp * ntimeintervalsspectral, 2_inttype) /= 0) &
                    fact = -one
                slicesfact = fact * slicesfact
 
                fact = one
                ii = pp * ntimeintervalsspectral
                do nn = 1, (ntimeintervalsspectral - 1)
 
                    ! Compute the coefficient multiplying the rotation matrix.
                    ! Again make a distinction between an odd and an even
                    ! number of time instances for the entire wheel.
 
                    angle = (nn + ii) * dangle
                    coef = one / sin(angle)
 
                    if (mod(ntot, 2_inttype) == 0) &
                        coef = coef * cos(angle)
 
                    coef = coef * fact * slicesfact
 
                    ! Update matrixCoefSpectral.
 
                    do j = 1, 3
                        do i = 1, 3
                            matrixcoefspectral(mm, nn, i, j) = &
                                matrixcoefspectral(mm, nn, i, j) + coef * rotmat(i, j)
                        end do
                    end do
 
                    fact = -fact
 
                end do
 
                ! Update diagMatCoefSpectral. Also here the distinction
                ! between odd and even number of time instances.
 
                angle = ii * dangle
                coef = one / sin(angle)
 
                if (mod(ntot, 2_inttype) == 0) &
                    coef = coef * cos(angle)
 
                coef = coef * slicesfact
 
                do j = 1, 3
                    do i = 1, 3
                        diagmatcoefspectral(mm, i, j) = &
                            diagmatcoefspectral(mm, i, j) - coef * rotmat(i, j)
                    end do
                end do
 
            end do slicesloop
 
            ! The matrix coefficients must be multiplied by the leading
            ! coefficient, which depends on the actual periodic time.
 
            coef = pi * timeref / sections(mm)%timePeriod
 
            do j = 1, 3
                do i = 1, 3
                    diagmatcoefspectral(mm, i, j) = &
                        coef * diagmatcoefspectral(mm, i, j)
                end do
            end do
 
            do nn = 1, (ntimeintervalsspectral - 1)
 
                do j = 1, 3
                    do i = 1, 3
                        matrixcoefspectral(mm, nn, i, j) = &
                            coef * matrixcoefspectral(mm, nn, i, j)
                    end do
                end do
 
            end do
 
        end do sectionloop
 
    end subroutine timespectralcoef
 
    subroutine timespectralmatrices
        !
        !       timeSpectralMatrices computes the matrices for the time
        !       derivative of the time spectral method for all sections. For
        !       scalar quantities these matrices only differ if sections have
        !       different periodic times. For vector quantities, such as
        !       momentum, these matrices can be different depending on whether
        !       the section is rotating or not and the number of slices
        !       present.
        !
        use constants
        use inputphysics, only: equationmode
        use inputtimespectral, only: ntimeintervalsspectral, dscalar, dvector, &
                                     rotmatrixspectral
        use section, only: nsections
        use utils, only: terminate
        implicit none
        !
        !      Local variables.
        !
        integer :: ierr
 
        integer(kind=intType) :: nn, mm, ll, kk, ii
        integer(kind=intType) :: i, j
 
        real(kind=realtype), dimension(3, 3) :: tmpmat
 
        real(kind=realtype), dimension(:, :), allocatable :: coefspectral
        real(kind=realtype), dimension(:, :, :, :), allocatable :: &
            matrixcoefspectral
        real(kind=realtype), dimension(:, :, :), allocatable :: &
            diagmatcoefspectral
        !
        ! 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 the matrices as well as the help
        ! variables needed to construct these matrices.
 
        !added to allow second call in mdUpdateRoutines
        if (allocated(dscalar)) deallocate (dscalar)
        if (allocated(dvector)) deallocate (dvector)
 
        nn = ntimeintervalsspectral
        mm = 3 * nn
        kk = nn - 1
 
        allocate (dscalar(nsections, nn, nn), &
                  dvector(nsections, mm, mm), &
                  coefspectral(nsections, kk), &
                  matrixcoefspectral(nsections, kk, 3, 3), &
                  diagmatcoefspectral(nsections, 3, 3), stat=ierr)
        if (ierr /= 0) &
             call terminate("timeSpectralMatrices", &
             "Memory allocation failure for the matrices of &
             &the spectral time derivatives.")
 
        ! Determine the help variables needed to construct the
        ! actual matrices.
 
        call timespectralcoef(coefspectral, matrixcoefspectral, &
                              diagmatcoefspectral)
        !
        !       Determine the time derivative matrices for the sections.
        !
        ! Loop over the number of sections.
 
        sectionloop: do ii = 1, nsections
            !
            !         Matrix for scalar quantities.
            !
            ! Loop over the number of rows.
 
            do nn = 1, ntimeintervalsspectral
 
                ! Set the diagonal element to zero, i.e. there is no
                ! contribution to the own time derivative.
 
                dscalar(ii, nn, nn) = zero
 
                ! Loop over the rest of the columns.
 
                do mm = 1, (ntimeintervalsspectral - 1)
 
                    ! Determine the corresponding column index.
 
                    ll = nn + mm
                    if (ll > ntimeintervalsspectral) &
                        ll = ll - ntimeintervalsspectral
 
                    ! Store the corresponding coefficient in dscalar.
 
                    dscalar(ii, nn, ll) = coefspectral(ii, mm)
 
                end do
            end do
            !
            !         Matrices for vector quantities.
            !
            ! Loop over the number of time intervals; the number of rows
            ! is 3 times this number.
 
            rowloop: do nn = 1, ntimeintervalsspectral
 
                ! Initialize the diagonal block to diagMatCoefSpectral,
                ! the additional diagonal entry needed for the rotational
                ! periodicity.
 
                kk = 3 * (nn - 1)
                do j = 1, 3
                    do i = 1, 3
                        dvector(ii, kk + i, kk + j) = diagmatcoefspectral(ii, i, j)
                    end do
                end do
 
                ! Loop over the other time intervals, which contribute to
                ! the time derivative.
 
                columnloop: do mm = 1, (ntimeintervalsspectral - 1)
 
                    ! Determine the corresponding column index and check the
                    ! situation we are having here.
 
                    ll = nn + mm
                    if (ll > ntimeintervalsspectral) then
 
                        ! Index is outside the range and a shift must be applied.
 
                        ll = ll - ntimeintervalsspectral
 
                        ! The vector must be rotated. This effect is incorporated
                        ! directly in the matrix of time derivatives.
 
                        do j = 1, 3
                            do i = 1, 3
                                tmpmat(i, j) = matrixcoefspectral(ii, mm, i, 1) &
                                               * rotmatrixspectral(ii, 1, j) &
                                               + matrixcoefspectral(ii, mm, i, 2) &
                                               * rotmatrixspectral(ii, 2, j) &
                                               + matrixcoefspectral(ii, mm, i, 3) &
                                               * rotmatrixspectral(ii, 3, j)
                            end do
                        end do
 
                    else
 
                        ! Index is in the range. Copy the matrix coefficient
                        ! into tmpMat.
 
                        do j = 1, 3
                            do i = 1, 3
                                tmpmat(i, j) = matrixcoefspectral(ii, mm, i, j)
                            end do
                        end do
 
                    end if
 
                    ! Determine the offset for the column index and store
                    ! this submatrix in the correct place of dvector.
 
                    ll = 3 * (ll - 1)
                    do j = 1, 3
                        do i = 1, 3
                            dvector(ii, kk + i, ll + j) = tmpmat(i, j)
                        end do
                    end do
 
                end do columnloop
            end do rowloop
        end do sectionloop
 
        ! Release the memory of the help variables needed to construct
        ! the matrices of the time derivatives.
 
        deallocate (coefspectral, matrixcoefspectral, &
                    diagmatcoefspectral, stat=ierr)
        if (ierr /= 0) &
            call terminate("timeSpectralMatrices", &
                           "Deallocation failure for the help variables.")
 
    end subroutine timespectralmatrices
 
    subroutine readrestartfile()
        !
        !       readRestartFile reads the fine grid solution(s) from the
        !       restart file(s). If the restart file(s) do not correspond to
        !       the current mesh, the solution(s) are interpolated onto this
        !       mesh. It is also allowed to change boundary conditions, e.g.
        !       an alpha and/or Mach sweep is possible. Furthermore there is
        !       some support when starting from a different turbulence model,
        !       although this should be used with care.
        !
        use constants
        use cgnsgrid
        use su_cgns
        use variablereading ! Full import since we need basically everything
        use blockpointers, only: ibegor, jbegor, kbegor, il, jl, kl, ndom, &
                                 nbkglobal, nx, ny, nz
        use communication, only: adflow_comm_world, myid
        use inputphysics, only: equationmode
        use inputtimespectral, only: ntimeintervalsspectral
        use monitor, only: ntimestepsrestart, timeunsteadyrestart
        use utils, only: terminate, setpointers
        use sorting, only: bsearchstrings
        use commonformats, only: strings, stringint1
        implicit none
        !
        !      Local variables.
        !
        integer :: nzones, celldim, physdim, ierr, nsols
 
        integer(cgsize_t), dimension(9) :: sizes
        integer, dimension(9) :: rindsizes
        integer, dimension(nSolsRead) :: fileids
 
        integer(kind=intType) :: ii, jj, nn
        integer(kind=intType) :: ntypemismatch
        integer(kind=intType) :: nhimin, nhjmin, nhkmin
        integer(kind=intType) :: nhimax, nhjmax, nhkmax
 
        character(len=7) :: integerstring
        character(len=maxCGNSNameLen) :: cgnsname
        character(len=2*maxStringLen) :: errormessage
 
        ! Initialize halosRead to .true. This will be overwritten if
        ! there is at least one block present for which the halo data
        ! cannot be read.
 
        halosread = .true.
 
        ! Initialize nTypeMismatch to 0.
 
        ntypemismatch = 0
 
        ! Loop over the number of files to be read and open them.
 
        fileopenloop: do solid = 1, nsolsread
 
            ! Open the restart file for reading.
 
            call cg_open_f(solfiles(solid), mode_read, cgnsind, ierr)
            if (ierr /= all_ok) then
                write (errormessage, *) "File ", trim(solfiles(solid)), &
                    " could not be opened for reading"
                call terminate("readRestartFile", errormessage)
            end if
 
            fileids(solid) = cgnsind
 
            ! Determine the number of bases in the cgns file.
            ! This must be at least 1.
 
            call cg_nbases_f(cgnsind, cgnsbase, ierr)
            if (ierr /= all_ok) &
                call terminate("readRestartFile", &
                               "Something wrong when calling cg_nbases_f")
 
            if (cgnsbase < 1) then
                write (errormessage, *) "CGNS file ", trim(solfiles(solid)), &
                    " does not contain a base"
                call terminate("readRestartFile", errormessage)
            end if
 
            ! Only data from the first base is read. Information from
            ! higher bases is ignored.
 
            cgnsbase = 1
 
            ! Read the cell and physical dimensions as well as the name for
            ! this base.
 
            call cg_base_read_f(cgnsind, cgnsbase, cgnsname, celldim, &
                                physdim, ierr)
            if (ierr /= all_ok) &
                call terminate("readRestartFile", &
                               "Something wrong when calling cg_base_read_f")
 
            ! Check the cell and physical dimensions. Both must be 3 for
            ! this code to work.
 
            if (celldim /= 3 .or. physdim /= 3) then
                write (errormessage, stringint1) "Both the number of cell and physical dimensions should be 3, not ", &
                    celldim, " and ", physdim
                call terminate("readRestartFile", errormessage)
            end if
 
        end do fileopenloop
 
        ! Broadcast nTimeStepsRestart and timeUnsteadyRestart to all
        ! processors. These values are needed to perform a consistent
        ! unsteady restart.
 
        call mpi_bcast(ntimestepsrestart, 1, adflow_integer, 0, &
                       adflow_comm_world, ierr)
        call mpi_bcast(timeunsteadyrestart, 1, adflow_real, 0, &
                       adflow_comm_world, ierr)
 
        ! Get the scaling factors for density, pressure and velocity
        ! by reading the reference state.
 
        call scalefactors(fileids)
 
        ! Loop over the number of files to be read and read the solution.
 
        solloop: do solid = 1, nsolsread
 
            ! Store the file index a bit easier and set the base to 1.
 
            cgnsind = fileids(solid)
            cgnsbase = 1
 
            ! Determine the number of zones (blocks) in the restart file
            ! and check if this is identical to the number in the grid file.
 
            call cg_nzones_f(cgnsind, cgnsbase, nzones, ierr)
            if (ierr /= all_ok) &
                call terminate("readRestartFile", &
                               "Something wrong when calling cg_nzones_f")
 
            if (nzones /= cgnsndom) &
                 call terminate("readRestartFile", &
                 "Number of blocks in grid file and restart &
                 &file differ")
 
            ! Create a sorted version of the zone names of the restart file
            ! and store its corresponding zone numbers in zoneNumbers.
 
            call getsortedzonenumbers
 
            ! Loop over the number of blocks stored on this processor.
 
            domains: do nn = 1, ndom
 
                ! Set the pointers for this block. Make sure that the
                ! correct data is set.
 
                ii = min(solid, ntimeintervalsspectral)
                call setpointers(nn, 1_inttype, ii)
 
                ! Store the zone name of the original grid a bit easier.
 
                cgnsname = cgnsdoms(nbkglobal)%zoneName
 
                ! Search in the sorted zone names of the restart file for
                ! cgnsName. The name must be found; otherwise the restart
                ! is pointless. If found, the zone number is set accordingly.
 
                jj = bsearchstrings(cgnsname, zonenames)
                if (jj == 0) then
                    write (errormessage, *) "Zone name ", trim(cgnsname), &
                        " not found in restart file ", &
                        trim(solfiles(solid))
                    call terminate("readRestartFile", errormessage)
                else
                    jj = zonenumbers(jj)
                end if
 
                cgnszone = jj
 
                ! Determine the dimensions of the zone and check if these are
                ! identical to the dimensions of the block in the grid file.
 
                call cg_zone_read_f(cgnsind, cgnsbase, cgnszone, &
                                    cgnsname, sizes, ierr)
                if (ierr /= all_ok) &
                     call terminate("readRestartFile", &
                     "Something wrong when calling &
                     &cg_zone_read_f")
 
                if (cgnsdoms(nbkglobal)%il /= sizes(1) .or. &
                     cgnsdoms(nbkglobal)%jl /= sizes(2) .or. &
                     cgnsdoms(nbkglobal)%kl /= sizes(3)) &
                     call terminate("readRestartFile", &
                     "Corresponding zones in restart file and &
                     &grid file have different dimensions")
 
                ! Determine the number of flow solutions in this zone and
                ! check if there is a solution stored.
 
                call cg_nsols_f(cgnsind, cgnsbase, cgnszone, nsols, ierr)
                if (ierr /= all_ok) &
                    call terminate("readRestartFile", &
                                   "Something wrong when calling cg_nsols_f")
 
                if (nsols == 0) &
                    call terminate("readRestartFile", &
                                   "No solution present in restart file")
 
                ! Check for multiple solutions. A distinction is needed for
                ! overset cases because there will be an extra solution node
                ! for the nodal iblanks.
 
                if ((nsols > 1) .or. nsols > 2) &
                    call terminate("readRestartFile", &
                                   "Multiple solutions present in restart file")
 
                ! Determine the location of the solution variables. A loop is
                ! done over the solution nodes which is either 1 or 2. In the
                ! latter case, pick the node not named "Nodal Blanks".
 
                do cgnssol = 1, nsols
                    call cg_sol_info_f(cgnsind, cgnsbase, cgnszone, cgnssol, &
                                       cgnsname, location, ierr)
                    if (ierr /= all_ok) &
                         call terminate("readRestartFile", &
                         "Something wrong when calling &
                         &cg_sol_info_f")
 
                    if (trim(cgnsname) /= "Nodal Blanks") exit
                end do
 
                ! Determine the rind info.
 
                call cg_goto_f(cgnsind, cgnsbase, ierr, "Zone_t", &
                               cgnszone, "FlowSolution_t", cgnssol, "end")
                if (ierr /= all_ok) &
                    call terminate("readRestartFile", &
                                   "Something wrong when calling cg_goto_f")
 
                call cg_rind_read_f(rindsizes, ierr)
                if (ierr /= all_ok) &
                     call terminate("readRestartFile", &
                     "Something wrong when calling &
                     &cg_rind_read_f")
 
                ! Check if halo's are present. If not, set halosRead to .false.
                ! This only needs to be done if this is not an older state for
                ! an unsteady computation.
 
                if (solid == 1 .or. equationmode == timespectral) then
                    if (rindsizes(1) == 0 .or. rindsizes(2) == 0 .or. rindsizes(3) == 0 .or. &
                        rindsizes(4) == 0 .or. rindsizes(5) == 0 .or. rindsizes(6) == 0) &
                        halosread = .false.
                end if
 
                ! Initialize the number of halo cells to read to 0.
 
                nhimin = 0; nhjmin = 0; nhkmin = 0
                nhimax = 0; nhjmax = 0; nhkmax = 0
 
                ! Determine the range which must be read. A few things must be
                ! taken into account: - in iBegor, iEndor, etc. The nodal
                !                       range is stored. As in CGNS the cell
                !                       range start at 1, 1 must be subtracted
                !                       from the upper bound.
                !                     - the rind info must be taken into
                !                       account, because only the upper bound
                !                       is changed in cgns; the lower bound
                !                       remains 1.
                !                     - in case the solution is stored in the
                !                       vertices one extra variable in each
                !                       direction is read. An averaging will
                !                       take place to obtain cell centered
                !                       values.
                ! Also when vertex data is present, set halosRead to .false.,
                ! because it is not possible to determine the halos.
 
                if (location == cellcenter) then
 
                    ! Correct the number of halo cells to be read.
                    ! Only if this is not an older state in time for an
                    ! unsteady computation.
 
                    if (solid == 1 .or. equationmode == timespectral) then
                        if (rindsizes(1) > 0) nhimin = 1; if (rindsizes(2) > 0) nhimax = 1
                        if (rindsizes(3) > 0) nhjmin = 1; if (rindsizes(4) > 0) nhjmax = 1
                        if (rindsizes(5) > 0) nhkmin = 1; if (rindsizes(6) > 0) nhkmax = 1
                    end if
 
                    ! Set the cell range to be read from the CGNS file.
 
                    rangemin(1) = ibegor - nhimin
                    rangemin(2) = jbegor - nhjmin
                    rangemin(3) = kbegor - nhkmin
 
                    rangemax(1) = rangemin(1) + nx - 1 + nhimin + nhimax
                    rangemax(2) = rangemin(2) + ny - 1 + nhjmin + nhjmax
                    rangemax(3) = rangemin(3) + nz - 1 + nhkmin + nhkmax
 
                else if (location == vertex) then
 
                    ! Set the nodal range such that enough info is present
                    ! to average the nodal data to create the cell centered
                    ! data in the owned cells. No halo cells will be
                    ! initialized.
 
                    halosread = .false.
 
                    rangemin(1) = ibegor
                    rangemin(2) = jbegor
                    rangemin(3) = kbegor
 
                    rangemax(1) = rangemin(1) + nx
                    rangemax(2) = rangemin(2) + ny
                    rangemax(3) = rangemin(3) + nz
                else
                    call terminate("readRestartFile", &
                         "Only CellCenter or Vertex data allowed in &
                         &restart file")
                end if
 
                ! Allocate the memory for buffer, needed to store the variable
                ! to be read, and bufferVertex in case the solution is stored
                ! in the vertices.
 
                allocate (buffer(2 - nhimin:il + nhimax, &
                                 2 - nhjmin:jl + nhjmax, &
                                 2 - nhkmin:kl + nhkmax), stat=ierr)
                if (ierr /= 0) &
                    call terminate("readRestartFile", &
                                   "Memory allocation failure for buffer")
 
                if (location == vertex) then
                    allocate (buffervertex(1:il, 1:jl, 1:kl), stat=ierr)
                    if (ierr /= 0) &
                        call terminate("readRestartFile", &
                                       "Memory allocation failure for bufferVertex")
                end if
 
                ! Create a sorted version of the variable names and store the
                ! corresponding type in varTypes.
 
                call getsortedvarnumbers
 
                ! Read the density and the turbulence variables.
 
                call readdensity(ntypemismatch)
                call readturbvar(ntypemismatch)
 
                ! Read the other variables, depending on the situation.
 
                testprim: if (solid == 1 .or. equationmode == timespectral) then
 
                    ! Either the first solution or time spectral mode. Read
                    ! the primitive variables from the restart file.
 
                    call readxvelocity(ntypemismatch)
                    call readyvelocity(ntypemismatch)
                    call readzvelocity(ntypemismatch)
                    call readpressure(ntypemismatch)
 
                    else testprim
 
                    ! Old solution in unsteady mode. Read the conservative
                    ! variables.
 
                    call readxmomentum(ntypemismatch)
                    call readymomentum(ntypemismatch)
                    call readzmomentum(ntypemismatch)
                    call readenergy(ntypemismatch)
 
                end if testprim
 
                ! Release the memory of buffer, varNames and varTypes.
 
                deallocate (buffer, varnames, vartypes, stat=ierr)
                if (ierr /= 0) &
                     call terminate("readRestartFile", &
                     "Deallocation error for buffer, varNames &
                     &and varTypes.")
 
                ! In case bufferVertex is allocated, release it.
 
                if (location == vertex) then
                    deallocate (buffervertex, stat=ierr)
                    if (ierr /= 0) &
                        call terminate("readRestartFile", &
                                       "Deallocation error for bufferVertex")
                end if
 
            end do domains
 
            ! Release the memory of zoneNames and zoneNumbers.
 
            deallocate (zonenames, zonenumbers, stat=ierr)
            if (ierr /= 0) &
                 call terminate("readRestartFile", &
                 "Deallocation failure for zoneNames &
                 &and zoneNumbers.")
 
            ! Close the cgns solution file.
 
            call cg_close_f(cgnsind, ierr)
            if (ierr /= all_ok) &
                call terminate("readRestartFile", &
                               "Something wrong when calling cg_close_f")
 
        end do solloop
 
        ! Determine the global sum of nTypeMismatch; the result only
        ! needs to be known on processor 0. Use ii as the global buffer
        ! to store the result. If a type mismatch occured,
        ! print a warning.
 
        call mpi_reduce(ntypemismatch, ii, 1, adflow_integer, &
                        mpi_sum, 0, adflow_comm_world, ierr)
        if (myid == 0 .and. ii > 0) then
 
            write (integerstring, "(i6)") ii
            integerstring = adjustl(integerstring)
 
            print "(a)", "#"
            print "(a)", "#                      Warning"
            print strings, "# ", trim(integerstring), " type mismatches occured when reading the solution of the blocks"
            print "(a)", "#"
        end if
 
    end subroutine readrestartfile
 
    subroutine getsortedzonenumbers
        !
        !       getSortedZoneNumbers reads the names of the zones of the
        !       cgns file given by cgnsInd and cgnsBase. Afterwards the
        !       zonenames are sorted in increasing order, such that a binary
        !       search algorithm can be employed. The original zone numbers
        !       are stored in zoneNumbers.
        !       If the zone contains a link to a zone containing the
        !       coordinates the name of the linked zone is taken.
        !
 
        use constants
        use cgnsgrid, only: cgnsndom
        use su_cgns
        use variablereading, only: zonenames, zonenumbers, cgnsind, cgnsbase
        use sorting, only: qsortstrings, bsearchstrings
        use utils, only: terminate
        use commonformats, only: strings
        implicit none
        !
        !      Local variables.
        !
        integer :: ierr
        integer :: zone, zonetype, ncoords, pathlength
        integer :: pos
 
        integer(kind=cgsize_t), dimension(9) :: sizesblock
 
        integer(kind=intType) :: nn, ii
 
        character(len=maxStringLen) :: errormessage, linkpath
        character(len=maxCGNSNameLen), dimension(cgnsNdom) :: tmpnames
 
        logical :: namefound
 
        character(len=7) :: int1string, int2string
 
        ! Allocate the memory for zoneNames and zoneNumbers.
 
        allocate (zonenames(cgnsndom), zonenumbers(cgnsndom), stat=ierr)
        if (ierr /= 0) &
             call terminate("getSortedZoneNumbers", &
             "Memory allocation failure for zoneNames &
             &and zoneNumbers")
 
        ! Loop over the number of zones in the file.
 
        cgnsdomains: do nn = 1, cgnsndom
 
            ! Initialize nameFound to .false.
 
            namefound = .false.
 
            ! Check if the zone is structured.
 
            zone = nn
            call cg_zone_type_f(cgnsind, cgnsbase, zone, zonetype, ierr)
            if (ierr /= all_ok) &
                call terminate("getSortedZoneNumbers", &
                               "Something wrong when calling cg_zone_type_f")
 
            if (zonetype /= structured) then
 
                write (int1string, "(i7)") cgnsbase
                int1string = adjustl(int1string)
                write (int2string, "(i7)") zone
                int2string = adjustl(int2string)
 
                write (errormessage, strings) "Base ", trim(int1string), ": Zone ", trim(int2string), &
                    " of the cgns restart file is not structured"
                call terminate("getSortedZoneNumbers", errormessage)
 
            end if
 
            ! Determine the number of grid coordinates of this zone.
 
            call cg_ncoords_f(cgnsind, cgnsbase, zone, ncoords, ierr)
            if (ierr /= all_ok) &
                call terminate("getSortedZoneNumbers", &
                               "Something wrong when calling cg_ncoords_f")
 
            ! If ncoords == 3, there are coordinates present. Then check if
            ! it is a link. If so, take the zone name of the link.
 
            if (ncoords == 3) then
 
                ! Go to the coordinates node.
 
                call cg_goto_f(cgnsind, cgnsbase, ierr, "Zone_t", zone, &
                               "GridCoordinates_t", 1, "end")
                if (ierr /= all_ok) &
                    call terminate("getSortedZoneNumbers", &
                                   "Something wrong when calling cg_goto_f")
 
                ! Check if this node is a link.
 
                call cg_is_link_f(pathlength, ierr)
                if (ierr /= all_ok) &
                    call terminate("getSortedZoneNumbers", &
                                   "Something wrong when calling cg_is_link_f")
 
                if (pathlength > 0) then
 
                    ! Determine the name of the linkPath.
 
                    call cg_link_read_f(errormessage, linkpath, ierr)
                    if (ierr /= all_ok) &
                         call terminate("getSortedZoneNumbers", &
                         "Something wrong when calling &
                         &cg_link_read_f")
 
                    ! Find the zone name.
                    ! Find, starting from the back, the forward slash.
 
                    pos = index(linkpath, "/", .true.)
                    if (pos > 0) then
                        linkpath = linkpath(:pos - 1)
 
                        ! Find the next forward slash from the back and
                        ! remove the leading part from the path name.
 
                        pos = index(linkpath, "/", .true.)
                        if (pos > 0) linkpath = linkpath(pos + 1:)
                    end if
 
                    ! Create the zone name and set nameFound to .true..
 
                    linkpath = adjustl(linkpath)
                    zonenames(nn) = trim(linkpath)
                    namefound = .true.
 
                end if
            end if
 
            ! If no name was yet found set it to the name of the
            ! current zone.
 
            if (.not. namefound) then
                call cg_zone_read_f(cgnsind, cgnsbase, zone, &
                                    zonenames(nn), sizesblock, ierr)
                if (ierr /= all_ok) &
                     call terminate("getSortedZoneNumbers", &
                     "Something wrong when calling &
                     &cg_zone_read_f")
            end if
 
        end do cgnsdomains
 
        ! Set tmpNames to zoneNames and sort the latter
        ! in increasing order.
 
        do nn = 1, cgnsndom
            tmpnames(nn) = zonenames(nn)
        end do
 
        ! Sort zoneNames in increasing order.
 
        call qsortstrings(zonenames, cgnsndom)
 
        ! Initialize zoneNumbers to -1. This serves as a check during
        ! the search.
 
        do nn = 1, cgnsndom
            zonenumbers(nn) = -1
        end do
 
        ! Find the original zone numbers for the sorted zone names.
 
        do nn = 1, cgnsndom
            ii = bsearchstrings(tmpnames(nn), zonenames)
 
            ! Check if the zone number is not already taken. If this is the
            ! case, this means that two identical zone names are present.
 
            if (zonenumbers(ii) /= -1) &
                 call terminate("getSortedZoneNumbers", &
                 "Error occurs only when two identical zone &
                 &names are present")
 
            ! And set the zone number.
 
            zonenumbers(ii) = nn
        end do
 
    end subroutine getsortedzonenumbers
 
    subroutine getsortedvarnumbers
        !
        !       getSortedVarNumbers reads the names of variables stored in
        !       the given solution node of the cgns file, indicated by
        !       cgnsInd, cgnsBase and cgnsZone. Afterwards the variable
        !       names are sorted in increasing order, such that they can be
        !       used in a binary search. Their original variable number and
        !       type is stored.
        !
        use constants
        use su_cgns
        use variablereading, only: vartypes, varnames, cgnsbase, cgnsind, &
                                   cgnszone, cgnssol, nvar
        use sorting, only: qsortstrings, bsearchstrings
        use utils, only: terminate
        implicit none
        !
        !      Local variables.
        !
        integer :: i, ierr
        integer, dimension(:), allocatable :: tmptypes
 
        integer(kind=intType) :: nn, ii
 
        integer(kind=intType), dimension(:), allocatable :: varnumbers
 
        character(len=maxCGNSNameLen), allocatable, dimension(:) :: &
            tmpnames
 
        ! Determine the number of solution variables stored.
 
        call cg_nfields_f(cgnsind, cgnsbase, cgnszone, cgnssol, &
                          nvar, ierr)
        if (ierr /= all_ok) &
            call terminate("getSortedVarNumbers", &
                           "Something wrong when calling cg_nfield_f")
 
        ! Allocate the memory for varnames, vartypes and varnumber
 
        allocate (varnames(nvar), vartypes(nvar), varnumbers(nvar), &
                  stat=ierr)
        if (ierr /= 0) &
            call terminate("getSortedVarNumbers", &
                           "Memory allocation failure for varNames, etc.")
 
        ! Loop over the number of variables and store their names and
        ! types.
 
        do i = 1, nvar
            call cg_field_info_f(cgnsind, cgnsbase, cgnszone, cgnssol, &
                                 i, vartypes(i), varnames(i), ierr)
            if (ierr /= 0) &
                call terminate("getSortedVarNumbers", &
                               "Something wrong when calling cg_field_info_f")
        end do
 
        ! Allocate the memory for tmpTypes and tmpNames and initialize
        ! their values.
 
        allocate (tmptypes(nvar), tmpnames(nvar), stat=ierr)
        if (ierr /= 0) &
            call terminate("getSortedVarNumbers", &
                           "Memory allocation failure for tmp variables")
 
        do i = 1, nvar
            tmptypes(i) = vartypes(i)
            tmpnames(i) = varnames(i)
        end do
 
        ! Sort varNames in increasing order.
 
        nn = nvar
        call qsortstrings(varnames, nn)
 
        ! Initialize varNumbers to -1. This serves as a check during
        ! the search.
 
        do i = 1, nvar
            varnumbers(i) = -1
        end do
 
        ! Find the original types and numbers for the just sorted
        ! variable names.
 
        do i = 1, nvar
            ii = bsearchstrings(tmpnames(i), varnames)
 
            ! Check if the variable number is not already taken. If this is
            ! the case, this means that two identical var names are present.
 
            if (varnumbers(ii) /= -1) &
                 call terminate("getSortedVarNumbers", &
                 "Error occurs only when two identical &
                 &variable names are present")
 
            ! And set the variable number and type.
 
            varnumbers(ii) = i
            vartypes(ii) = tmptypes(i)
        end do
 
        ! Release the memory of varNumbers, tmpNames and tmpTypes.
 
        deallocate (varnumbers, tmptypes, tmpnames, stat=ierr)
        if (ierr /= 0) &
            call terminate("getSortedVarNumbers", &
                           "Deallocation error for tmp variables")
 
    end subroutine getsortedvarnumbers
#endif
end module initializeflow