solverUtils.F90

Go to the documentation of this file.

module solverutils
contains
#ifndef USE_TAPENADE
    subroutine timestep(onlyRadii)
        !
        ! Shell function to call timeStep_block on all blocks
        !
        use constants
        use blockpointers, only: ndom
        use inputtimespectral, only: ntimeintervalsspectral
        use iteration, only: currentlevel
        use utils, only: setpointers
        implicit none
        !
        !      Subroutine argument.
        !
        logical, intent(in) :: onlyRadii
        !
        !      Local variables.
        !
        integer(kind=intType) :: sps, nn
 
        ! Loop over the number of spectral solutions.
 
        spectralloop: do sps = 1, ntimeintervalsspectral
 
            ! Loop over the number of blocks.
 
            domains: do nn = 1, ndom
 
                ! Set the pointers for this block.
 
                call setpointers(nn, currentlevel, sps)
 
                call timestep_block(onlyradii)
 
            end do domains
 
        end do spectralloop
 
    end subroutine timestep
#endif
    subroutine timestep_block(onlyRadii)
        !
        !       timeStep computes the time step, or more precisely the time
        !       step divided by the volume per unit CFL, in the owned cells.
        !       However, for the artificial dissipation schemes, the spectral
        !       radIi in the halo's are needed. Therefore the loop is taken
        !       over the the first level of halo cells. The spectral radIi are
        !       stored and possibly modified for high aspect ratio cells.
        !
        use constants
        use blockpointers, only: ie, je, ke, il, jl, kl, w, p, rlv, rev, &
                                 radi, radj, radk, si, sj, sk, sfacei, sfacej, sfacek, dtl, gamma, vol, &
                                 addgridvelocities, sectionid
        use flowvarrefstate, only: timeref, eddymodel, gammainf, pinfcorr, &
                                   viscous, rhoinf
        use inputdiscretization, only: adis, dirscaling, radiineededcoarse, &
                                       radiineededfine, precond, acousticscalefactor
        use inputphysics, only: equationmode
        use iteration, only: groundlevel, currentlevel
        use section, only: sections
        use inputtimespectral, only: ntimeintervalsspectral
        use utils, only: terminate
        implicit none
        !
        !      Subroutine argument.
        !
        logical, intent(in) :: onlyRadii
        !
        !      Local parameters.
        !
        real(kind=realtype), parameter :: b = 2.0_realtype
        !
        !      Local variables.
        !
        integer(kind=intType) :: i, j, k, ii
 
        real(kind=realtype) :: plim, rlim, clim2
        real(kind=realtype) :: uux, uuy, uuz, cc2, qsi, qsj, qsk, sx, sy, sz, rmu
        real(kind=realtype) :: ri, rj, rk, rij, rjk, rki
        real(kind=realtype) :: vsi, vsj, vsk, rfl, dpi, dpj, dpk
        real(kind=realtype) :: sface, tmp
 
        logical :: radiiNeeded, doScaling
 
        ! Determine whether or not the spectral radii are needed for the
        ! flux computation.
 
        radiineeded = radiineededcoarse
        if (currentlevel <= groundlevel) radiineeded = radiineededfine
 
        ! Return immediately if only the spectral radii must be computed
        ! and these are not needed for the flux computation.
 
        if (onlyradii .and. (.not. radiineeded)) return
 
        ! Set the value of plim. To be fully consistent this must have
        ! the dimension of a pressure. Therefore a fraction of pInfCorr
        ! is used. Idem for rlim; compute clim2 as well.
 
        plim = 0.001_realtype * pinfcorr
        rlim = 0.001_realtype * rhoinf
        clim2 = 0.000001_realtype * gammainf * pinfcorr / rhoinf
 
        doscaling = (dirscaling .and. currentlevel <= groundlevel)
 
        ! Initialize sFace to zero. This value will be used if the
        ! block is not moving.
 
        sface = zero
        !
        !           Inviscid contribution, depending on the preconditioner.
        !           Compute the cell centered values of the spectral radii.
        !
        select case (precond)
 
        case (noprecond)
 
            ! No preconditioner. Simply the standard spectral radius.
            ! Loop over the cells, including the first level halo.
 
#ifdef TAPENADE_REVERSE
            !$AD II-LOOP
            do ii = 0, ie * je * ke - 1
                i = mod(ii, ie) + 1
                j = mod(ii / ie, je) + 1
                k = ii / (ie * je) + 1
#else
                do k = 1, ke
                    do j = 1, je
                        do i = 1, ie
#endif
                            ! Compute the velocities and speed of sound squared.
 
                            uux = w(i, j, k, ivx)
                            uuy = w(i, j, k, ivy)
                            uuz = w(i, j, k, ivz)
                            cc2 = gamma(i, j, k) * p(i, j, k) / w(i, j, k, irho)
                            cc2 = max(cc2, clim2)
 
                            ! Set the dot product of the grid velocity and the
                            ! normal in i-direction for a moving face. To avoid
                            ! a number of multiplications by 0.5 simply the sum
                            ! is taken.
 
                            if (addgridvelocities) &
                                sface = sfacei(i - 1, j, k) + sfacei(i, j, k)
 
                            ! Spectral radius in i-direction.
 
                            sx = si(i - 1, j, k, 1) + si(i, j, k, 1)
                            sy = si(i - 1, j, k, 2) + si(i, j, k, 2)
                            sz = si(i - 1, j, k, 3) + si(i, j, k, 3)
 
                            qsi = uux * sx + uuy * sy + uuz * sz - sface
 
                            ri = half * (abs(qsi) &
                                         + acousticscalefactor * sqrt(cc2 * (sx**2 + sy**2 + sz**2)))
 
                            ! The grid velocity in j-direction.
 
                            if (addgridvelocities) &
                                sface = sfacej(i, j - 1, k) + sfacej(i, j, k)
 
                            ! Spectral radius in j-direction.
 
                            sx = sj(i, j - 1, k, 1) + sj(i, j, k, 1)
                            sy = sj(i, j - 1, k, 2) + sj(i, j, k, 2)
                            sz = sj(i, j - 1, k, 3) + sj(i, j, k, 3)
 
                            qsj = uux * sx + uuy * sy + uuz * sz - sface
 
                            rj = half * (abs(qsj) &
                                         + acousticscalefactor * sqrt(cc2 * (sx**2 + sy**2 + sz**2)))
 
                            ! The grid velocity in k-direction.
 
                            if (addgridvelocities) &
                                sface = sfacek(i, j, k - 1) + sfacek(i, j, k)
 
                            ! Spectral radius in k-direction.
 
                            sx = sk(i, j, k - 1, 1) + sk(i, j, k, 1)
                            sy = sk(i, j, k - 1, 2) + sk(i, j, k, 2)
                            sz = sk(i, j, k - 1, 3) + sk(i, j, k, 3)
 
                            qsk = uux * sx + uuy * sy + uuz * sz - sface
 
                            rk = half * (abs(qsk) &
                                         + acousticscalefactor * sqrt(cc2 * (sx**2 + sy**2 + sz**2)))
 
                            ! Compute the inviscid contribution to the time step.
 
                            if (.not. onlyradii) dtl(i, j, k) = ri + rj + rk
 
                            !
                            !           Adapt the spectral radii if directional scaling must be
                            !           applied.
                            !
                            if (doscaling) then
 
                                ! Avoid division by zero by clipping radi, radJ and
                                ! radK.
 
                                ri = max(ri, eps)
                                rj = max(rj, eps)
                                rk = max(rk, eps)
 
                                ! Compute the scaling in the three coordinate
                                ! directions.
 
                                rij = (ri / rj)**adis
                                rjk = (rj / rk)**adis
                                rki = (rk / ri)**adis
 
                                ! Create the scaled versions of the aspect ratios.
                                ! Note that the multiplication is done with radi, radJ
                                ! and radK, such that the influence of the clipping
                                ! is negligible.
 
                                radi(i, j, k) = ri * (one + one / rij + rki)
                                radj(i, j, k) = rj * (one + one / rjk + rij)
                                radk(i, j, k) = rk * (one + one / rki + rjk)
                            else
                                radi(i, j, k) = ri
                                radj(i, j, k) = rj
                                radk(i, j, k) = rk
                            end if
#ifdef TAPENADE_REVERSE
                        end do
#else
                    end do
                end do
            end do
#endif
 
        case (turkel)
            call terminate("timeStep", "Turkel preconditioner not implemented yet")
 
        case (choimerkle)
            call terminate("timeStep", &
                           "choi merkle preconditioner not implemented yet")
        end select
 
        ! The rest of this file can be skipped if only the spectral
        ! radii need to be computed.
        testradiionly: if (.not. onlyradii) then
 
            ! The viscous contribution, if needed.
 
            viscousterm: if (viscous) then
 
                ! Loop over the owned cell centers.
 
                do k = 2, kl
                    do j = 2, jl
                        do i = 2, il
 
                            ! Compute the effective viscosity coefficient. The
                            ! factor 0.5 is a combination of two things. In the
                            ! standard central discretization of a second
                            ! derivative there is a factor 2 multiplying the
                            ! central node. However in the code below not the
                            ! average but the sum of the left and the right face
                            ! is taken and squared. This leads to a factor 4.
                            ! Combining both effects leads to 0.5. Furthermore,
                            ! it is divided by the volume and density to obtain
                            ! the correct dimensions and multiplied by the
                            ! non-dimensional factor factVis.
 
                            rmu = rlv(i, j, k)
                            if (eddymodel) rmu = rmu + rev(i, j, k)
                            rmu = half * rmu / (w(i, j, k, irho) * vol(i, j, k))
 
                            ! Add the viscous contribution in i-direction to the
                            ! (inverse) of the time step.
 
                            sx = si(i, j, k, 1) + si(i - 1, j, k, 1)
                            sy = si(i, j, k, 2) + si(i - 1, j, k, 2)
                            sz = si(i, j, k, 3) + si(i - 1, j, k, 3)
 
                            vsi = rmu * (sx * sx + sy * sy + sz * sz)
                            dtl(i, j, k) = dtl(i, j, k) + vsi
 
                            ! Add the viscous contribution in j-direction to the
                            ! (inverse) of the time step.
 
                            sx = sj(i, j, k, 1) + sj(i, j - 1, k, 1)
                            sy = sj(i, j, k, 2) + sj(i, j - 1, k, 2)
                            sz = sj(i, j, k, 3) + sj(i, j - 1, k, 3)
 
                            vsj = rmu * (sx * sx + sy * sy + sz * sz)
                            dtl(i, j, k) = dtl(i, j, k) + vsj
 
                            ! Add the viscous contribution in k-direction to the
                            ! (inverse) of the time step.
 
                            sx = sk(i, j, k, 1) + sk(i, j, k - 1, 1)
                            sy = sk(i, j, k, 2) + sk(i, j, k - 1, 2)
                            sz = sk(i, j, k, 3) + sk(i, j, k - 1, 3)
 
                            vsk = rmu * (sx * sx + sy * sy + sz * sz)
                            dtl(i, j, k) = dtl(i, j, k) + vsk
 
                        end do
                    end do
                end do
 
            end if viscousterm
 
            ! For the spectral mode an additional term term must be
            ! taken into account, which corresponds to the contribution
            ! of the highest frequency.
 
            if (equationmode == timespectral) then
 
                tmp = ntimeintervalsspectral * pi * timeref &
                      / sections(sectionid)%timePeriod
 
                ! Loop over the owned cell centers and add the term.
 
                do k = 2, kl
                    do j = 2, jl
                        do i = 2, il
                            dtl(i, j, k) = dtl(i, j, k) + tmp * vol(i, j, k)
                        end do
                    end do
                end do
 
            end if
 
            ! Currently the inverse of dt/vol is stored in dtl. Invert
            ! this value such that the time step per unit cfl number is
            ! stored and correct in cases of high gradients.
 
            do k = 2, kl
                do j = 2, jl
                    do i = 2, il
                        dpi = abs(p(i + 1, j, k) - two * p(i, j, k) + p(i - 1, j, k)) &
                              / (p(i + 1, j, k) + two * p(i, j, k) + p(i - 1, j, k) + plim)
                        dpj = abs(p(i, j + 1, k) - two * p(i, j, k) + p(i, j - 1, k)) &
                              / (p(i, j + 1, k) + two * p(i, j, k) + p(i, j - 1, k) + plim)
                        dpk = abs(p(i, j, k + 1) - two * p(i, j, k) + p(i, j, k - 1)) &
                              / (p(i, j, k + 1) + two * p(i, j, k) + p(i, j, k - 1) + plim)
                        rfl = one / (one + b * (dpi + dpj + dpk))
 
                        dtl(i, j, k) = rfl / dtl(i, j, k)
                    end do
                end do
            end do
 
        end if testradiionly
 
    end subroutine timestep_block
 
#ifndef USE_TAPENADE
    subroutine gridvelocitiesfinelevel(useOldCoor, t, sps)
        !
        ! Shell function to call gridVelocitiesFineLevel on all blocks
        !
        use blockpointers
        use constants
        use inputtimespectral
        use iteration
        use utils, only: setpointers
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: sps
        logical, intent(in) :: useOldCoor
        real(kind=realtype), dimension(*), intent(in) :: t  !
        !      Local variables.
        !
        integer(kind=intType) :: nn
 
        ! Loop over the number of blocks.
 
        domains: do nn = 1, ndom
 
            ! Set the pointers for this block.
 
            call setpointers(nn, groundlevel, sps)
            if (.NOT. usetsinterpolatedgridvelocity) then
                call gridvelocitiesfinelevel_block(useoldcoor, t, sps, nn)
            else
                call gridvelocitiesfinelevel_ts_block(nn, sps)
            end if
 
        end do domains
 
    end subroutine gridvelocitiesfinelevel
#endif
#ifndef USE_TAPENADE
    subroutine gridvelocitiesfinelevel_ts_block(nn, sps)
 
        use precision
        use constants
        use blockpointers
        use inputphysics, only: machgrid, veldirfreestream
        use flowvarrefstate, only: gammainf, pinf, rhoinf
        use inputtimespectral, only: dscalar, ntimeintervalsspectral
 
        integer(kind=intType), intent(in) :: nn, sps
        integer :: i, j, k, mm, ii, ie_l, je_l, ke_l
        real(kind=realtype) :: x_vc, y_vc, z_vc
        real(kind=realtype) :: x_fc, y_fc, z_fc
        real(kind=realtype) :: ainf
        real(kind=realtype) :: velxfreestream, velyfreestream, velzfreestream
 
        ! get the grid free stream velocity
        ainf = sqrt(gammainf * pinf / rhoinf)
        velxfreestream = (ainf * machgrid) * (-veldirfreestream(1))
        velyfreestream = (ainf * machgrid) * (-veldirfreestream(2))
        velzfreestream = (ainf * machgrid) * (-veldirfreestream(3))
 
        ! Grid velocities of the cell centers, including the
        ! 1st level halo cells.
 
        ! Initialize with free stream velocity
        ie_l = flowdoms(nn, 1, sps)%ie
        je_l = flowdoms(nn, 1, sps)%je
        ke_l = flowdoms(nn, 1, sps)%ke
 
        do k = 1, ke_l
            do j = 1, je_l
                do i = 1, ie_l
 
                    s(i, j, k, 1) = velxfreestream
                    s(i, j, k, 2) = velyfreestream
                    s(i, j, k, 3) = velzfreestream
 
                end do
            end do
        end do
 
        ! The velocity contributed from mesh deformation
        do mm = 1, ntimeintervalsspectral
 
            ie_l = flowdoms(nn, 1, mm)%ie
            je_l = flowdoms(nn, 1, mm)%je
            ke_l = flowdoms(nn, 1, mm)%ke
 
            do k = 1, ke_l
                do j = 1, je_l
                    do i = 1, ie_l
 
                        x_vc = eighth * (flowdoms(nn, 1, mm)%x(i - 1, j - 1, k - 1, 1) + &
                                         flowdoms(nn, 1, mm)%x(i, j - 1, k - 1, 1) &
                                         + flowdoms(nn, 1, mm)%x(i - 1, j, k - 1, 1) + &
                                         flowdoms(nn, 1, mm)%x(i, j, k - 1, 1) &
                                         + flowdoms(nn, 1, mm)%x(i - 1, j - 1, k, 1) + &
                                         flowdoms(nn, 1, mm)%x(i, j - 1, k, 1) &
                                         + flowdoms(nn, 1, mm)%x(i - 1, j, k, 1) + &
                                         flowdoms(nn, 1, mm)%x(i, j, k, 1))
 
                        y_vc = eighth * (flowdoms(nn, 1, mm)%x(i - 1, j - 1, k - 1, 2) + &
                                         flowdoms(nn, 1, mm)%x(i, j - 1, k - 1, 2) &
                                         + flowdoms(nn, 1, mm)%x(i - 1, j, k - 1, 2) + &
                                         flowdoms(nn, 1, mm)%x(i, j, k - 1, 2) &
                                         + flowdoms(nn, 1, mm)%x(i - 1, j - 1, k, 2) + &
                                         flowdoms(nn, 1, mm)%x(i, j - 1, k, 2) &
                                         + flowdoms(nn, 1, mm)%x(i - 1, j, k, 2) + &
                                         flowdoms(nn, 1, mm)%x(i, j, k, 2))
 
                        z_vc = eighth * (flowdoms(nn, 1, mm)%x(i - 1, j - 1, k - 1, 3) + &
                                         flowdoms(nn, 1, mm)%x(i, j - 1, k - 1, 3) &
                                         + flowdoms(nn, 1, mm)%x(i - 1, j, k - 1, 3) + &
                                         flowdoms(nn, 1, mm)%x(i, j, k - 1, 3) &
                                         + flowdoms(nn, 1, mm)%x(i - 1, j - 1, k, 3) + &
                                         flowdoms(nn, 1, mm)%x(i, j - 1, k, 3) &
                                         + flowdoms(nn, 1, mm)%x(i - 1, j, k, 3) + &
                                         flowdoms(nn, 1, mm)%x(i, j, k, 3))
 
                        s(i, j, k, 1) = s(i, j, k, 1) + dscalar(1, sps, mm) * x_vc
                        s(i, j, k, 2) = s(i, j, k, 2) + dscalar(1, sps, mm) * y_vc
                        s(i, j, k, 3) = s(i, j, k, 3) + dscalar(1, sps, mm) * z_vc
 
                    end do
                end do
            end do
 
        end do
 
        ! Normal grid velocities of the faces.
 
        ! sFaceI=dot(sI, v)
        ! =dot(sI, v_freestream + v_meshmotion)
        ! =dot(sI, v_freestream) + dot(sI, v_meshmotion)
 
        ! sFaceJ, sFaceK follow the same rule.
 
        ! dot(sI, v_freestream)
        ie_l = flowdoms(nn, 1, sps)%ie
        je_l = flowdoms(nn, 1, sps)%je
        ke_l = flowdoms(nn, 1, sps)%ke
 
        ! i
        do k = 1, ke_l
            do j = 1, je_l
                do i = 0, ie_l
 
                    sfacei(i, j, k) = velxfreestream * si(i, j, k, 1) + velyfreestream * si(i, j, k, 2) &
                                      + velzfreestream * si(i, j, k, 3)
 
                end do
            end do
        end do
 
        ! j
        do k = 1, ke_l
            do j = 0, je_l
                do i = 1, ie_l
 
                    sfacej(i, j, k) = velxfreestream * sj(i, j, k, 1) + velyfreestream * sj(i, j, k, 2) &
                                      + velzfreestream * sj(i, j, k, 3)
 
                end do
            end do
        end do
 
        ! k
        do k = 0, ke_l
            do j = 1, je_l
                do i = 1, ie_l
 
                    sfacek(i, j, k) = velxfreestream * sk(i, j, k, 1) + velyfreestream * sk(i, j, k, 2) &
                                      + velzfreestream * sk(i, j, k, 3)
 
                end do
            end do
        end do
 
        !  dot(sI, v_meshmotion)
 
        do mm = 1, ntimeintervalsspectral
 
            ie_l = flowdoms(nn, 1, mm)%ie
            je_l = flowdoms(nn, 1, mm)%je
            ke_l = flowdoms(nn, 1, mm)%ke
 
            ! i
            do k = 1, ke_l
                do j = 1, je_l
                    do i = 0, ie_l
 
                        x_fc = fourth * (flowdoms(nn, 1, mm)%x(i, j - 1, k - 1, 1) + &
                                         flowdoms(nn, 1, mm)%x(i, j, k, 1) &
                                         + flowdoms(nn, 1, mm)%x(i, j - 1, k, 1) + &
                                         flowdoms(nn, 1, mm)%x(i, j, k - 1, 1))
                        y_fc = fourth * (flowdoms(nn, 1, mm)%x(i, j - 1, k - 1, 2) + &
                                         flowdoms(nn, 1, mm)%x(i, j, k, 2) &
                                         + flowdoms(nn, 1, mm)%x(i, j - 1, k, 2) + &
                                         flowdoms(nn, 1, mm)%x(i, j, k - 1, 2))
                        z_fc = fourth * (flowdoms(nn, 1, mm)%x(i, j - 1, k - 1, 3) + &
                                         flowdoms(nn, 1, mm)%x(i, j, k, 3) &
                                         + flowdoms(nn, 1, mm)%x(i, j - 1, k, 3) + &
                                         flowdoms(nn, 1, mm)%x(i, j, k - 1, 3))
 
                        sfacei(i, j, k) = sfacei(i, j, k) &
                                          + dscalar(1, sps, mm) * x_fc * si(i, j, k, 1) &
                                          + dscalar(1, sps, mm) * y_fc * si(i, j, k, 2) &
                                          + dscalar(1, sps, mm) * z_fc * si(i, j, k, 3)
 
                    end do
                end do
            end do
 
            ! j
            do k = 1, ke_l
                do j = 0, je_l
                    do i = 1, ie_l
 
                        x_fc = fourth * (flowdoms(nn, 1, mm)%x(i - 1, j, k - 1, 1) + &
                                         flowdoms(nn, 1, mm)%x(i, j, k, 1) &
                                         + flowdoms(nn, 1, mm)%x(i - 1, j, k, 1) + &
                                         flowdoms(nn, 1, mm)%x(i, j, k - 1, 1))
                        y_fc = fourth * (flowdoms(nn, 1, mm)%x(i - 1, j, k - 1, 2) + &
                                         flowdoms(nn, 1, mm)%x(i, j, k, 2) &
                                         + flowdoms(nn, 1, mm)%x(i - 1, j, k, 2) + &
                                         flowdoms(nn, 1, mm)%x(i, j, k - 1, 2))
                        z_fc = fourth * (flowdoms(nn, 1, mm)%x(i - 1, j, k - 1, 3) + &
                                         flowdoms(nn, 1, mm)%x(i, j, k, 3) &
                                         + flowdoms(nn, 1, mm)%x(i - 1, j, k, 3) + &
                                         flowdoms(nn, 1, mm)%x(i, j, k - 1, 3))
 
                        sfacej(i, j, k) = sfacej(i, j, k) &
                                          + dscalar(1, sps, mm) * x_fc * sj(i, j, k, 1) &
                                          + dscalar(1, sps, mm) * y_fc * sj(i, j, k, 2) &
                                          + dscalar(1, sps, mm) * z_fc * sj(i, j, k, 3)
 
                    end do
                end do
            end do
 
            ! k
            do k = 0, ke_l
                do j = 1, je_l
                    do i = 1, ie_l
 
                        x_fc = fourth * (flowdoms(nn, 1, mm)%x(i - 1, j - 1, k, 1) + &
                                         flowdoms(nn, 1, mm)%x(i, j, k, 1) &
                                         + flowdoms(nn, 1, mm)%x(i, j - 1, k, 1) + &
                                         flowdoms(nn, 1, mm)%x(i - 1, j, k, 1))
                        y_fc = fourth * (flowdoms(nn, 1, mm)%x(i - 1, j - 1, k, 2) + &
                                         flowdoms(nn, 1, mm)%x(i, j, k, 2) &
                                         + flowdoms(nn, 1, mm)%x(i, j - 1, k, 2) + &
                                         flowdoms(nn, 1, mm)%x(i - 1, j, k, 2))
                        z_fc = fourth * (flowdoms(nn, 1, mm)%x(i - 1, j - 1, k, 3) + &
                                         flowdoms(nn, 1, mm)%x(i, j, k, 3) &
                                         + flowdoms(nn, 1, mm)%x(i, j - 1, k, 3) + &
                                         flowdoms(nn, 1, mm)%x(i - 1, j, k, 3))
 
                        sfacek(i, j, k) = sfacek(i, j, k) &
                                          + dscalar(1, sps, mm) * x_fc * sk(i, j, k, 1) &
                                          + dscalar(1, sps, mm) * y_fc * sk(i, j, k, 2) &
                                          + dscalar(1, sps, mm) * z_fc * sk(i, j, k, 3)
 
                    end do
                end do
            end do
 
        end do
 
    end subroutine gridvelocitiesfinelevel_ts_block
#endif
    subroutine gridvelocitiesfinelevel_block(useOldCoor, t, sps, nn)
        !
        !       gridVelocitiesFineLevel computes the grid velocities for
        !       the cell centers and the normal grid velocities for the faces
        !       of moving blocks for the currently finest grid, i.e.
        !       groundLevel. The velocities are computed at time t for
        !       spectral mode sps. If useOldCoor is .true. the velocities
        !       are determined using the unsteady time integrator in
        !       combination with the old coordinates; otherwise the analytic
        !       form is used.
        !
        use blockpointers
        use cgnsgrid
        use flowvarrefstate
        use inputmotion
        use inputunsteady
        use iteration
        use inputphysics
        use inputtsstabderiv
        use monitor
        use communication
        use flowutils, only: derivativerotmatrixrigid, getdirvector
        use utils, only: setcoeftimeintegrator, tsalpha, tsbeta, tsmach, terminate, &
                         rotmatrixrigidbody, getdirangle
 
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: sps, nn
        logical, intent(in) :: useOldCoor
 
        real(kind=realtype), dimension(*), intent(in) :: t
        !
        !      Local variables.
        !
        integer(kind=intType) :: mm
        integer(kind=intType) :: i, j, k, ii, iie, jje, kke
 
        real(kind=realtype) :: oneover4dt, oneover8dt
        real(kind=realtype) :: velxgrid, velygrid, velzgrid, ainf
        real(kind=realtype) :: velxgrid0, velygrid0, velzgrid0
 
        real(kind=realtype), dimension(3) :: sc, xc, xxc
        real(kind=realtype), dimension(3) :: rotcenter, rotrate
 
        real(kind=realtype), dimension(3) :: rotationpoint
        real(kind=realtype), dimension(3, 3) :: rotationmatrix, &
                                                derivrotationmatrix
 
        real(kind=realtype) :: tnew, told
        real(kind=realtype), dimension(:, :), pointer :: sface
 
        real(kind=realtype), dimension(:, :, :), pointer :: xx, ss
        real(kind=realtype), dimension(:, :, :, :), pointer :: xxold
 
        real(kind=realtype) :: intervalmach, alphats, alphaincrement, &
                               betats, betaincrement
        real(kind=realtype), dimension(3) :: veldir
        real(kind=realtype), dimension(3) :: refdirection
 
        ! Compute the mesh velocity from the given mesh Mach number.
 
        ! vel{x,y,z}Grid0 is the ACTUAL velocity you want at the
        ! geometry.
        ainf = sqrt(gammainf * pinf / rhoinf)
        velxgrid0 = (ainf * machgrid) * (-veldirfreestream(1))
        velygrid0 = (ainf * machgrid) * (-veldirfreestream(2))
        velzgrid0 = (ainf * machgrid) * (-veldirfreestream(3))
 
        ! Compute the derivative of the rotation matrix and the rotation
        ! point; needed for velocity due to the rigid body rotation of
        ! the entire grid. It is assumed that the rigid body motion of
        ! the grid is only specified if there is only 1 section present.
 
        call derivativerotmatrixrigid(derivrotationmatrix, rotationpoint, t(1))
 
        !compute the rotation matrix to update the velocities for the time
        !spectral stability derivative case...
 
#ifndef USE_TAPENADE
! We do not differentiate the time spectral stability stuff for now
        if (tsstability) then
            ! Determine the time values of the old and new time level.
            ! It is assumed that the rigid body rotation of the mesh is only
            ! used when only 1 section is present.
 
            tnew = timeunsteady + timeunsteadyrestart
            told = tnew - t(1)
 
            if (tspmode .or. tsqmode .or. tsrmode) then
                ! Compute the rotation matrix of the rigid body rotation as
                ! well as the rotation point; the latter may vary in time due
                ! to rigid body translation.
 
                call rotmatrixrigidbody(tnew, told, rotationmatrix, rotationpoint)
 
                if (tsalphafollowing) then
 
                    velxgrid0 = rotationmatrix(1, 1) * velxgrid0 &
                                + rotationmatrix(1, 2) * velygrid0 &
                                + rotationmatrix(1, 3) * velzgrid0
                    velygrid0 = rotationmatrix(2, 1) * velxgrid0 &
                                + rotationmatrix(2, 2) * velygrid0 &
                                + rotationmatrix(2, 3) * velzgrid0
                    velzgrid0 = rotationmatrix(3, 1) * velxgrid0 &
                                + rotationmatrix(3, 2) * velygrid0 &
                                + rotationmatrix(3, 3) * velzgrid0
 
                end if
 
            elseif (tsalphamode) then
 
                !Determine the alpha for this time instance
                alphaincrement = tsalpha(degreepolalpha, coefpolalpha, &
                                         degreefouralpha, omegafouralpha, &
                                         coscoeffouralpha, sincoeffouralpha, t(1))
 
                alphats = alpha + alphaincrement
                !Determine the grid velocity for this alpha
                refdirection(:) = zero
                refdirection(1) = one
                call getdirvector(refdirection, alphats, beta, veldir, liftindex)
 
                !do I need to update the lift direction and drag direction as well?
                !set the effictive grid velocity for this time interval
                velxgrid0 = (ainf * machgrid) * (-veldir(1))
                velygrid0 = (ainf * machgrid) * (-veldir(2))
                velzgrid0 = (ainf * machgrid) * (-veldir(3))
 
            elseif (tsbetamode) then
 
                !Determine the alpha for this time instance
                betaincrement = tsbeta(degreepolbeta, coefpolbeta, &
                                       degreefourbeta, omegafourbeta, &
                                       coscoeffourbeta, sincoeffourbeta, t(1))
 
                betats = beta + betaincrement
                !Determine the grid velocity for this alpha
                refdirection(:) = zero
                refdirection(1) = one
                call getdirvector(refdirection, alpha, betats, veldir, liftindex)
 
                !do I need to update the lift direction and drag direction as well?
                !set the effictive grid velocity for this time interval
                velxgrid0 = (ainf * machgrid) * (-veldir(1))
                velygrid0 = (ainf * machgrid) * (-veldir(2))
                velzgrid0 = (ainf * machgrid) * (-veldir(3))
            elseif (tsmachmode) then
                !determine the mach number at this time interval
                intervalmach = tsmach(degreepolmach, coefpolmach, &
                                      degreefourmach, omegafourmach, &
                                      coscoeffourmach, sincoeffourmach, t(1))
                !set the effective grid velocity
                velxgrid0 = (ainf * (intervalmach + machgrid)) * (-veldirfreestream(1))
                velygrid0 = (ainf * (intervalmach + machgrid)) * (-veldirfreestream(2))
                velzgrid0 = (ainf * (intervalmach + machgrid)) * (-veldirfreestream(3))
 
            elseif (tsaltitudemode) then
                call terminate('gridVelocityFineLevel', 'altitude motion not yet implemented...')
            else
                call terminate('gridVelocityFineLevel', 'Not a recognized Stability Motion')
            end if
        end if
#endif
 
        testmoving: if (blockismoving) then
            ! Determine the situation we are having here.
 
            testuseoldcoor: if (useoldcoor) then
#ifndef USE_TAPENADE
! We do not consider the finite difference based velocity computation for now
                !
                !             The velocities must be determined via a finite
                !             difference formula using the coordinates of the old
                !             levels.
                !
                ! Set the coefficients for the time integrator and store
                ! the inverse of the physical nonDimensional time step,
                ! divided by 4 and 8, a bit easier.
 
                call setcoeftimeintegrator
                oneover4dt = fourth * timeref / deltat
                oneover8dt = half * oneover4dt
                !
                !             Grid velocities of the cell centers, including the
                !             1st level halo cells.
                !
                ! Loop over the cells, including the 1st level halo's.
 
                do k = 1, ke
                    do j = 1, je
                        do i = 1, ie
 
                            ! The velocity of the cell center is determined
                            ! by a finite difference formula. First store
                            ! the current coordinate, multiplied by 8 and
                            ! coefTime(0) in sc.
 
                            sc(1) = (x(i - 1, j - 1, k - 1, 1) + x(i, j - 1, k - 1, 1) &
                                     + x(i - 1, j, k - 1, 1) + x(i, j, k - 1, 1) &
                                     + x(i - 1, j - 1, k, 1) + x(i, j - 1, k, 1) &
                                     + x(i - 1, j, k, 1) + x(i, j, k, 1)) &
                                    * coeftime(0)
                            sc(2) = (x(i - 1, j - 1, k - 1, 2) + x(i, j - 1, k - 1, 2) &
                                     + x(i - 1, j, k - 1, 2) + x(i, j, k - 1, 2) &
                                     + x(i - 1, j - 1, k, 2) + x(i, j - 1, k, 2) &
                                     + x(i - 1, j, k, 2) + x(i, j, k, 2)) &
                                    * coeftime(0)
                            sc(3) = (x(i - 1, j - 1, k - 1, 3) + x(i, j - 1, k - 1, 3) &
                                     + x(i - 1, j, k - 1, 3) + x(i, j, k - 1, 3) &
                                     + x(i - 1, j - 1, k, 3) + x(i, j - 1, k, 3) &
                                     + x(i - 1, j, k, 3) + x(i, j, k, 3)) &
                                    * coeftime(0)
 
                            ! Loop over the older levels to complete the
                            ! finite difference formula.
 
                            do ii = 1, noldlevels
                                sc(1) = sc(1) + (xold(ii, i - 1, j - 1, k - 1, 1) &
                                                 + xold(ii, i, j - 1, k - 1, 1) &
                                                 + xold(ii, i - 1, j, k - 1, 1) &
                                                 + xold(ii, i, j, k - 1, 1) &
                                                 + xold(ii, i - 1, j - 1, k, 1) &
                                                 + xold(ii, i, j - 1, k, 1) &
                                                 + xold(ii, i - 1, j, k, 1) &
                                                 + xold(ii, i, j, k, 1)) &
                                        * coeftime(ii)
                                sc(2) = sc(2) + (xold(ii, i - 1, j - 1, k - 1, 2) &
                                                 + xold(ii, i, j - 1, k - 1, 2) &
                                                 + xold(ii, i - 1, j, k - 1, 2) &
                                                 + xold(ii, i, j, k - 1, 2) &
                                                 + xold(ii, i - 1, j - 1, k, 2) &
                                                 + xold(ii, i, j - 1, k, 2) &
                                                 + xold(ii, i - 1, j, k, 2) &
                                                 + xold(ii, i, j, k, 2)) &
                                        * coeftime(ii)
                                sc(3) = sc(3) + (xold(ii, i - 1, j - 1, k - 1, 3) &
                                                 + xold(ii, i, j - 1, k - 1, 3) &
                                                 + xold(ii, i - 1, j, k - 1, 3) &
                                                 + xold(ii, i, j, k - 1, 3) &
                                                 + xold(ii, i - 1, j - 1, k, 3) &
                                                 + xold(ii, i, j - 1, k, 3) &
                                                 + xold(ii, i - 1, j, k, 3) &
                                                 + xold(ii, i, j, k, 3)) &
                                        * coeftime(ii)
                            end do
 
                            ! Divide by 8 delta t to obtain the correct
                            ! velocities.
 
                            s(i, j, k, 1) = sc(1) * oneover8dt
                            s(i, j, k, 2) = sc(2) * oneover8dt
                            s(i, j, k, 3) = sc(3) * oneover8dt
                        end do
                    end do
                end do
                !
                !             Normal grid velocities of the faces.
                !
                ! Loop over the three directions.
 
                loopdir: do mm = 1, 3
 
                    ! Set the upper boundaries depending on the direction.
 
                    select case (mm)
                    case (1_inttype)       ! normals in i-direction
                        iie = ie; jje = je; kke = ke
 
                    case (2_inttype)       ! normals in j-direction
                        iie = je; jje = ie; kke = ke
 
                    case (3_inttype)       ! normals in k-direction
                        iie = ke; jje = ie; kke = je
                    end select
                    !
                    !               Normal grid velocities in generalized i-direction.
                    !               Mm == 1: i-direction
                    !               mm == 2: j-direction
                    !               mm == 3: k-direction
                    !
                    do i = 0, iie
 
                        ! Set the pointers for the coordinates, normals and
                        ! normal velocities for this generalized i-plane.
                        ! This depends on the value of mm.
 
                        select case (mm)
                        case (1_inttype)       ! normals in i-direction
                            xx => x(i, :, :, :); xxold => xold(:, i, :, :, :)
                            ss => si(i, :, :, :); sface => sfacei(i, :, :)
 
                        case (2_inttype)       ! normals in j-direction
                            xx => x(:, i, :, :); xxold => xold(:, :, i, :, :)
                            ss => sj(:, i, :, :); sface => sfacej(:, i, :)
 
                        case (3_inttype)       ! normals in k-direction
                            xx => x(:, :, i, :); xxold => xold(:, :, :, i, :)
                            ss => sk(:, :, i, :); sface => sfacek(:, :, i)
                        end select
 
                        ! Loop over the k and j-direction of this
                        ! generalized i-face. Note that due to the usage of
                        ! the pointers xx and xxOld an offset of +1 must be
                        ! used in the coordinate arrays, because x and xOld
                        ! originally start at 0 for the i, j and k indices.
 
                        do k = 1, kke
                            do j = 1, jje
 
                                ! The velocity of the face center is determined
                                ! by a finite difference formula. First store
                                ! the current coordinate, multiplied by 4 and
                                ! coefTime(0) in sc.
 
                                sc(1) = coeftime(0) * (xx(j + 1, k + 1, 1) + xx(j, k + 1, 1) &
                                                       + xx(j + 1, k, 1) + xx(j, k, 1))
                                sc(2) = coeftime(0) * (xx(j + 1, k + 1, 2) + xx(j, k + 1, 2) &
                                                       + xx(j + 1, k, 2) + xx(j, k, 2))
                                sc(3) = coeftime(0) * (xx(j + 1, k + 1, 3) + xx(j, k + 1, 3) &
                                                       + xx(j + 1, k, 3) + xx(j, k, 3))
 
                                ! Loop over the older levels to complete the
                                ! finite difference.
 
                                do ii = 1, noldlevels
 
                                    sc(1) = sc(1) + coeftime(ii) &
                                            * (xxold(ii, j + 1, k + 1, 1) &
                                               + xxold(ii, j, k + 1, 1) &
                                               + xxold(ii, j + 1, k, 1) &
                                               + xxold(ii, j, k, 1))
                                    sc(2) = sc(2) + coeftime(ii) &
                                            * (xxold(ii, j + 1, k + 1, 2) &
                                               + xxold(ii, j, k + 1, 2) &
                                               + xxold(ii, j + 1, k, 2) &
                                               + xxold(ii, j, k, 2))
                                    sc(3) = sc(3) + coeftime(ii) &
                                            * (xxold(ii, j + 1, k + 1, 3) &
                                               + xxold(ii, j, k + 1, 3) &
                                               + xxold(ii, j + 1, k, 3) &
                                               + xxold(ii, j, k, 3))
                                end do
 
                                ! Determine the dot product of sc and the normal
                                ! and divide by 4 deltaT to obtain the correct
                                ! value of the normal velocity.
 
                                sface(j, k) = sc(1) * ss(j, k, 1) + sc(2) * ss(j, k, 2) &
                                              + sc(3) * ss(j, k, 3)
                                sface(j, k) = sface(j, k) * oneover4dt
 
                            end do
                        end do
                    end do
 
                end do loopdir
#endif
                else testuseoldcoor
                !
                !             The velocities must be determined analytically.
                !
                ! Store the rotation center and determine the
                ! nonDimensional rotation rate of this block. As the
                ! reference length is 1 timeRef == 1/uRef and at the end
                ! the nonDimensional velocity is computed.
 
                j = nbkglobal
 
                rotcenter = cgnsdoms(j)%rotCenter
                rotrate = timeref * cgnsdoms(j)%rotRate
 
                velxgrid = velxgrid0
                velygrid = velygrid0
                velzgrid = velzgrid0
                !
                !             Grid velocities of the cell centers, including the
                !             1st level halo cells.
                !
                ! Loop over the cells, including the 1st level halo's.
 
                do k = 1, ke
                    do j = 1, je
                        do i = 1, ie
 
                            ! Determine the coordinates of the cell center,
                            ! which are stored in xc.
 
                            xc(1) = eighth * (flowdoms(nn, groundlevel, sps)%x(i - 1, j - 1, k - 1, 1) &
                                              + flowdoms(nn, groundlevel, sps)%x(i, j - 1, k - 1, 1) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, j, k - 1, 1) &
                                              + flowdoms(nn, groundlevel, sps)%x(i, j, k - 1, 1) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, j - 1, k, 1) &
                                              + flowdoms(nn, groundlevel, sps)%x(i, j - 1, k, 1) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, j, k, 1) &
                                              + flowdoms(nn, groundlevel, sps)%x(i, j, k, 1))
                            xc(2) = eighth * (flowdoms(nn, groundlevel, sps)%x(i - 1, j - 1, k - 1, 2) &
                                              + flowdoms(nn, groundlevel, sps)%x(i, j - 1, k - 1, 2) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, j, k - 1, 2) &
                                              + flowdoms(nn, groundlevel, sps)%x(i, j, k - 1, 2) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, j - 1, k, 2) &
                                              + flowdoms(nn, groundlevel, sps)%x(i, j - 1, k, 2) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, j, k, 2) &
                                              + flowdoms(nn, groundlevel, sps)%x(i, j, k, 2))
                            xc(3) = eighth * (flowdoms(nn, groundlevel, sps)%x(i - 1, j - 1, k - 1, 3) &
                                              + flowdoms(nn, groundlevel, sps)%x(i, j - 1, k - 1, 3) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, j, k - 1, 3) &
                                              + flowdoms(nn, groundlevel, sps)%x(i, j, k - 1, 3) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, j - 1, k, 3) &
                                              + flowdoms(nn, groundlevel, sps)%x(i, j - 1, k, 3) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, j, k, 3) &
                                              + flowdoms(nn, groundlevel, sps)%x(i, j, k, 3))
 
                            ! Determine the coordinates relative to the
                            ! center of rotation.
 
                            xxc(1) = xc(1) - rotcenter(1)
                            xxc(2) = xc(2) - rotcenter(2)
                            xxc(3) = xc(3) - rotcenter(3)
 
                            ! Determine the rotation speed of the cell center,
                            ! which is omega*r.
 
                            sc(1) = rotrate(2) * xxc(3) - rotrate(3) * xxc(2)
                            sc(2) = rotrate(3) * xxc(1) - rotrate(1) * xxc(3)
                            sc(3) = rotrate(1) * xxc(2) - rotrate(2) * xxc(1)
 
                            ! Determine the coordinates relative to the
                            ! rigid body rotation point.
 
                            xxc(1) = xc(1) - rotationpoint(1)
                            xxc(2) = xc(2) - rotationpoint(2)
                            xxc(3) = xc(3) - rotationpoint(3)
 
                            ! Determine the total velocity of the cell center.
                            ! This is a combination of rotation speed of this
                            ! block and the entire rigid body rotation.
 
                            s(i, j, k, 1) = sc(1) + velxgrid &
                                            + derivrotationmatrix(1, 1) * xxc(1) &
                                            + derivrotationmatrix(1, 2) * xxc(2) &
                                            + derivrotationmatrix(1, 3) * xxc(3)
                            s(i, j, k, 2) = sc(2) + velygrid &
                                            + derivrotationmatrix(2, 1) * xxc(1) &
                                            + derivrotationmatrix(2, 2) * xxc(2) &
                                            + derivrotationmatrix(2, 3) * xxc(3)
                            s(i, j, k, 3) = sc(3) + velzgrid &
                                            + derivrotationmatrix(3, 1) * xxc(1) &
                                            + derivrotationmatrix(3, 2) * xxc(2) &
                                            + derivrotationmatrix(3, 3) * xxc(3)
                        end do
                    end do
                end do
                !
                !             Normal grid velocities of the faces.
                !
                ! Loop over the three directions.
 
                ! The original code is elegant but the Tapenade has a difficult time
                ! to understand it. Thus, we unfold it and make it easier for the
                ! Tapenade.
 
                ! i-direction
                do k = 1, ke
                    do j = 1, je
                        do i = 0, ie
 
                            ! Determine the coordinates of the face center,
                            ! which are stored in xc.
 
                            xc(1) = fourth * (flowdoms(nn, groundlevel, sps)%x(i, j - 1, k - 1, 1) + &
                                              flowdoms(nn, groundlevel, sps)%x(i, j, k - 1, 1) &
                                              + flowdoms(nn, groundlevel, sps)%x(i, j - 1, k, 1) + &
                                              flowdoms(nn, groundlevel, sps)%x(i, j, k, 1))
                            xc(2) = fourth * (flowdoms(nn, groundlevel, sps)%x(i, j - 1, k - 1, 2) + &
                                              flowdoms(nn, groundlevel, sps)%x(i, j, k - 1, 2) &
                                              + flowdoms(nn, groundlevel, sps)%x(i, j - 1, k, 2) + &
                                              flowdoms(nn, groundlevel, sps)%x(i, j, k, 2))
                            xc(3) = fourth * (flowdoms(nn, groundlevel, sps)%x(i, j - 1, k - 1, 3) + &
                                              flowdoms(nn, groundlevel, sps)%x(i, j, k - 1, 3) &
                                              + flowdoms(nn, groundlevel, sps)%x(i, j - 1, k, 3) + &
                                              flowdoms(nn, groundlevel, sps)%x(i, j, k, 3))
 
                            call cellfacevelocities(xc, rotcenter, rotrate, &
                                                    velxgrid, velygrid, velzgrid, derivrotationmatrix, sc)
 
                            ! Store the dot product of grid velocity sc and
                            ! the normal ss in sFace.
 
                            sfacei(i, j, k) = sc(1) * si(i, j, k, 1) + sc(2) * si(i, j, k, 2) &
                                              + sc(3) * si(i, j, k, 3)
                        end do
                    end do
                end do
 
                ! j-direction
                do k = 1, ke
                    do j = 0, je
                        do i = 1, ie
 
                            ! Determine the coordinates of the face center,
                            ! which are stored in xc.
 
                            xc(1) = fourth * (flowdoms(nn, groundlevel, sps)%x(i - 1, j, k, 1) + &
                                              flowdoms(nn, groundlevel, sps)%x(i, j, k - 1, 1) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, j, k - 1, 1) + &
                                              flowdoms(nn, groundlevel, sps)%x(i, j, k, 1))
                            xc(2) = fourth * (flowdoms(nn, groundlevel, sps)%x(i - 1, j, k, 2) + &
                                              flowdoms(nn, groundlevel, sps)%x(i, j, k - 1, 2) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, j, k - 1, 2) + &
                                              flowdoms(nn, groundlevel, sps)%x(i, j, k, 2))
                            xc(3) = fourth * (flowdoms(nn, groundlevel, sps)%x(i - 1, j, k, 3) + &
                                              flowdoms(nn, groundlevel, sps)%x(i, j, k - 1, 3) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, j, k - 1, 3) + &
                                              flowdoms(nn, groundlevel, sps)%x(i, j, k, 3))
 
                            call cellfacevelocities(xc, rotcenter, rotrate, &
                                                    velxgrid, velygrid, velzgrid, derivrotationmatrix, sc)
 
                            ! Store the dot product of grid velocity sc and
                            ! the normal ss in sFace.
 
                            sfacej(i, j, k) = sc(1) * sj(i, j, k, 1) + sc(2) * sj(i, j, k, 2) &
                                              + sc(3) * sj(i, j, k, 3)
                        end do
                    end do
                end do
 
                ! k-direction
                do k = 0, ke
                    do j = 1, je
                        do i = 1, ie
 
                            ! Determine the coordinates of the face center,
                            ! which are stored in xc.
 
                            xc(1) = fourth * (flowdoms(nn, groundlevel, sps)%x(i, j - 1, k, 1) + &
                                              flowdoms(nn, groundlevel, sps)%x(i - 1, j, k, 1) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, j - 1, k, 1) + &
                                              flowdoms(nn, groundlevel, sps)%x(i, j, k, 1))
                            xc(2) = fourth * (flowdoms(nn, groundlevel, sps)%x(i, j - 1, k, 2) + &
                                              flowdoms(nn, groundlevel, sps)%x(i - 1, j, k, 2) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, j - 1, k, 2) + &
                                              flowdoms(nn, groundlevel, sps)%x(i, j, k, 2))
                            xc(3) = fourth * (flowdoms(nn, groundlevel, sps)%x(i, j - 1, k, 3) + &
                                              flowdoms(nn, groundlevel, sps)%x(i - 1, j, k, 3) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, j - 1, k, 3) + &
                                              flowdoms(nn, groundlevel, sps)%x(i, j, k, 3))
 
                            call cellfacevelocities(xc, rotcenter, rotrate, &
                                                    velxgrid, velygrid, velzgrid, derivrotationmatrix, sc)
 
                            ! Store the dot product of grid velocity sc and
                            ! the normal ss in sFace.
 
                            sfacek(i, j, k) = sc(1) * sk(i, j, k, 1) + sc(2) * sk(i, j, k, 2) &
                                              + sc(3) * sk(i, j, k, 3)
                        end do
                    end do
                end do
 
            end if testuseoldcoor
        end if testmoving
 
    end subroutine gridvelocitiesfinelevel_block
 
    subroutine cellfacevelocities(xc, rotCenter, rotRate, velxGrid, velyGrid, velzGrid, derivRotationMatrix, sc)
        !
        !  Returns the cell face velocities for a given face center
        !
        use constants
 
        implicit none
        !
        !      Subroutine arguments.
        !
        real(kind=realtype), dimension(3), intent(in) :: xc, rotcenter, rotrate
        real(kind=realtype), intent(in) :: velxgrid, velygrid, velzgrid
        real(kind=realtype), dimension(3, 3), intent(in) :: derivrotationmatrix
        real(kind=realtype), dimension(3), intent(out) :: sc
        !
        !      Local variables.
        !
        real(kind=realtype), dimension(3) :: rotationpoint, xxc
 
        ! Determine the coordinates relative to the
        ! center of rotation.
 
        xxc(1) = xc(1) - rotcenter(1)
        xxc(2) = xc(2) - rotcenter(2)
        xxc(3) = xc(3) - rotcenter(3)
 
        ! Determine the rotation speed of the face center,
        ! which is omega*r.
 
        sc(1) = rotrate(2) * xxc(3) - rotrate(3) * xxc(2)
        sc(2) = rotrate(3) * xxc(1) - rotrate(1) * xxc(3)
        sc(3) = rotrate(1) * xxc(2) - rotrate(2) * xxc(1)
 
        ! Determine the coordinates relative to the
        ! rigid body rotation point.
 
        xxc(1) = xc(1) - rotationpoint(1)
        xxc(2) = xc(2) - rotationpoint(2)
        xxc(3) = xc(3) - rotationpoint(3)
 
        ! Determine the total velocity of the cell face.
        ! This is a combination of rotation speed of this
        ! block and the entire rigid body rotation.
 
        sc(1) = sc(1) + velxgrid &
                + derivrotationmatrix(1, 1) * xxc(1) &
                + derivrotationmatrix(1, 2) * xxc(2) &
                + derivrotationmatrix(1, 3) * xxc(3)
        sc(2) = sc(2) + velygrid &
                + derivrotationmatrix(2, 1) * xxc(1) &
                + derivrotationmatrix(2, 2) * xxc(2) &
                + derivrotationmatrix(2, 3) * xxc(3)
        sc(3) = sc(3) + velzgrid &
                + derivrotationmatrix(3, 1) * xxc(1) &
                + derivrotationmatrix(3, 2) * xxc(2) &
                + derivrotationmatrix(3, 3) * xxc(3)
 
    end subroutine cellfacevelocities
 
#ifndef USE_TAPENADE
    subroutine slipvelocitiesfinelevel(useOldCoor, t, sps)
        !
        ! Shell function to call slipVelocitiesFineLevel on all blocks
        !
        use constants
        use blockpointers, only: ndom
        use inputtimespectral, only: ntimeintervalsspectral
        use iteration, only: groundlevel
        use utils, only: setpointers
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: sps
        logical, intent(in) :: useOldCoor
        real(kind=realtype), dimension(*), intent(in) :: t  !
        !      Local variables.
        !
        integer(kind=intType) :: nn
 
        ! Loop over the number of blocks.
 
        domains: do nn = 1, ndom
 
            ! Set the pointers for this block.
 
            call setpointers(nn, groundlevel, sps)
 
            call slipvelocitiesfinelevel_block(useoldcoor, t, sps, nn)
 
        end do domains
 
    end subroutine slipvelocitiesfinelevel
#endif
 
    subroutine slipvelocitiesfinelevel_block(useOldCoor, t, sps, nn)
        !
        !       slipVelocitiesFineLevel computes the slip velocities for
        !       viscous subfaces on all viscous boundaries on groundLevel for
        !       the given spectral solution. If useOldCoor is .true. the
        !       velocities are determined using the unsteady time integrator;
        !       otherwise the analytic form is used.
        !
        use constants
        use inputtimespectral
        use blockpointers
        use cgnsgrid
        use flowvarrefstate
        use inputmotion
        use inputunsteady
        use iteration
        use inputphysics
        use inputtsstabderiv
        use monitor
        use communication
        use flowutils, only: derivativerotmatrixrigid, getdirvector
        use utils, only: tsalpha, tsbeta, tsmach, terminate, rotmatrixrigidbody, &
                         setcoeftimeintegrator, getdirangle
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: sps, nn
        logical, intent(in) :: useOldCoor
 
        real(kind=realtype), dimension(*), intent(in) :: t
        !
        !      Local variables.
        !
        integer(kind=intType) :: mm, i, j, level, ii
 
        real(kind=realtype) :: oneover4dt
        real(kind=realtype) :: velxgrid, velygrid, velzgrid, ainf
        real(kind=realtype) :: velxgrid0, velygrid0, velzgrid0
 
        real(kind=realtype), dimension(3) :: xc, xxc
        real(kind=realtype), dimension(3) :: rotcenter, rotrate
 
        real(kind=realtype), dimension(3) :: rotationpoint
        real(kind=realtype), dimension(3, 3) :: rotationmatrix, &
                                                derivrotationmatrix
 
        real(kind=realtype) :: tnew, told
 
        real(kind=realtype), dimension(:, :, :), pointer :: uslip
        real(kind=realtype), dimension(:, :, :), pointer :: xface
        real(kind=realtype), dimension(:, :, :, :), pointer :: xfaceold
 
        real(kind=realtype) :: intervalmach, alphats, alphaincrement, &
                               betats, betaincrement
        real(kind=realtype), dimension(3) :: veldir
        real(kind=realtype), dimension(3) :: refdirection
 
        ! Determine the situation we are having here.
 
        testuseoldcoor: if (useoldcoor) then
 
#ifndef USE_TAPENADE
! the time-stepping portion is not ADed.
 
            ! The velocities must be determined via a finite difference
            ! formula using the coordinates of the old levels.
 
            ! Set the coefficients for the time integrator and store the
            ! inverse of the physical nonDimensional time step, divided
            ! by 4, a bit easier.
 
            call setcoeftimeintegrator
            oneover4dt = fourth * timeref / deltat
 
            ! Loop over the number of viscous subfaces.
 
            bocoloop1: do mm = 1, nviscbocos
 
                ! Set the pointer for uSlip to make the code more
                ! readable.
 
                uslip => bcdata(mm)%uSlip
 
                ! Determine the grid face on which the subface is located
                ! and set some variables accordingly.
 
                select case (bcfaceid(mm))
 
                case (imin)
                    xface => x(1, :, :, :); xfaceold => xold(:, 1, :, :, :)
 
                case (imax)
                    xface => x(il, :, :, :); xfaceold => xold(:, il, :, :, :)
 
                case (jmin)
                    xface => x(:, 1, :, :); xfaceold => xold(:, :, 1, :, :)
 
                case (jmax)
                    xface => x(:, jl, :, :); xfaceold => xold(:, :, jl, :, :)
 
                case (kmin)
                    xface => x(:, :, 1, :); xfaceold => xold(:, :, :, 1, :)
 
                case (kmax)
                    xface => x(:, :, kl, :); xfaceold => xold(:, :, :, kl, :)
 
                end select
 
                ! Some boundary faces have a different rotation speed than
                ! the corresponding block. This happens e.g. in the tip gap
                ! region of turboMachinary problems where the casing does
                ! not rotate. As the coordinate difference corresponds to
                ! the rotation rate of the block, a correction must be
                ! computed. Therefore compute the difference in rotation
                ! rate and store the rotation center a bit easier. Note that
                ! the rotation center of subface is taken, because if there
                ! is a difference in rotation rate this info for the subface
                ! must always be specified.
 
                j = nbkglobal
                i = cgnssubface(mm)
 
                rotcenter = cgnsdoms(j)%bocoInfo(i)%rotCenter
                rotrate = timeref * (cgnsdoms(j)%bocoInfo(i)%rotRate &
                                     - cgnsdoms(j)%rotRate)
 
                ! Loop over the quadrilateral faces of the viscous subface.
                ! Note that due to the usage of the pointers xFace and
                ! xFaceOld an offset of +1 must be used in the coordinate
                ! arrays, because x and xOld originally start at 0 for the
                ! i, j and k indices.
 
                do j = bcdata(mm)%jcBeg, bcdata(mm)%jcEnd
                    do i = bcdata(mm)%icBeg, bcdata(mm)%icEnd
 
                        ! Determine the coordinates of the centroid of the
                        ! face, multiplied by 4.
 
                        xc(1) = xface(i + 1, j + 1, 1) + xface(i + 1, j, 1) &
                                + xface(i, j + 1, 1) + xface(i, j, 1)
                        xc(2) = xface(i + 1, j + 1, 2) + xface(i + 1, j, 2) &
                                + xface(i, j + 1, 2) + xface(i, j, 2)
                        xc(3) = xface(i + 1, j + 1, 3) + xface(i + 1, j, 3) &
                                + xface(i, j + 1, 3) + xface(i, j, 3)
 
                        ! Multiply the sum of the 4 vertex coordinates with
                        ! coefTime(0) to obtain the contribution for the
                        ! current time level. The division by 4*deltaT will
                        ! take place later. This is both more efficient and
                        ! more accurate for extremely small time steps.
 
                        uslip(i, j, 1) = coeftime(0) * xc(1)
                        uslip(i, j, 2) = coeftime(0) * xc(2)
                        uslip(i, j, 3) = coeftime(0) * xc(3)
 
                        ! Loop over the older time levels and take their
                        ! contribution into account.
 
                        do level = 1, noldlevels
 
                            uslip(i, j, 1) = uslip(i, j, 1) + coeftime(level) &
                                             * (xfaceold(level, i + 1, j + 1, 1) &
                                                + xfaceold(level, i + 1, j, 1) &
                                                + xfaceold(level, i, j + 1, 1) &
                                                + xfaceold(level, i, j, 1))
 
                            uslip(i, j, 2) = uslip(i, j, 2) + coeftime(level) &
                                             * (xfaceold(level, i + 1, j + 1, 2) &
                                                + xfaceold(level, i + 1, j, 2) &
                                                + xfaceold(level, i, j + 1, 2) &
                                                + xfaceold(level, i, j, 2))
 
                            uslip(i, j, 3) = uslip(i, j, 3) + coeftime(level) &
                                             * (xfaceold(level, i + 1, j + 1, 3) &
                                                + xfaceold(level, i + 1, j, 3) &
                                                + xfaceold(level, i, j + 1, 3) &
                                                + xfaceold(level, i, j, 3))
                        end do
 
                        ! Divide by 4 times the time step to obtain the
                        ! correct velocity.
 
                        uslip(i, j, 1) = uslip(i, j, 1) * oneover4dt
                        uslip(i, j, 2) = uslip(i, j, 2) * oneover4dt
                        uslip(i, j, 3) = uslip(i, j, 3) * oneover4dt
 
                        ! Determine the correction due to the difference
                        ! in rotation rate between the block and subface.
 
                        ! First determine the coordinates relative to the
                        ! rotation center. Remember that 4 times this value
                        ! is currently stored in xc.
 
                        xc(1) = fourth * xc(1) - rotcenter(1)
                        xc(2) = fourth * xc(2) - rotcenter(2)
                        xc(3) = fourth * xc(3) - rotcenter(3)
 
                        ! Compute the velocity, which is the cross product
                        ! of rotRate and xc and add it to uSlip.
 
                        uslip(i, j, 1) = uslip(i, j, 1) &
                                         + rotrate(2) * xc(3) - rotrate(3) * xc(2)
                        uslip(i, j, 2) = uslip(i, j, 2) &
                                         + rotrate(3) * xc(1) - rotrate(1) * xc(3)
                        uslip(i, j, 3) = uslip(i, j, 3) &
                                         + rotrate(1) * xc(2) - rotrate(2) * xc(1)
 
                    end do
                end do
 
            end do bocoloop1
#else
            continue
#endif
        else
 
            ! The velocities must be determined analytically.
 
            ! Compute the mesh velocity from the given mesh Mach number.
 
            !  aInf = sqrt(gammaInf*pInf/rhoInf)
            !  velxGrid = aInf*MachGrid(1)
            !  velyGrid = aInf*MachGrid(2)
            !  velzGrid = aInf*MachGrid(3)
 
            ainf = sqrt(gammainf * pinf / rhoinf)
            velxgrid0 = (ainf * machgrid) * (-veldirfreestream(1))
            velygrid0 = (ainf * machgrid) * (-veldirfreestream(2))
            velzgrid0 = (ainf * machgrid) * (-veldirfreestream(3))
 
            ! Compute the derivative of the rotation matrix and the rotation
            ! point; needed for velocity due to the rigid body rotation of
            ! the entire grid. It is assumed that the rigid body motion of
            ! the grid is only specified if there is only 1 section present.
 
            call derivativerotmatrixrigid(derivrotationmatrix, rotationpoint, &
                                          t(1))
 
            !compute the rotation matrix to update the velocities for the time
            !spectral stability derivative case...
 
#ifndef USE_TAPENADE
! the stability portion is not ADed.
 
            if (tsstability) then
                ! Determine the time values of the old and new time level.
                ! It is assumed that the rigid body rotation of the mesh is only
                ! used when only 1 section is present.
 
                tnew = timeunsteady + timeunsteadyrestart
                told = tnew - t(1)
 
                if (tspmode .or. tsqmode .or. tsrmode) then
                    ! Compute the rotation matrix of the rigid body rotation as
                    ! well as the rotation point; the latter may vary in time due
                    ! to rigid body translation.
 
                    call rotmatrixrigidbody(tnew, told, rotationmatrix, rotationpoint)
 
                    if (tsalphafollowing) then
 
                        velxgrid0 = rotationmatrix(1, 1) * velxgrid0 &
                                    + rotationmatrix(1, 2) * velygrid0 &
                                    + rotationmatrix(1, 3) * velzgrid0
                        velygrid0 = rotationmatrix(2, 1) * velxgrid0 &
                                    + rotationmatrix(2, 2) * velygrid0 &
                                    + rotationmatrix(2, 3) * velzgrid0
                        velzgrid0 = rotationmatrix(3, 1) * velxgrid0 &
                                    + rotationmatrix(3, 2) * velygrid0 &
                                    + rotationmatrix(3, 3) * velzgrid0
 
                    end if
                elseif (tsalphamode) then
                    !Determine the alpha for this time instance
                    alphaincrement = tsalpha(degreepolalpha, coefpolalpha, &
                                             degreefouralpha, omegafouralpha, &
                                             coscoeffouralpha, sincoeffouralpha, t(1))
 
                    alphats = alpha + alphaincrement
                    !Determine the grid velocity for this alpha
                    refdirection(:) = zero
                    refdirection(1) = one
                    call getdirvector(refdirection, alphats, beta, veldir, liftindex)
 
                    !do I need to update the lift direction and drag direction as well?
                    !set the effictive grid velocity for this time interval
                    velxgrid0 = (ainf * machgrid) * (-veldir(1))
                    velygrid0 = (ainf * machgrid) * (-veldir(2))
                    velzgrid0 = (ainf * machgrid) * (-veldir(3))
 
                elseif (tsbetamode) then
 
                    !Determine the alpha for this time instance
                    betaincrement = tsbeta(degreepolbeta, coefpolbeta, &
                                           degreefourbeta, omegafourbeta, &
                                           coscoeffourbeta, sincoeffourbeta, t(1))
 
                    betats = beta + betaincrement
                    !Determine the grid velocity for this alpha
                    refdirection(:) = zero
                    refdirection(1) = one
                    call getdirvector(refdirection, alpha, betats, veldir, liftindex)
 
                    !do I need to update the lift direction and drag direction as well?
                    !set the effictive grid velocity for this time interval
                    velxgrid0 = (ainf * machgrid) * (-veldir(1))
                    velygrid0 = (ainf * machgrid) * (-veldir(2))
                    velzgrid0 = (ainf * machgrid) * (-veldir(3))
                elseif (tsmachmode) then
                    !determine the mach number at this time interval
                    intervalmach = tsmach(degreepolmach, coefpolmach, &
                                          degreefourmach, omegafourmach, &
                                          coscoeffourmach, sincoeffourmach, t(1))
                    !set the effective grid velocity
                    velxgrid0 = (ainf * (intervalmach + machgrid)) * (-veldirfreestream(1))
                    velygrid0 = (ainf * (intervalmach + machgrid)) * (-veldirfreestream(2))
                    velzgrid0 = (ainf * (intervalmach + machgrid)) * (-veldirfreestream(3))
 
                elseif (tsaltitudemode) then
                    call terminate('gridVelocityFineLevel', 'altitude motion not yet implemented...')
                else
                    call terminate('gridVelocityFineLevel', 'Not a recognized Stability Motion')
                end if
            end if
#endif
 
            ! Loop over the number of viscous subfaces.
 
            bocoloop2: do mm = 1, nviscbocos
 
                ! Store the rotation center and the rotation rate
                ! for this subface.
 
                ii = cgnssubface(mm)
 
                rotcenter = cgnsdoms(nbkglobal)%bocoInfo(ii)%rotCenter
                rotrate = timeref * cgnsdoms(nbkglobal)%bocoInfo(ii)%rotRate
 
                ! useWindAxis should go back here!
                velxgrid = velxgrid0
                velygrid = velygrid0
                velzgrid = velzgrid0
 
                ! Loop over the quadrilateral faces of the viscous
                ! subface.
 
                ! The new procedure is less elegant as the previous one.
                ! But the new stands up to Tapenade.
                if (bcfaceid(mm) == imin) then
 
                    do j = bcdata(mm)%jcBeg, bcdata(mm)%jcEnd
                        do i = bcdata(mm)%icBeg, bcdata(mm)%icEnd
 
                            ! Compute the coordinates of the centroid of the face.
                            ! Normally this would be an average of i-1 and i, but
                            ! due to the usage of the pointer xFace and the fact
                            ! that x starts at index 0 this is shifted 1 index.
 
                            xc(1) = fourth * (flowdoms(nn, groundlevel, sps)%x(1, i, j, 1) &
                                              + flowdoms(nn, groundlevel, sps)%x(1, i, j - 1, 1) &
                                              + flowdoms(nn, groundlevel, sps)%x(1, i - 1, j, 1) &
                                              + flowdoms(nn, groundlevel, sps)%x(1, i - 1, j - 1, 1))
                            xc(2) = fourth * (flowdoms(nn, groundlevel, sps)%x(1, i, j, 2) &
                                              + flowdoms(nn, groundlevel, sps)%x(1, i, j - 1, 2) &
                                              + flowdoms(nn, groundlevel, sps)%x(1, i - 1, j, 2) &
                                              + flowdoms(nn, groundlevel, sps)%x(1, i - 1, j - 1, 2))
                            xc(3) = fourth * (flowdoms(nn, groundlevel, sps)%x(1, i, j, 3) &
                                              + flowdoms(nn, groundlevel, sps)%x(1, i, j - 1, 3) &
                                              + flowdoms(nn, groundlevel, sps)%x(1, i - 1, j, 3) &
                                              + flowdoms(nn, groundlevel, sps)%x(1, i - 1, j - 1, 3))
 
                            ! Determine the coordinates relative to the center
                            ! of rotation.
 
                            xxc(1) = xc(1) - rotcenter(1)
                            xxc(2) = xc(2) - rotcenter(2)
                            xxc(3) = xc(3) - rotcenter(3)
 
                            ! Compute the velocity, which is the cross product
                            ! of rotRate and xc.
 
                            bcdata(mm)%uSlip(i, j, 1) = rotrate(2) * xxc(3) - rotrate(3) * xxc(2)
                            bcdata(mm)%uSlip(i, j, 2) = rotrate(3) * xxc(1) - rotrate(1) * xxc(3)
                            bcdata(mm)%uSlip(i, j, 3) = rotrate(1) * xxc(2) - rotrate(2) * xxc(1)
 
                            ! Determine the coordinates relative to the
                            ! rigid body rotation point.
 
                            xxc(1) = xc(1) - rotationpoint(1)
                            xxc(2) = xc(2) - rotationpoint(2)
                            xxc(3) = xc(3) - rotationpoint(3)
 
                            ! Determine the total velocity of the cell center.
                            ! This is a combination of rotation speed of this
                            ! block and the entire rigid body rotation.
 
                            bcdata(mm)%uSlip(i, j, 1) = bcdata(mm)%uSlip(i, j, 1) + velxgrid &
                                                        + derivrotationmatrix(1, 1) * xxc(1) &
                                                        + derivrotationmatrix(1, 2) * xxc(2) &
                                                        + derivrotationmatrix(1, 3) * xxc(3)
                            bcdata(mm)%uSlip(i, j, 2) = bcdata(mm)%uSlip(i, j, 2) + velygrid &
                                                        + derivrotationmatrix(2, 1) * xxc(1) &
                                                        + derivrotationmatrix(2, 2) * xxc(2) &
                                                        + derivrotationmatrix(2, 3) * xxc(3)
                            bcdata(mm)%uSlip(i, j, 3) = bcdata(mm)%uSlip(i, j, 3) + velzgrid &
                                                        + derivrotationmatrix(3, 1) * xxc(1) &
                                                        + derivrotationmatrix(3, 2) * xxc(2) &
                                                        + derivrotationmatrix(3, 3) * xxc(3)
                        end do
                    end do
 
                else if (bcfaceid(mm) == imax) then
 
                    do j = bcdata(mm)%jcBeg, bcdata(mm)%jcEnd
                        do i = bcdata(mm)%icBeg, bcdata(mm)%icEnd
 
                            ! Compute the coordinates of the centroid of the face.
                            ! Normally this would be an average of i-1 and i, but
                            ! due to the usage of the pointer xFace and the fact
                            ! that x starts at index 0 this is shifted 1 index.
 
                            xc(1) = fourth * (flowdoms(nn, groundlevel, sps)%x(il, i, j, 1) &
                                              + flowdoms(nn, groundlevel, sps)%x(il, i, j - 1, 1) &
                                              + flowdoms(nn, groundlevel, sps)%x(il, i - 1, j, 1) &
                                              + flowdoms(nn, groundlevel, sps)%x(il, i - 1, j - 1, 1))
                            xc(2) = fourth * (flowdoms(nn, groundlevel, sps)%x(il, i, j, 2) &
                                              + flowdoms(nn, groundlevel, sps)%x(il, i, j - 1, 2) &
                                              + flowdoms(nn, groundlevel, sps)%x(il, i - 1, j, 2) &
                                              + flowdoms(nn, groundlevel, sps)%x(il, i - 1, j - 1, 2))
                            xc(3) = fourth * (flowdoms(nn, groundlevel, sps)%x(il, i, j, 3) &
                                              + flowdoms(nn, groundlevel, sps)%x(il, i, j - 1, 3) &
                                              + flowdoms(nn, groundlevel, sps)%x(il, i - 1, j, 3) &
                                              + flowdoms(nn, groundlevel, sps)%x(il, i - 1, j - 1, 3))
 
                            ! Determine the coordinates relative to the center
                            ! of rotation.
 
                            xxc(1) = xc(1) - rotcenter(1)
                            xxc(2) = xc(2) - rotcenter(2)
                            xxc(3) = xc(3) - rotcenter(3)
 
                            ! Compute the velocity, which is the cross product
                            ! of rotRate and xc.
 
                            bcdata(mm)%uSlip(i, j, 1) = rotrate(2) * xxc(3) - rotrate(3) * xxc(2)
                            bcdata(mm)%uSlip(i, j, 2) = rotrate(3) * xxc(1) - rotrate(1) * xxc(3)
                            bcdata(mm)%uSlip(i, j, 3) = rotrate(1) * xxc(2) - rotrate(2) * xxc(1)
 
                            ! Determine the coordinates relative to the
                            ! rigid body rotation point.
 
                            xxc(1) = xc(1) - rotationpoint(1)
                            xxc(2) = xc(2) - rotationpoint(2)
                            xxc(3) = xc(3) - rotationpoint(3)
 
                            ! Determine the total velocity of the cell center.
                            ! This is a combination of rotation speed of this
                            ! block and the entire rigid body rotation.
 
                            bcdata(mm)%uSlip(i, j, 1) = bcdata(mm)%uSlip(i, j, 1) + velxgrid &
                                                        + derivrotationmatrix(1, 1) * xxc(1) &
                                                        + derivrotationmatrix(1, 2) * xxc(2) &
                                                        + derivrotationmatrix(1, 3) * xxc(3)
                            bcdata(mm)%uSlip(i, j, 2) = bcdata(mm)%uSlip(i, j, 2) + velygrid &
                                                        + derivrotationmatrix(2, 1) * xxc(1) &
                                                        + derivrotationmatrix(2, 2) * xxc(2) &
                                                        + derivrotationmatrix(2, 3) * xxc(3)
                            bcdata(mm)%uSlip(i, j, 3) = bcdata(mm)%uSlip(i, j, 3) + velzgrid &
                                                        + derivrotationmatrix(3, 1) * xxc(1) &
                                                        + derivrotationmatrix(3, 2) * xxc(2) &
                                                        + derivrotationmatrix(3, 3) * xxc(3)
                        end do
                    end do
 
                else if (bcfaceid(mm) == jmin) then
 
                    do j = bcdata(mm)%jcBeg, bcdata(mm)%jcEnd
                        do i = bcdata(mm)%icBeg, bcdata(mm)%icEnd
 
                            ! Compute the coordinates of the centroid of the face.
                            ! Normally this would be an average of i-1 and i, but
                            ! due to the usage of the pointer xFace and the fact
                            ! that x starts at index 0 this is shifted 1 index.
 
                            xc(1) = fourth * (flowdoms(nn, groundlevel, sps)%x(i, 1, j, 1) &
                                              + flowdoms(nn, groundlevel, sps)%x(i, 1, j - 1, 1) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, 1, j, 1) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, 1, j - 1, 1))
                            xc(2) = fourth * (flowdoms(nn, groundlevel, sps)%x(i, 1, j, 2) &
                                              + flowdoms(nn, groundlevel, sps)%x(i, 1, j - 1, 2) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, 1, j, 2) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, 1, j - 1, 2))
                            xc(3) = fourth * (flowdoms(nn, groundlevel, sps)%x(i, 1, j, 3) &
                                              + flowdoms(nn, groundlevel, sps)%x(i, 1, j - 1, 3) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, 1, j, 3) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, 1, j - 1, 3))
 
                            ! Determine the coordinates relative to the center
                            ! of rotation.
 
                            xxc(1) = xc(1) - rotcenter(1)
                            xxc(2) = xc(2) - rotcenter(2)
                            xxc(3) = xc(3) - rotcenter(3)
 
                            ! Compute the velocity, which is the cross product
                            ! of rotRate and xc.
 
                            bcdata(mm)%uSlip(i, j, 1) = rotrate(2) * xxc(3) - rotrate(3) * xxc(2)
                            bcdata(mm)%uSlip(i, j, 2) = rotrate(3) * xxc(1) - rotrate(1) * xxc(3)
                            bcdata(mm)%uSlip(i, j, 3) = rotrate(1) * xxc(2) - rotrate(2) * xxc(1)
 
                            ! Determine the coordinates relative to the
                            ! rigid body rotation point.
 
                            xxc(1) = xc(1) - rotationpoint(1)
                            xxc(2) = xc(2) - rotationpoint(2)
                            xxc(3) = xc(3) - rotationpoint(3)
 
                            ! Determine the total velocity of the cell center.
                            ! This is a combination of rotation speed of this
                            ! block and the entire rigid body rotation.
 
                            bcdata(mm)%uSlip(i, j, 1) = bcdata(mm)%uSlip(i, j, 1) + velxgrid &
                                                        + derivrotationmatrix(1, 1) * xxc(1) &
                                                        + derivrotationmatrix(1, 2) * xxc(2) &
                                                        + derivrotationmatrix(1, 3) * xxc(3)
                            bcdata(mm)%uSlip(i, j, 2) = bcdata(mm)%uSlip(i, j, 2) + velygrid &
                                                        + derivrotationmatrix(2, 1) * xxc(1) &
                                                        + derivrotationmatrix(2, 2) * xxc(2) &
                                                        + derivrotationmatrix(2, 3) * xxc(3)
                            bcdata(mm)%uSlip(i, j, 3) = bcdata(mm)%uSlip(i, j, 3) + velzgrid &
                                                        + derivrotationmatrix(3, 1) * xxc(1) &
                                                        + derivrotationmatrix(3, 2) * xxc(2) &
                                                        + derivrotationmatrix(3, 3) * xxc(3)
                        end do
                    end do
 
                else if (bcfaceid(mm) == jmax) then
 
                    do j = bcdata(mm)%jcBeg, bcdata(mm)%jcEnd
                        do i = bcdata(mm)%icBeg, bcdata(mm)%icEnd
 
                            ! Compute the coordinates of the centroid of the face.
                            ! Normally this would be an average of i-1 and i, but
                            ! due to the usage of the pointer xFace and the fact
                            ! that x starts at index 0 this is shifted 1 index.
 
                            xc(1) = fourth * (flowdoms(nn, groundlevel, sps)%x(i, jl, j, 1) &
                                              + flowdoms(nn, groundlevel, sps)%x(i, jl, j - 1, 1) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, jl, j, 1) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, jl, j - 1, 1))
                            xc(2) = fourth * (flowdoms(nn, groundlevel, sps)%x(i, jl, j, 2) &
                                              + flowdoms(nn, groundlevel, sps)%x(i, jl, j - 1, 2) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, jl, j, 2) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, jl, j - 1, 2))
                            xc(3) = fourth * (flowdoms(nn, groundlevel, sps)%x(i, jl, j, 3) &
                                              + flowdoms(nn, groundlevel, sps)%x(i, jl, j - 1, 3) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, jl, j, 3) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, jl, j - 1, 3))
 
                            ! Determine the coordinates relative to the center
                            ! of rotation.
 
                            xxc(1) = xc(1) - rotcenter(1)
                            xxc(2) = xc(2) - rotcenter(2)
                            xxc(3) = xc(3) - rotcenter(3)
 
                            ! Compute the velocity, which is the cross product
                            ! of rotRate and xc.
 
                            bcdata(mm)%uSlip(i, j, 1) = rotrate(2) * xxc(3) - rotrate(3) * xxc(2)
                            bcdata(mm)%uSlip(i, j, 2) = rotrate(3) * xxc(1) - rotrate(1) * xxc(3)
                            bcdata(mm)%uSlip(i, j, 3) = rotrate(1) * xxc(2) - rotrate(2) * xxc(1)
 
                            ! Determine the coordinates relative to the
                            ! rigid body rotation point.
 
                            xxc(1) = xc(1) - rotationpoint(1)
                            xxc(2) = xc(2) - rotationpoint(2)
                            xxc(3) = xc(3) - rotationpoint(3)
 
                            ! Determine the total velocity of the cell center.
                            ! This is a combination of rotation speed of this
                            ! block and the entire rigid body rotation.
 
                            bcdata(mm)%uSlip(i, j, 1) = bcdata(mm)%uSlip(i, j, 1) + velxgrid &
                                                        + derivrotationmatrix(1, 1) * xxc(1) &
                                                        + derivrotationmatrix(1, 2) * xxc(2) &
                                                        + derivrotationmatrix(1, 3) * xxc(3)
                            bcdata(mm)%uSlip(i, j, 2) = bcdata(mm)%uSlip(i, j, 2) + velygrid &
                                                        + derivrotationmatrix(2, 1) * xxc(1) &
                                                        + derivrotationmatrix(2, 2) * xxc(2) &
                                                        + derivrotationmatrix(2, 3) * xxc(3)
                            bcdata(mm)%uSlip(i, j, 3) = bcdata(mm)%uSlip(i, j, 3) + velzgrid &
                                                        + derivrotationmatrix(3, 1) * xxc(1) &
                                                        + derivrotationmatrix(3, 2) * xxc(2) &
                                                        + derivrotationmatrix(3, 3) * xxc(3)
                        end do
                    end do
 
                else if (bcfaceid(mm) == kmin) then
 
                    do j = bcdata(mm)%jcBeg, bcdata(mm)%jcEnd
                        do i = bcdata(mm)%icBeg, bcdata(mm)%icEnd
 
                            ! Compute the coordinates of the centroid of the face.
                            ! Normally this would be an average of i-1 and i, but
                            ! due to the usage of the pointer xFace and the fact
                            ! that x starts at index 0 this is shifted 1 index.
 
                            xc(1) = fourth * (flowdoms(nn, groundlevel, sps)%x(i, j, 1, 1) &
                                              + flowdoms(nn, groundlevel, sps)%x(i, j - 1, 1, 1) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, j, 1, 1) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, j - 1, 1, 1))
                            xc(2) = fourth * (flowdoms(nn, groundlevel, sps)%x(i, j, 1, 2) &
                                              + flowdoms(nn, groundlevel, sps)%x(i, j - 1, 1, 2) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, j, 1, 2) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, j - 1, 1, 2))
                            xc(3) = fourth * (flowdoms(nn, groundlevel, sps)%x(i, j, 1, 3) &
                                              + flowdoms(nn, groundlevel, sps)%x(i, j - 1, 1, 3) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, j, 1, 3) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, j - 1, 1, 3))
 
                            ! Determine the coordinates relative to the center
                            ! of rotation.
 
                            xxc(1) = xc(1) - rotcenter(1)
                            xxc(2) = xc(2) - rotcenter(2)
                            xxc(3) = xc(3) - rotcenter(3)
 
                            ! Compute the velocity, which is the cross product
                            ! of rotRate and xc.
 
                            bcdata(mm)%uSlip(i, j, 1) = rotrate(2) * xxc(3) - rotrate(3) * xxc(2)
                            bcdata(mm)%uSlip(i, j, 2) = rotrate(3) * xxc(1) - rotrate(1) * xxc(3)
                            bcdata(mm)%uSlip(i, j, 3) = rotrate(1) * xxc(2) - rotrate(2) * xxc(1)
 
                            ! Determine the coordinates relative to the
                            ! rigid body rotation point.
 
                            xxc(1) = xc(1) - rotationpoint(1)
                            xxc(2) = xc(2) - rotationpoint(2)
                            xxc(3) = xc(3) - rotationpoint(3)
 
                            ! Determine the total velocity of the cell center.
                            ! This is a combination of rotation speed of this
                            ! block and the entire rigid body rotation.
 
                            bcdata(mm)%uSlip(i, j, 1) = bcdata(mm)%uSlip(i, j, 1) + velxgrid &
                                                        + derivrotationmatrix(1, 1) * xxc(1) &
                                                        + derivrotationmatrix(1, 2) * xxc(2) &
                                                        + derivrotationmatrix(1, 3) * xxc(3)
                            bcdata(mm)%uSlip(i, j, 2) = bcdata(mm)%uSlip(i, j, 2) + velygrid &
                                                        + derivrotationmatrix(2, 1) * xxc(1) &
                                                        + derivrotationmatrix(2, 2) * xxc(2) &
                                                        + derivrotationmatrix(2, 3) * xxc(3)
                            bcdata(mm)%uSlip(i, j, 3) = bcdata(mm)%uSlip(i, j, 3) + velzgrid &
                                                        + derivrotationmatrix(3, 1) * xxc(1) &
                                                        + derivrotationmatrix(3, 2) * xxc(2) &
                                                        + derivrotationmatrix(3, 3) * xxc(3)
                        end do
                    end do
 
                else if (bcfaceid(mm) == kmax) then
 
                    do j = bcdata(mm)%jcBeg, bcdata(mm)%jcEnd
                        do i = bcdata(mm)%icBeg, bcdata(mm)%icEnd
 
                            ! Compute the coordinates of the centroid of the face.
                            ! Normally this would be an average of i-1 and i, but
                            ! due to the usage of the pointer xFace and the fact
                            ! that x starts at index 0 this is shifted 1 index.
 
                            xc(1) = fourth * (flowdoms(nn, groundlevel, sps)%x(i, j, kl, 1) &
                                              + flowdoms(nn, groundlevel, sps)%x(i, j - 1, kl, 1) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, j, kl, 1) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, j - 1, kl, 1))
                            xc(2) = fourth * (flowdoms(nn, groundlevel, sps)%x(i, j, kl, 2) &
                                              + flowdoms(nn, groundlevel, sps)%x(i, j - 1, kl, 2) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, j, kl, 2) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, j - 1, kl, 2))
                            xc(3) = fourth * (flowdoms(nn, groundlevel, sps)%x(i, j, kl, 3) &
                                              + flowdoms(nn, groundlevel, sps)%x(i, j - 1, kl, 3) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, j, kl, 3) &
                                              + flowdoms(nn, groundlevel, sps)%x(i - 1, j - 1, kl, 3))
 
                            ! Determine the coordinates relative to the center
                            ! of rotation.
 
                            xxc(1) = xc(1) - rotcenter(1)
                            xxc(2) = xc(2) - rotcenter(2)
                            xxc(3) = xc(3) - rotcenter(3)
 
                            ! Compute the velocity, which is the cross product
                            ! of rotRate and xc.
 
                            bcdata(mm)%uSlip(i, j, 1) = rotrate(2) * xxc(3) - rotrate(3) * xxc(2)
                            bcdata(mm)%uSlip(i, j, 2) = rotrate(3) * xxc(1) - rotrate(1) * xxc(3)
                            bcdata(mm)%uSlip(i, j, 3) = rotrate(1) * xxc(2) - rotrate(2) * xxc(1)
 
                            ! Determine the coordinates relative to the
                            ! rigid body rotation point.
 
                            xxc(1) = xc(1) - rotationpoint(1)
                            xxc(2) = xc(2) - rotationpoint(2)
                            xxc(3) = xc(3) - rotationpoint(3)
 
                            ! Determine the total velocity of the cell center.
                            ! This is a combination of rotation speed of this
                            ! block and the entire rigid body rotation.
 
                            bcdata(mm)%uSlip(i, j, 1) = bcdata(mm)%uSlip(i, j, 1) + velxgrid &
                                                        + derivrotationmatrix(1, 1) * xxc(1) &
                                                        + derivrotationmatrix(1, 2) * xxc(2) &
                                                        + derivrotationmatrix(1, 3) * xxc(3)
                            bcdata(mm)%uSlip(i, j, 2) = bcdata(mm)%uSlip(i, j, 2) + velygrid &
                                                        + derivrotationmatrix(2, 1) * xxc(1) &
                                                        + derivrotationmatrix(2, 2) * xxc(2) &
                                                        + derivrotationmatrix(2, 3) * xxc(3)
                            bcdata(mm)%uSlip(i, j, 3) = bcdata(mm)%uSlip(i, j, 3) + velzgrid &
                                                        + derivrotationmatrix(3, 1) * xxc(1) &
                                                        + derivrotationmatrix(3, 2) * xxc(2) &
                                                        + derivrotationmatrix(3, 3) * xxc(3)
                        end do
                    end do
 
                end if
 
            end do bocoloop2
 
        end if testuseoldcoor
 
    end subroutine slipvelocitiesfinelevel_block
 
#ifndef USE_TAPENADE
    subroutine normalvelocitiesalllevels(sps)
        !
        ! Shell function to call normalVelocities_block on all blocks/levels
        !
        use constants
        use blockpointers, only: ndom, flowdoms
        use inputtimespectral, only: ntimeintervalsspectral
        use iteration, only: groundlevel
        use utils, only: setpointers
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: sps
        !      Local variables.
        !
        integer(kind=intType) :: nn, level, nLevels
 
        nlevels = ubound(flowdoms, 2)
        levelloop: do level = groundlevel, nlevels
            domains: do nn = 1, ndom
 
                ! Set the pointers for this block.
 
                call setpointers(nn, level, sps)
 
                call normalvelocities_block(sps)
 
            end do domains
        end do levelloop
    end subroutine normalvelocitiesalllevels
#endif
 
    subroutine normalvelocities_block(sps)
        !
        !       normalVelocitiesAllLevels computes the normal grid
        !       velocities of some boundary faces of the moving blocks for
        !       spectral mode sps. All grid levels from ground level to the
        !       coarsest level are considered.
        !
        use constants
        use blockpointers, only: il, jl, kl, addgridvelocities, nbocos, bcdata, &
                                 sfacei, sfacej, sfacek, bcfaceid, si, sj, sk
        !use iteration
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: sps
        !
        !      Local variables.
        !
        integer(kind=intType) :: mm
        integer(kind=intType) :: i, j
 
        real(kind=realtype) :: weight, mult
 
        real(kind=realtype), dimension(:, :), pointer :: sface
        real(kind=realtype), dimension(:, :, :), pointer :: ss
 
        ! Check for a moving block. As it is possible that in a
        ! multidisicplinary environment additional grid velocities
        ! are set, the test should be done on addGridVelocities
        ! and not on blockIsMoving.
 
        testmoving: if (addgridvelocities) then
            !
            !             Determine the normal grid velocities of the boundaries.
            !             As these values are based on the unit normal. A division
            !             by the length of the normal is needed.
            !             Furthermore the boundary unit normals are per definition
            !             outward pointing, while on the iMin, jMin and kMin
            !             boundaries the face normals are inward pointing. This
            !             is taken into account by the factor mult.
            !
            ! Loop over the boundary subfaces.
 
            bocoloop: do mm = 1, nbocos
 
                ! Check whether rFace is allocated.
 
                testassoc: if (associated(bcdata(mm)%rFace)) then
 
                    ! Determine the block face on which the subface is
                    ! located and set some variables accordingly.
 
                    ! The new procedure is less elegant as the previous one.
                    ! But the new stands up to Tapenade.
                    if (bcfaceid(mm) == imin) then
 
                        mult = -one
 
                        do j = bcdata(mm)%jcBeg, bcdata(mm)%jcEnd
                            do i = bcdata(mm)%icBeg, bcdata(mm)%icEnd
 
                                ! Compute the inverse of the length of the normal
                                ! vector and possibly correct for inward pointing.
 
                                weight = sqrt(si(1, i, j, 1)**2 + si(1, i, j, 2)**2 &
                                              + si(1, i, j, 3)**2)
                                if (weight > zero) weight = mult / weight
 
                                ! Compute the normal velocity based on the outward
                                ! pointing unit normal.
 
                                bcdata(mm)%rFace(i, j) = weight * sfacei(1, i, j)
                            end do
                        end do
 
                    else if (bcfaceid(mm) == imax) then
 
                        mult = one
 
                        do j = bcdata(mm)%jcBeg, bcdata(mm)%jcEnd
                            do i = bcdata(mm)%icBeg, bcdata(mm)%icEnd
 
                                ! Compute the inverse of the length of the normal
                                ! vector and possibly correct for inward pointing.
 
                                weight = sqrt(si(il, i, j, 1)**2 + si(il, i, j, 2)**2 &
                                              + si(il, i, j, 3)**2)
                                if (weight > zero) weight = mult / weight
 
                                ! Compute the normal velocity based on the outward
                                ! pointing unit normal.
 
                                bcdata(mm)%rFace(i, j) = weight * sfacei(il, i, j)
                            end do
                        end do
 
                    else if (bcfaceid(mm) == jmin) then
 
                        mult = -one
 
                        do j = bcdata(mm)%jcBeg, bcdata(mm)%jcEnd
                            do i = bcdata(mm)%icBeg, bcdata(mm)%icEnd
 
                                ! Compute the inverse of the length of the normal
                                ! vector and possibly correct for inward pointing.
 
                                weight = sqrt(sj(i, 1, j, 1)**2 + sj(i, 1, j, 2)**2 &
                                              + sj(i, 1, j, 3)**2)
                                if (weight > zero) weight = mult / weight
 
                                ! Compute the normal velocity based on the outward
                                ! pointing unit normal.
 
                                bcdata(mm)%rFace(i, j) = weight * sfacej(i, 1, j)
                            end do
                        end do
 
                    else if (bcfaceid(mm) == jmax) then
 
                        mult = one
 
                        do j = bcdata(mm)%jcBeg, bcdata(mm)%jcEnd
                            do i = bcdata(mm)%icBeg, bcdata(mm)%icEnd
 
                                ! Compute the inverse of the length of the normal
                                ! vector and possibly correct for inward pointing.
 
                                weight = sqrt(sj(i, jl, j, 1)**2 + sj(i, jl, j, 2)**2 &
                                              + sj(i, jl, j, 3)**2)
                                if (weight > zero) weight = mult / weight
 
                                ! Compute the normal velocity based on the outward
                                ! pointing unit normal.
 
                                bcdata(mm)%rFace(i, j) = weight * sfacej(i, jl, j)
                            end do
                        end do
 
                    else if (bcfaceid(mm) == kmin) then
 
                        mult = -one
 
                        do j = bcdata(mm)%jcBeg, bcdata(mm)%jcEnd
                            do i = bcdata(mm)%icBeg, bcdata(mm)%icEnd
 
                                ! Compute the inverse of the length of the normal
                                ! vector and possibly correct for inward pointing.
 
                                weight = sqrt(sk(i, j, 1, 1)**2 + sk(i, j, 1, 2)**2 &
                                              + sk(i, j, 1, 3)**2)
                                if (weight > zero) weight = mult / weight
 
                                ! Compute the normal velocity based on the outward
                                ! pointing unit normal.
 
                                bcdata(mm)%rFace(i, j) = weight * sfacek(i, j, 1)
                            end do
                        end do
 
                    else if (bcfaceid(mm) == kmax) then
 
                        mult = one
 
                        do j = bcdata(mm)%jcBeg, bcdata(mm)%jcEnd
                            do i = bcdata(mm)%icBeg, bcdata(mm)%icEnd
 
                                ! Compute the inverse of the length of the normal
                                ! vector and possibly correct for inward pointing.
 
                                weight = sqrt(sk(i, j, kl, 1)**2 + sk(i, j, kl, 2)**2 &
                                              + sk(i, j, kl, 3)**2)
                                if (weight > zero) weight = mult / weight
 
                                ! Compute the normal velocity based on the outward
                                ! pointing unit normal.
 
                                bcdata(mm)%rFace(i, j) = weight * sfacek(i, j, kl)
                            end do
                        end do
 
                    end if
 
                end if testassoc
            end do bocoloop
 
            else testmoving
 
            ! Block is not moving. Loop over the boundary faces and set
            ! the normal grid velocity to zero if allocated.
 
            do mm = 1, nbocos
                if (associated(bcdata(mm)%rFace)) &
                    bcdata(mm)%rFace = zero
            end do
 
        end if testmoving
 
    end subroutine normalvelocities_block
 
    ! ----------------------------------------------------------------------
    !                                                                      |
    !                    No Tapenade Routine below this line               |
    !                                                                      |
    ! ----------------------------------------------------------------------
 
#ifndef USE_TAPENADE
 
    subroutine shiftsolution
        !
        !       shiftSolution shifts the solution of the older time levels,
        !       such that a new time step can be started.
        !
        use constants
        use blockpointers, only: il, jl, kl, nbkglobal, wold, w, ndom
        use cgnsgrid, only: cgnsdoms
        use flowvarrefstate, only: nw
        use iteration, only: groundlevel, noldlevels
        use inputtimespectral, only: ntimeintervalsspectral
        use inputunsteady, only: deltat
        use monitor, only: timeunsteadyrestart, timeunsteady
        use utils, only: setpointers, rotmatrixrigidbody
        implicit none
        !
        !      Local variables.
        !
        integer(kind=intType) :: i, j, k, l, sps, nn, mm, ll
 
        real(kind=realtype) :: told, tnew
        real(kind=realtype) :: vvx, vvy, vvz, vxi, veta, vzeta
        real(kind=realtype) :: t, anglex, angley, anglez
        real(kind=realtype) :: phi, cosphi, sinphi
        real(kind=realtype) :: xix, xiy, xiz, etax, etay, etaz
        real(kind=realtype) :: zetax, zetay, zetaz
 
        real(kind=realtype), dimension(3) :: rotationpoint
        real(kind=realtype), dimension(3, 3) :: rotationmatrix
 
        ! Compute the rotation matrix of the rigid body rotation as well
        ! as the rotation point; the latter is not needed to correct the
        ! velocities, but the routine rotMatrixRigidBody is also used
        ! for coordinates.
 
        tnew = timeunsteady + timeunsteadyrestart
        told = tnew - deltat
 
        call rotmatrixrigidbody(tnew, told, rotationmatrix, rotationpoint)
 
        ! Loop over the number of spectral solutions and local blocks.
        ! Although this routine is only called in unsteady mode where the
        ! number of spectral solutions is 1, this loop is there just for
        ! consistency.
 
        spectralloop: do sps = 1, ntimeintervalsspectral
            domains: do nn = 1, ndom
 
                ! Set the pointers for this block on the ground level.
 
                call setpointers(nn, groundlevel, sps)
 
                ! Shift the solution already stored in wOld.
 
                loopoldlevels: do mm = noldlevels, 2, -1
 
                    ! Shift the owned solution variables from level mm-1 to mm.
 
                    ll = mm - 1
 
                    do l = 1, nw
                        do k = 2, kl
                            do j = 2, jl
                                do i = 2, il
                                    wold(mm, i, j, k, l) = wold(ll, i, j, k, l)
                                end do
                            end do
                        end do
                    end do
 
                end do loopoldlevels
 
                ! Shift the current solution into the 1st level of wOld.
                ! Note that in wOld the conservative flow variables are stored,
                ! while in w the velocity components are stored and not
                ! the momentum. Therefore this must be corrected.
                ! Also the turbulent primitive variables are stored, but this
                ! is okay, because the quasi-linear form of the turbulent
                ! transport equations is solved and not the conservative one.
 
                do l = 1, nw
                    do k = 2, kl
                        do j = 2, jl
                            do i = 2, il
                                wold(1, i, j, k, l) = w(i, j, k, l)
                            end do
                        end do
                    end do
                end do
 
                ! Make sure that the momentum variables are stored in wOld.
 
                do k = 2, kl
                    do j = 2, jl
                        do i = 2, il
                            wold(1, i, j, k, ivx) = wold(1, i, j, k, ivx) * wold(1, i, j, k, irho)
                            wold(1, i, j, k, ivy) = wold(1, i, j, k, ivy) * wold(1, i, j, k, irho)
                            wold(1, i, j, k, ivz) = wold(1, i, j, k, ivz) * wold(1, i, j, k, irho)
                        end do
                    end do
                end do
 
                ! To improve the initial guess of the velocity field the
                ! velocity of rotating parts is rotated. First the rigid
                ! body motion.
 
                do k = 2, kl
                    do j = 2, jl
                        do i = 2, il
                            vvx = w(i, j, k, ivx)
                            vvy = w(i, j, k, ivy)
                            vvz = w(i, j, k, ivz)
 
                            w(i, j, k, ivx) = rotationmatrix(1, 1) * vvx &
                                              + rotationmatrix(1, 2) * vvy &
                                              + rotationmatrix(1, 3) * vvz
                            w(i, j, k, ivy) = rotationmatrix(2, 1) * vvx &
                                              + rotationmatrix(2, 2) * vvy &
                                              + rotationmatrix(2, 3) * vvz
                            w(i, j, k, ivz) = rotationmatrix(3, 1) * vvx &
                                              + rotationmatrix(3, 2) * vvy &
                                              + rotationmatrix(3, 3) * vvz
                        end do
                    end do
                end do
 
                ! Apply an additional correction for the velocity components
                ! if a rotation rate is prescribed for this block.
 
                rottest: if (cgnsdoms(nbkglobal)%rotatingFrameSpecified) then
 
                    ! Compute the rotation angles.
 
                    anglex = deltat * cgnsdoms(nbkglobal)%rotRate(1)
                    angley = deltat * cgnsdoms(nbkglobal)%rotRate(2)
                    anglez = deltat * cgnsdoms(nbkglobal)%rotRate(3)
 
                    ! Compute the unit vector in the direction of the rotation
                    ! axis, which will be called the xi-direction.
 
                    t = one / max(eps, sqrt(anglex**2 + angley**2 + anglez**2))
                    xix = t * anglex
                    xiy = t * angley
                    xiz = t * anglez
 
                    ! Determine the rotation angle in xi-direction and its sine
                    ! and cosine. Due to the definition of the xi-direction this
                    ! angle will always be positive.
 
                    phi = xix * anglex + xiy * angley + xiz * anglez
                    cosphi = cos(phi)
                    sinphi = sin(phi)
 
                    ! Loop over the cell centers.
 
                    do k = 2, kl
                        do j = 2, jl
                            do i = 2, il
 
                                ! Abbreviate the velocity components a bit easier.
 
                                vvx = w(i, j, k, ivx)
                                vvy = w(i, j, k, ivy)
                                vvz = w(i, j, k, ivz)
 
                                ! Determine the component of the velocity vector
                                ! in xi direction and determine the direction eta,
                                ! the direction of the velocity when the xi component
                                ! is substracted.
 
                                vxi = vvx * xix + vvy * xiy + vvz * xiz
 
                                etax = vvx - vxi * xix
                                etay = vvy - vxi * xiy
                                etaz = vvz - vxi * xiz
 
                                t = one / max(eps, sqrt(etax**2 + etay**2 + etaz**2))
                                etax = t * etax
                                etay = t * etay
                                etaz = t * etaz
 
                                ! Determine the velocity component in eta direction.
 
                                veta = vvx * etax + vvy * etay + vvz * etaz
 
                                ! Determine the unit vector in zeta-direction. This is
                                ! the cross product of the unit vectors in xi and in
                                ! eta-direction.
 
                                zetax = xiy * etaz - xiz * etay
                                zetay = xiz * etax - xix * etaz
                                zetaz = xix * etay - xiy * etax
 
                                ! Determine the velocity components in eta and zeta
                                ! direction after the rotation.
 
                                vzeta = veta * sinphi
                                veta = veta * cosphi
 
                                ! Compute the new Cartesian velocity components.
 
                                w(i, j, k, ivx) = vxi * xix + veta * etax + vzeta * zetax
                                w(i, j, k, ivy) = vxi * xiy + veta * etay + vzeta * zetay
                                w(i, j, k, ivz) = vxi * xiz + veta * etaz + vzeta * zetaz
 
                            end do
                        end do
                    end do
 
                end if rottest
 
            end do domains
        end do spectralloop
 
    end subroutine shiftsolution
 
    subroutine computeutau
        !
        ! Shell function to call computUTau  on all blocks
        !
        use constants
        use blockpointers, only: ndom
        use inputtimespectral, only: ntimeintervalsspectral
        use iteration, only: groundlevel
        use utils, only: setpointers
        implicit none
        !
        !      Local variables.
        !
        integer(kind=intType) :: sps, nn
 
        ! Loop over the number of spectral solutions.
 
        spectralloop: do sps = 1, ntimeintervalsspectral
 
            ! Loop over the number of blocks.
 
            domains: do nn = 1, ndom
 
                ! Set the pointers for this block.
 
                call setpointers(nn, groundlevel, sps)
 
                call computeutau_block
 
            end do domains
 
        end do spectralloop
 
    end subroutine computeutau
 
    subroutine computeutau_block
        !
        !       computeUtau computes the skin friction velocity for the
        !       viscous subfaces. This data is only needed if wall functions
        !       are used.
        !
        use constants
        use blockpointers
        use inputphysics
        use inputtimespectral
        use iteration
        use turbcurvefits, only: curveupre
 
        implicit none
        !
        !      Local variables.
        !
        integer(kind=intType) :: mm, i, j
 
        real(kind=realtype) :: re, vvx, vvy, vvz, veln, veltmag
 
        real(kind=realtype), dimension(:, :, :), pointer :: ww, norm, uslip
        real(kind=realtype), dimension(:, :), pointer :: dd2wall, rrlv
        real(kind=realtype), dimension(:, :), pointer :: utau
 
        ! Return immediately if no wall functions must be used.
 
        if (.not. wallfunctions) return
 
        ! Loop over the viscous subfaces of this block.
 
        viscsubfaces: do mm = 1, nviscbocos
 
            ! Set a bunch of pointers depending on the face id to make
            ! a generic treatment possible.
 
            select case (bcfaceid(mm))
 
            case (imin)
                ww => w(2, 1:, 1:, :);
                dd2wall => d2wall(2, :, :); rrlv => rlv(2, 1:, 1:)
 
                !=========================================================
 
            case (imax)
                ww => w(il, 1:, 1:, :)
                dd2wall => d2wall(il, :, :); rrlv => rlv(il, 1:, 1:)
 
                !=========================================================
 
            case (jmin)
                ww => w(1:, 2, 1:, :)
                dd2wall => d2wall(:, 2, :); rrlv => rlv(1:, 2, 1:)
 
                !=========================================================
 
            case (jmax)
                ww => w(1:, jl, 1:, :)
                dd2wall => d2wall(:, jl, :); rrlv => rlv(1:, jl, 1:)
 
                !=========================================================
 
            case (kmin)
                ww => w(1:, 1:, 2, :)
                dd2wall => d2wall(:, :, 2); rrlv => rlv(1:, 1:, 2)
 
                !=========================================================
 
            case (kmax)
                ww => w(1:, 1:, kl, :)
                dd2wall => d2wall(:, :, kl); rrlv => rlv(1:, 1:, kl)
 
            end select
 
            ! Set the pointers for the unit outward normals, uSlip
            ! and utau to make the code more readable.
 
            norm => bcdata(mm)%norm
            uslip => bcdata(mm)%uSlip
            utau => viscsubface(mm)%utau
 
            ! Loop over the quadrilateral faces of the subface. Note
            ! that the nodal range of BCData must be used and not the
            ! cell range, because the latter may include the halo's in i
            ! and j-direction. The offset +1 is there, because inBeg and
            ! jnBeg refer to nodal ranges and not to cell ranges.
            ! Note that an offset of -1 must be used in dd2Wall, because
            ! the original array d2Wall starts at 2.
 
            do j = (bcdata(mm)%jnBeg + 1), bcdata(mm)%jnEnd
                do i = (bcdata(mm)%inBeg + 1), bcdata(mm)%inEnd
 
                    ! Compute the velocity difference between the internal
                    ! cell and the wall.
 
                    vvx = ww(i, j, ivx) - uslip(i, j, 1)
                    vvy = ww(i, j, ivy) - uslip(i, j, 2)
                    vvz = ww(i, j, ivz) - uslip(i, j, 3)
 
                    ! Compute the normal velocity of the internal cell.
 
                    veln = vvx * norm(i, j, 1) + vvy * norm(i, j, 2) + vvz * norm(i, j, 3)
 
                    ! Compute the magnitude of the tangential velocity.
 
                    veltmag = max(eps, sqrt(vvx * vvx + vvy * vvy + vvz * vvz - veln * veln))
 
                    ! Compute the Reynolds number. Note that an offset of -1
                    ! must be used in dd2Wall, because the original array
                    ! d2Wall starts at 2.
                    ! Afterwards compute utau.
 
                    re = ww(i, j, irho) * veltmag * dd2wall(i - 1, j - 1) / rrlv(i, j)
                    utau(i, j) = veltmag / max(curveupre(re), eps)
 
                end do
            end do
 
        end do viscsubfaces
 
    end subroutine computeutau_block
 
    subroutine gridvelocitiesfinelevelpart1(useOldCoor, t, sps)
        !
        ! Shell function to call gridVelocitiesFineLevel on all blocks
        !
        use blockpointers
        use constants
        use inputtimespectral
        use iteration
        use utils, only: setpointers
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: sps
        logical, intent(in) :: useOldCoor
        real(kind=realtype), dimension(*), intent(in) :: t  !
        !      Local variables.
        !
        integer(kind=intType) :: nn
 
        ! Loop over the number of blocks.
 
        domains: do nn = 1, ndom
 
            ! Set the pointers for this block.
 
            call setpointers(nn, groundlevel, sps)
            call gridvelocitiesfinelevelpart1_block(useoldcoor, t, sps)
 
        end do domains
 
    end subroutine gridvelocitiesfinelevelpart1
 
    subroutine gridvelocitiesfinelevelpart1_block(useOldCoor, t, sps)
        !
        !       gridVelocitiesFineLevel computes the grid velocities for
        !       the cell centers and the normal grid velocities for the faces
        !       of moving blocks for the currently finest grid, i.e.
        !       groundLevel. The velocities are computed at time t for
        !       spectral mode sps. If useOldCoor is .true. the velocities
        !       are determined using the unsteady time integrator in
        !       combination with the old coordinates; otherwise the analytic
        !       form is used.
        !       Now it is split up into two parts.
        !       First part calculate the grid velocity using FIRST order BDF.
        !       Second part calculate the surface normal and normal velocity.
        !
        use blockpointers
        use cgnsgrid
        use flowvarrefstate
        use inputmotion
        use inputunsteady
        use iteration
        use inputphysics
        use inputtsstabderiv
        use monitor
        use communication
        use utils, only: tsalpha, tsbeta, tsmach, terminate, rotmatrixrigidbody, &
                         setcoeftimeintegrator, getdirangle
        use flowutils, only: derivativerotmatrixrigid, getdirvector
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: sps
        logical, intent(in) :: useOldCoor
 
        real(kind=realtype), dimension(*), intent(in) :: t
        !
        !      Local variables.
        !
        integer(kind=intType) :: nn, mm
        integer(kind=intType) :: i, j, k, ii, iie, jje, kke
 
        real(kind=realtype) :: oneover4dt, oneover8dt
        real(kind=realtype) :: velxgrid, velygrid, velzgrid, ainf
        real(kind=realtype) :: velxgrid0, velygrid0, velzgrid0
 
        real(kind=realtype), dimension(3) :: sc, xc, xxc
        real(kind=realtype), dimension(3) :: rotcenter, rotrate
 
        real(kind=realtype), dimension(3) :: rotationpoint
        real(kind=realtype), dimension(3, 3) :: rotationmatrix, &
                                                derivrotationmatrix
 
        real(kind=realtype) :: tnew, told
        real(kind=realtype), dimension(:, :), pointer :: sface
        real(kind=realtype), dimension(:, :, :), pointer :: svelo
 
        real(kind=realtype), dimension(:, :, :), pointer :: xx, ss
        real(kind=realtype), dimension(:, :, :, :), pointer :: xxold
 
        real(kind=realtype) :: intervalmach, alphats, alphaincrement, &
                               betats, betaincrement
        real(kind=realtype), dimension(3) :: veldir
        real(kind=realtype), dimension(3) :: refdirection
 
        ! Compute the mesh velocity from the given mesh Mach number.
 
        ! vel{x,y,z}Grid0 is the ACTUAL velocity you want at the
        ! geometry.
        ainf = sqrt(gammainf * pinf / rhoinf)
        velxgrid0 = (ainf * machgrid) * (-veldirfreestream(1))
        velygrid0 = (ainf * machgrid) * (-veldirfreestream(2))
        velzgrid0 = (ainf * machgrid) * (-veldirfreestream(3))
 
        ! Compute the derivative of the rotation matrix and the rotation
        ! point; needed for velocity due to the rigid body rotation of
        ! the entire grid. It is assumed that the rigid body motion of
        ! the grid is only specified if there is only 1 section present.
 
        call derivativerotmatrixrigid(derivrotationmatrix, rotationpoint, t(1))
 
        !compute the rotation matrix to update the velocities for the time
        !spectral stability derivative case...
 
        if (tsstability) then
            ! Determine the time values of the old and new time level.
            ! It is assumed that the rigid body rotation of the mesh is only
            ! used when only 1 section is present.
 
            tnew = timeunsteady + timeunsteadyrestart
            told = tnew - t(1)
 
            if (tspmode .or. tsqmode .or. tsrmode) then
                ! Compute the rotation matrix of the rigid body rotation as
                ! well as the rotation point; the latter may vary in time due
                ! to rigid body translation.
 
                call rotmatrixrigidbody(tnew, told, rotationmatrix, rotationpoint)
 
                velxgrid0 = rotationmatrix(1, 1) * velxgrid0 &
                            + rotationmatrix(1, 2) * velygrid0 &
                            + rotationmatrix(1, 3) * velzgrid0
                velygrid0 = rotationmatrix(2, 1) * velxgrid0 &
                            + rotationmatrix(2, 2) * velygrid0 &
                            + rotationmatrix(2, 3) * velzgrid0
                velzgrid0 = rotationmatrix(3, 1) * velxgrid0 &
                            + rotationmatrix(3, 2) * velygrid0 &
                            + rotationmatrix(3, 3) * velzgrid0
 
            elseif (tsalphamode) then
                !Determine the alpha for this time instance
                alphaincrement = tsalpha(degreepolalpha, coefpolalpha, &
                                         degreefouralpha, omegafouralpha, &
                                         coscoeffouralpha, sincoeffouralpha, t(1))
 
                alphats = alpha + alphaincrement
                !Determine the grid velocity for this alpha
                refdirection(:) = zero
                refdirection(1) = one
                call getdirvector(refdirection, alphats, beta, veldir, liftindex)
 
                !do I need to update the lift direction and drag direction as well?
                !set the effictive grid velocity for this time interval
                velxgrid0 = (ainf * machgrid) * (-veldir(1))
                velygrid0 = (ainf * machgrid) * (-veldir(2))
                velzgrid0 = (ainf * machgrid) * (-veldir(3))
 
            elseif (tsbetamode) then
 
                !Determine the alpha for this time instance
                betaincrement = tsbeta(degreepolbeta, coefpolbeta, &
                                       degreefourbeta, omegafourbeta, &
                                       coscoeffourbeta, sincoeffourbeta, t(1))
 
                betats = beta + betaincrement
                !Determine the grid velocity for this alpha
                refdirection(:) = zero
                refdirection(1) = one
                call getdirvector(refdirection, alpha, betats, veldir, liftindex)
 
                !do I need to update the lift direction and drag direction as well?
                !set the effictive grid velocity for this time interval
                velxgrid0 = (ainf * machgrid) * (-veldir(1))
                velygrid0 = (ainf * machgrid) * (-veldir(2))
                velzgrid0 = (ainf * machgrid) * (-veldir(3))
            elseif (tsmachmode) then
                !determine the mach number at this time interval
                intervalmach = tsmach(degreepolmach, coefpolmach, &
                                      degreefourmach, omegafourmach, &
                                      coscoeffourmach, sincoeffourmach, t(1))
                !set the effective grid velocity
                velxgrid0 = (ainf * (intervalmach + machgrid)) * (-veldirfreestream(1))
                velygrid0 = (ainf * (intervalmach + machgrid)) * (-veldirfreestream(2))
                velzgrid0 = (ainf * (intervalmach + machgrid)) * (-veldirfreestream(3))
 
            elseif (tsaltitudemode) then
                call terminate('gridVelocityFineLevel', 'altitude motion not yet implemented...')
            else
                call terminate('gridVelocityFineLevel', 'Not a recognized Stability Motion')
            end if
        end if
 
        testmoving: if (blockismoving) then
            ! REMOVED the rigid body rotation part for simplicity
 
            !
            !             The velocities must be determined via a finite
            !             difference formula using the coordinates of the old
            !             levels.
            !
            ! Set the coefficients for the time integrator and store
            ! the inverse of the physical nonDimensional time step,
            ! divided by 4 and 8, a bit easier.
 
            call setcoeftimeintegrator
            oneover4dt = fourth * timeref / deltat
            oneover8dt = half * oneover4dt
            !
            !             Grid velocities of the cell centers, including the
            !             1st level halo cells.
            !
            ! Loop over the cells, including the 1st level halo's.
 
            do k = 1, ke
                do j = 1, je
                    do i = 1, ie
 
                        ! Using FIRST order BDF for all cases
                        ! Refer to eq. 11b, found paper by C.Farhat http://dx.doi.org/10.1016/S0021-9991(03)00311-5
                        ! Same applies for the velocities of the faces below. Theta(n+1) = 1, Theta(n) = -1 therfore
                        ! it becoms a first order scheme.
 
                        ! The velocity of the cell center is determined
                        ! by a finite difference formula. First store
                        ! the current coordinate, multiplied by 8 and
                        ! coefTime(0) in sc.
 
                        sc(1) = (x(i - 1, j - 1, k - 1, 1) + x(i, j - 1, k - 1, 1) &
                                 + x(i - 1, j, k - 1, 1) + x(i, j, k - 1, 1) &
                                 + x(i - 1, j - 1, k, 1) + x(i, j - 1, k, 1) &
                                 + x(i - 1, j, k, 1) + x(i, j, k, 1))
                        sc(2) = (x(i - 1, j - 1, k - 1, 2) + x(i, j - 1, k - 1, 2) &
                                 + x(i - 1, j, k - 1, 2) + x(i, j, k - 1, 2) &
                                 + x(i - 1, j - 1, k, 2) + x(i, j - 1, k, 2) &
                                 + x(i - 1, j, k, 2) + x(i, j, k, 2))
                        sc(3) = (x(i - 1, j - 1, k - 1, 3) + x(i, j - 1, k - 1, 3) &
                                 + x(i - 1, j, k - 1, 3) + x(i, j, k - 1, 3) &
                                 + x(i - 1, j - 1, k, 3) + x(i, j - 1, k, 3) &
                                 + x(i - 1, j, k, 3) + x(i, j, k, 3))
 
                        ! Loop over the older levels to complete the
                        ! finite difference formula.
 
                        ii = 1 ! There was a loop over all old levels
                        sc(1) = sc(1) + (xold(ii, i - 1, j - 1, k - 1, 1) &
                                         + xold(ii, i, j - 1, k - 1, 1) &
                                         + xold(ii, i - 1, j, k - 1, 1) &
                                         + xold(ii, i, j, k - 1, 1) &
                                         + xold(ii, i - 1, j - 1, k, 1) &
                                         + xold(ii, i, j - 1, k, 1) &
                                         + xold(ii, i - 1, j, k, 1) &
                                         + xold(ii, i, j, k, 1)) &
                                * (-1.0_realtype)
                        sc(2) = sc(2) + (xold(ii, i - 1, j - 1, k - 1, 2) &
                                         + xold(ii, i, j - 1, k - 1, 2) &
                                         + xold(ii, i - 1, j, k - 1, 2) &
                                         + xold(ii, i, j, k - 1, 2) &
                                         + xold(ii, i - 1, j - 1, k, 2) &
                                         + xold(ii, i, j - 1, k, 2) &
                                         + xold(ii, i - 1, j, k, 2) &
                                         + xold(ii, i, j, k, 2)) &
                                * (-1.0_realtype)
                        sc(3) = sc(3) + (xold(ii, i - 1, j - 1, k - 1, 3) &
                                         + xold(ii, i, j - 1, k - 1, 3) &
                                         + xold(ii, i - 1, j, k - 1, 3) &
                                         + xold(ii, i, j, k - 1, 3) &
                                         + xold(ii, i - 1, j - 1, k, 3) &
                                         + xold(ii, i, j - 1, k, 3) &
                                         + xold(ii, i - 1, j, k, 3) &
                                         + xold(ii, i, j, k, 3)) &
                                * (-1.0_realtype)
 
                        ! Divide by 8 delta t to obtain the correct
                        ! velocities.
 
                        s(i, j, k, 1) = sc(1) * oneover8dt
                        s(i, j, k, 2) = sc(2) * oneover8dt
                        s(i, j, k, 3) = sc(3) * oneover8dt
                    end do
                end do
            end do
 
            !
            !             Velocities of the faces, vector.
            !
            ! Loop over the three directions.
 
            loopdir: do mm = 1, 3
 
                ! Set the upper boundaries depending on the direction.
 
                select case (mm)
                case (1_inttype)       ! normals in i-direction
                    iie = ie; jje = je; kke = ke
 
                case (2_inttype)       ! normals in j-direction
                    iie = je; jje = ie; kke = ke
 
                case (3_inttype)       ! normals in k-direction
                    iie = ke; jje = ie; kke = je
                end select
                !
                !               Face velocities in generalized i-direction.
                !               mm == 1: i-direction
                !               mm == 2: j-direction
                !               mm == 3: k-direction
                !
                do i = 0, iie
 
                    ! Set the pointers for the coordinates, normals and
                    ! normal velocities for this generalized i-plane.
                    ! This depends on the value of mm.
 
                    select case (mm)
                    case (1_inttype)       ! normals in i-direction
                        xx => x(i, :, :, :); xxold => xold(:, i, :, :, :)
                        svelo => sveloiale(i, :, :, :)
 
                    case (2_inttype)       ! normals in j-direction
                        xx => x(:, i, :, :); xxold => xold(:, :, i, :, :)
                        svelo => svelojale(:, i, :, :)
 
                    case (3_inttype)       ! normals in k-direction
                        xx => x(:, :, i, :); xxold => xold(:, :, :, i, :)
                        svelo => svelokale(:, :, i, :)
                    end select
 
                    ! Loop over the k and j-direction of this
                    ! generalized i-face. Note that due to the usage of
                    ! the pointers xx and xxOld an offset of +1 must be
                    ! used in the coordinate arrays, because x and xOld
                    ! originally start at 0 for the i, j and k indices.
                    ! print *, mm
                    do k = 1, kke
                        do j = 1, jje
 
                            ! The velocity of the face center is determined
                            ! by a finite difference formula. First store
                            ! the current coordinate, multiplied by 4 and
                            ! coefTime(0) in sc.
 
                            sc(1) = (xx(j + 1, k + 1, 1) + xx(j, k + 1, 1) &
                                     + xx(j + 1, k, 1) + xx(j, k, 1))
                            sc(2) = (xx(j + 1, k + 1, 2) + xx(j, k + 1, 2) &
                                     + xx(j + 1, k, 2) + xx(j, k, 2))
                            sc(3) = (xx(j + 1, k + 1, 3) + xx(j, k + 1, 3) &
                                     + xx(j + 1, k, 3) + xx(j, k, 3))
 
                            ii = 1 ! There was a loop who looped over nOldLevels
                            sc(1) = sc(1) + (xxold(ii, j + 1, k + 1, 1) &
                                             + xxold(ii, j, k + 1, 1) &
                                             + xxold(ii, j + 1, k, 1) &
                                             + xxold(ii, j, k, 1)) &
                                    * (-1.0_realtype)
                            sc(2) = sc(2) + (xxold(ii, j + 1, k + 1, 2) &
                                             + xxold(ii, j, k + 1, 2) &
                                             + xxold(ii, j + 1, k, 2) &
                                             + xxold(ii, j, k, 2)) &
                                    * (-1.0_realtype)
                            sc(3) = sc(3) + (xxold(ii, j + 1, k + 1, 3) &
                                             + xxold(ii, j, k + 1, 3) &
                                             + xxold(ii, j + 1, k, 3) &
                                             + xxold(ii, j, k, 3)) &
                                    * (-1.0_realtype)
 
                            ! Determine the dot product of sc and the normal
                            ! and divide by 4 deltaT to obtain the correct
                            ! value of the normal velocity.
 
                            svelo(j, k, 1) = sc(1) * oneover4dt
                            svelo(j, k, 2) = sc(2) * oneover4dt
                            svelo(j, k, 3) = sc(3) * oneover4dt
 
                            ! if ((i.ge.2) .and. (i.le.3) .and. (j.ge.2) .and. (j.le.3) .and. (k.ge.2) .and. (k.le.3)) then
                            !    print *, i,j,k, sVelo(j,k,:)
                            !    print *, '                                   ', xx(j,k,:)
                            !    print *, '                                   ', xxOld(1,j,k,:)
                            ! end if
 
                        end do
                    end do
                end do
 
            end do loopdir
 
        end if testmoving
 
    end subroutine gridvelocitiesfinelevelpart1_block
 
    !
    !      Here begins the second part
    !
 
    subroutine gridvelocitiesfinelevelpart2(useOldCoor, t, sps)
        !
        ! Shell function to call gridVelocitiesFineLevel on all blocks
        !
        use blockpointers
        use constants
        use inputtimespectral
        use iteration
        use utils, only: setpointers
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: sps
        logical, intent(in) :: useOldCoor
        real(kind=realtype), dimension(*), intent(in) :: t  !
        !      Local variables.
        !
        integer(kind=intType) :: nn
 
        ! Loop over the number of blocks.
 
        domains: do nn = 1, ndom
 
            ! Set the pointers for this block.
 
            call setpointers(nn, groundlevel, sps)
            call gridvelocitiesfinelevelpart2_block(useoldcoor, t, sps)
 
        end do domains
 
    end subroutine gridvelocitiesfinelevelpart2
 
    subroutine gridvelocitiesfinelevelpart2_block(useOldCoor, t, sps)
        !
        use blockpointers
        use cgnsgrid
        use flowvarrefstate
        use inputmotion
        use inputunsteady
        use iteration
        use inputphysics
        use inputtsstabderiv
        use monitor
        use communication
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: sps
        logical, intent(in) :: useOldCoor
        real(kind=realtype), dimension(*), intent(in) :: t
        !
        !      Local variables.
        !
        integer(kind=intType) :: nn, mm
        integer(kind=intType) :: i, j, k, ii, iie, jje, kke
 
        real(kind=realtype) :: oneover4dt, oneover8dt
        real(kind=realtype), dimension(3) :: sc, xc, xxc
        real(kind=realtype), dimension(:, :), pointer :: sface
        real(kind=realtype), dimension(:, :, :), pointer :: svelo
        real(kind=realtype), dimension(:, :, :), pointer :: xx, ss
        real(kind=realtype), dimension(:, :, :, :), pointer :: xxold
 
        testmoving: if (blockismoving) then
            !
            !             Normal grid velocities of the faces.
            !
            ! Loop over the three directions.
 
            loopdir: do mm = 1, 3
 
                ! Set the upper boundaries depending on the direction.
 
                select case (mm)
                case (1_inttype)       ! normals in i-direction
                    iie = ie; jje = je; kke = ke
 
                case (2_inttype)       ! normals in j-direction
                    iie = je; jje = ie; kke = ke
 
                case (3_inttype)       ! normals in k-direction
                    iie = ke; jje = ie; kke = je
                end select
                !
                !               Normal grid velocities in generalized i-direction.
                !               Mm == 1: i-direction
                !               mm == 2: j-direction
                !               mm == 3: k-direction
                !
                do i = 0, iie
 
                    ! Set the pointers for the coordinates, normals and
                    ! normal velocities for this generalized i-plane.
                    ! This depends on the value of mm.
 
                    select case (mm)
                    case (1_inttype)       ! normals in i-direction
                        ss => si(i, :, :, :); sface => sfacei(i, :, :)
                        svelo => sveloiale(i, :, :, :)
 
                    case (2_inttype)       ! normals in j-direction
                        ss => sj(:, i, :, :); sface => sfacej(:, i, :)
                        svelo => svelojale(:, i, :, :)
 
                    case (3_inttype)       ! normals in k-direction
                        ss => sk(:, :, i, :); sface => sfacek(:, :, i)
                        svelo => svelokale(:, :, i, :)
                    end select
 
                    ! Loop over the k and j-direction of this
                    ! generalized i-face. Note that due to the usage of
                    ! the pointers xx and xxOld an offset of +1 must be
                    ! used in the coordinate arrays, because x and xOld
                    ! originally start at 0 for the i, j and k indices.
 
                    do k = 1, kke
                        do j = 1, jje
 
                            ! Determine the dot product of sc and the normal
                            ! and divide by 4 deltaT to obtain the correct
                            ! value of the normal velocity.
 
                            sface(j, k) = svelo(j, k, 1) * ss(j, k, 1) &
                                          + svelo(j, k, 2) * ss(j, k, 2) &
                                          + svelo(j, k, 3) * ss(j, k, 3)
 
                        end do
                    end do
                end do
 
            end do loopdir
        end if testmoving
 
    end subroutine gridvelocitiesfinelevelpart2_block
 
    subroutine utauwf(rFilv)
        !
        !       utauWF substitutes the wall shear stress with values from a
        !       look-up table, if desired.
        !
        use constants
        use blockpointers, only: si, sj, sk, fw, rlv, d2wall, w, bcdata, viscsubface, &
                                 ie, je, ke, il, jl, kl, nviscbocos, bcfaceid
        use inputphysics, only: wallfunctions
        use turbcurvefits, only: curveupre
        implicit none
        !
        !      Subroutine argument.
        !
        real(kind=realtype), intent(in) :: rfilv
        !
        !      Local variables.
        !
        integer(kind=intType) :: i, j, nn
 
        real(kind=realtype) :: fact
        real(kind=realtype) :: tauxx, tauyy, tauzz
        real(kind=realtype) :: tauxy, tauxz, tauyz
        real(kind=realtype) :: rbar, ubar, vbar, wbar, vvx, vvy, vvz
        real(kind=realtype) :: fmx, fmy, fmz, frhoe
        real(kind=realtype) :: veln, velnx, velny, velnz, tx, ty, tz
        real(kind=realtype) :: veltx, velty, veltz, veltmag
        real(kind=realtype) :: txnx, txny, txnz, tynx, tyny, tynz
        real(kind=realtype) :: tznx, tzny, tznz
        real(kind=realtype) :: tautn, tauwall, utau, re
 
        real(kind=realtype), dimension(:, :, :), pointer :: ww1, ww2
        real(kind=realtype), dimension(:, :, :), pointer :: ss, rres
        real(kind=realtype), dimension(:, :, :), pointer :: norm
        real(kind=realtype), dimension(:, :), pointer :: rrlv2, dd2wall2
 
        ! Return immediately if no wall functions must be used.
 
        if (.not. wallfunctions) return
 
        ! Loop over the viscous subfaces of this block.
 
        viscsubfaces: do nn = 1, nviscbocos
 
            ! Set a bunch of variables depending on the face id to make
            ! a generic treatment possible.
 
            select case (bcfaceid(nn))
 
            case (imin)
                fact = -one
 
                ss => si(1, :, :, :); rres => fw(2, 1:, 1:, :)
                ww2 => w(2, 1:, 1:, :); ww1 => w(1, 1:, 1:, :)
                dd2wall2 => d2wall(2, :, :); rrlv2 => rlv(2, 1:, 1:)
 
                !===========================================================
 
            case (imax)
                fact = one
 
                ss => si(il, :, :, :); rres => fw(il, 1:, 1:, :)
                ww2 => w(il, 1:, 1:, :); ww1 => w(ie, 1:, 1:, :)
                dd2wall2 => d2wall(il, :, :); rrlv2 => rlv(il, 1:, 1:)
 
                !===========================================================
 
            case (jmin)
                fact = -one
 
                ss => sj(:, 1, :, :); rres => fw(1:, 2, 1:, :)
                ww2 => w(1:, 2, 1:, :); ww1 => w(1:, 1, 1:, :)
                dd2wall2 => d2wall(:, 2, :); rrlv2 => rlv(1:, 2, 1:)
 
                !===========================================================
 
            case (jmax)
                fact = one
 
                ss => sj(:, jl, :, :); rres => fw(1:, jl, 1:, :)
                ww2 => w(1:, jl, 1:, :); ww1 => w(1:, je, 1:, :)
                dd2wall2 => d2wall(:, jl, :); rrlv2 => rlv(1:, jl, 1:)
 
                !===========================================================
 
            case (kmin)
                fact = -one
 
                ss => sk(:, :, 1, :); rres => fw(1:, 1:, 2, :)
                ww2 => w(1:, 1:, 2, :); ww1 => w(1:, 1:, 1, :)
                dd2wall2 => d2wall(:, :, 2); rrlv2 => rlv(1:, 1:, 2)
 
                !===========================================================
 
            case (kmax)
                fact = one
 
                ss => sk(:, :, kl, :); rres => fw(1:, 1:, kl, :)
                ww2 => w(1:, 1:, kl, :); ww1 => w(1:, 1:, ke, :)
                dd2wall2 => d2wall(:, :, kl); rrlv2 => rlv(1:, 1:, kl)
 
            end select
 
            ! Set the pointer for the unit outward normals.
 
            norm => bcdata(nn)%norm
 
            ! Loop over the quadrilateral faces of the subface. Note
            ! that the nodal range of BCData must be used and not the
            ! cell range, because the latter may include the halo's in i
            ! and j-direction. The offset +1 is there, because inBeg and
            ! jnBeg refer to nodal ranges and not to cell ranges.
 
            do j = (bcdata(nn)%jnBeg + 1), bcdata(nn)%jnEnd
                do i = (bcdata(nn)%inBeg + 1), bcdata(nn)%inEnd
 
                    ! Store the viscous stress tensor a bit easier.
 
                    tauxx = viscsubface(nn)%tau(i, j, 1)
                    tauyy = viscsubface(nn)%tau(i, j, 2)
                    tauzz = viscsubface(nn)%tau(i, j, 3)
                    tauxy = viscsubface(nn)%tau(i, j, 4)
                    tauxz = viscsubface(nn)%tau(i, j, 5)
                    tauyz = viscsubface(nn)%tau(i, j, 6)
 
                    ! Compute the velocities at the wall face; these are only
                    ! non-zero for moving a block. Also compute the density,
                    ! which is needed to compute the wall shear stress via
                    ! wall functions.
 
                    rbar = half * (ww2(i, j, irho) + ww1(i, j, irho))
                    ubar = half * (ww2(i, j, ivx) + ww1(i, j, ivx))
                    vbar = half * (ww2(i, j, ivy) + ww1(i, j, ivy))
                    wbar = half * (ww2(i, j, ivz) + ww1(i, j, ivz))
 
                    ! Compute the velocity difference between the internal cell
                    ! and the wall.
 
                    vvx = ww2(i, j, ivx) - ubar
                    vvy = ww2(i, j, ivy) - vbar
                    vvz = ww2(i, j, ivz) - wbar
 
                    ! Compute the normal velocity of the internal cell.
 
                    veln = vvx * norm(i, j, 1) + vvy * norm(i, j, 2) + vvz * norm(i, j, 3)
                    velnx = veln * norm(i, j, 1)
                    velny = veln * norm(i, j, 2)
                    velnz = veln * norm(i, j, 3)
 
                    ! Compute the tangential velocity, its magnitude and its
                    ! unit vector of the internal cell.
 
                    veltx = vvx - velnx
                    velty = vvy - velny
                    veltz = vvz - velnz
 
                    veltmag = max(eps, sqrt(veltx**2 + velty**2 + veltz**2))
 
                    tx = veltx / veltmag
                    ty = velty / veltmag
                    tz = veltz / veltmag
 
                    ! Compute some coefficients needed for the transformation
                    ! between the cartesian frame and the frame defined by the
                    ! tangential direction (tx,ty,tz) and the normal direction.
                    ! The minus sign is present, because for this transformation
                    ! the normal direction should be inward pointing and norm
                    ! is outward pointing.
 
                    txnx = -tx * norm(i, j, 1)
                    txny = -tx * norm(i, j, 2)
                    txnz = -tx * norm(i, j, 3)
 
                    tynx = -ty * norm(i, j, 1)
                    tyny = -ty * norm(i, j, 2)
                    tynz = -ty * norm(i, j, 3)
 
                    tznx = -tz * norm(i, j, 1)
                    tzny = -tz * norm(i, j, 2)
                    tznz = -tz * norm(i, j, 3)
 
                    ! Compute the tn component of the wall shear stress
                    ! tensor. Normally this is the only nonzero shear
                    ! stress component in the t-n frame.
 
                    tautn = tauxx * txnx + tauyy * tyny + tauzz * tznz &
                            + tauxy * (txny + tynx) &
                            + tauxz * (txnz + tznx) &
                            + tauyz * (tynz + tzny)
 
                    ! Compute the Reynolds number using the velocity, density,
                    ! laminar viscosity and wall distance. Note that an offset
                    ! of -1 must be used in dd2Wall2, because the original array
                    ! d2Wall starts at 2.
 
                    re = ww2(i, j, irho) * veltmag * dd2wall2(i - 1, j - 1) / rrlv2(i, j)
 
                    ! Determine the friction velocity from the table and
                    ! compute the wall shear stress from it.
 
                    utau = veltmag / max(curveupre(re), eps)
                    tauwall = rbar * utau * utau
 
                    ! Compute the correction to the wall shear stress tautn and
                    ! transform this correction back to the cartesian frame.
                    ! Take rFilv into account, such that the correction to the
                    ! stress tensor is computed correctly.
 
                    tautn = rfilv * tauwall - tautn
 
                    tauxx = two * tautn * txnx
                    tauyy = two * tautn * tyny
                    tauzz = two * tautn * tznz
 
                    tauxy = tautn * (txny + tynx)
                    tauxz = tautn * (txnz + tznx)
                    tauyz = tautn * (tynz + tzny)
 
                    ! Compute the correction to the viscous flux at the wall.
 
                    fmx = tauxx * ss(i, j, 1) + tauxy * ss(i, j, 2) &
                          + tauxz * ss(i, j, 3)
                    fmy = tauxy * ss(i, j, 1) + tauyy * ss(i, j, 2) &
                          + tauyz * ss(i, j, 3)
                    fmz = tauxz * ss(i, j, 1) + tauyz * ss(i, j, 2) &
                          + tauzz * ss(i, j, 3)
                    frhoe = (ubar * tauxx + vbar * tauxy + wbar * tauxz) * ss(i, j, 1) &
                            + (ubar * tauxy + vbar * tauyy + wbar * tauyz) * ss(i, j, 2) &
                            + (ubar * tauxz + vbar * tauyz + wbar * tauzz) * ss(i, j, 3)
 
                    ! Add them to the residual. Note that now the factor rFilv
                    ! is already taken into account via tau. Fact is present to
                    ! take inward/outward pointing normals into account
 
                    rres(i, j, imx) = rres(i, j, imx) - fact * fmx
                    rres(i, j, imy) = rres(i, j, imy) - fact * fmy
                    rres(i, j, imz) = rres(i, j, imz) - fact * fmz
                    rres(i, j, irhoe) = rres(i, j, irhoe) - fact * frhoe
 
                    ! Store the friction velocity for later use.
 
                    viscsubface(nn)%utau(i, j) = utau
 
                    ! Also add the correction to the wall stress tensor.
 
                    viscsubface(nn)%tau(i, j, 1) = &
                        viscsubface(nn)%tau(i, j, 1) + tauxx
                    viscsubface(nn)%tau(i, j, 2) = &
                        viscsubface(nn)%tau(i, j, 2) + tauyy
                    viscsubface(nn)%tau(i, j, 3) = &
                        viscsubface(nn)%tau(i, j, 3) + tauzz
                    viscsubface(nn)%tau(i, j, 4) = &
                        viscsubface(nn)%tau(i, j, 4) + tauxy
                    viscsubface(nn)%tau(i, j, 5) = &
                        viscsubface(nn)%tau(i, j, 5) + tauxz
                    viscsubface(nn)%tau(i, j, 6) = &
                        viscsubface(nn)%tau(i, j, 6) + tauyz
                end do
            end do
 
        end do viscsubfaces
 
    end subroutine utauwf
 
#ifndef USE_TAPENADE
    subroutine slipvelocitiescoarselevels(sps)
        !
        !       slipVelocitiesCoarseLevels determines the slip velocities
        !       for the given spectral solution starting from the known
        !       velocities on the finer level.
        !
        use constants
        use blockpointers
        use iteration
        use utils, only: setpointers
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: sps
        !
        !      Local variables.
        !
        integer(kind=intType) :: i, j, iiMax, jjMax
        integer(kind=intType) :: if1, if2, jf1, jf2
        integer(kind=intType) :: nLevels, level, levm1, nn, mm
 
        integer(kind=intType), dimension(:, :), pointer :: iFine, jFine
 
        real(kind=realtype), dimension(:, :, :), pointer :: uslip
        real(kind=realtype), dimension(:, :, :), pointer :: uslipfine
 
        ! Determine the number of multigrid levels.
 
        nlevels = ubound(flowdoms, 2)
 
        ! Loop over coarser grid levels, where ground level is considered
        ! as the finest grid.
 
        levelloop: do level = (groundlevel + 1), nlevels
 
            ! Set levm1 for the finer level.
 
            levm1 = level - 1
 
            ! Loop over the number of local blocks.
 
            domains: do nn = 1, ndom
 
                ! Set the pointers to the coarse block.
 
                call setpointers(nn, level, sps)
 
                ! Loop over the number of viscous subfaces.
 
                bocoloop: do mm = 1, nviscbocos
 
                    ! Set the pointers for uSlip and uSlipFine to make the
                    ! code more readable.
 
                    uslip => bcdata(mm)%uSlip
                    uslipfine => flowdoms(nn, levm1, sps)%BCData(mm)%uSlip
 
                    ! Determine the grid face on which the subface is located
                    ! and set some variables accordingly.
 
                    select case (bcfaceid(mm))
 
                    case (imin, imax)
                        iimax = jl; jjmax = kl
                        ifine => mgjfine; jfine => mgkfine
 
                    case (jmin, jmax)
                        iimax = il; jjmax = kl
                        ifine => mgifine; jfine => mgkfine
 
                    case (kmin, kmax)
                        iimax = il; jjmax = jl
                        ifine => mgifine; jfine => mgjfine
 
                    end select
 
                    ! Loop over the number of faces of the viscous subface.
                    ! First in the generalized j-direction.
 
                    do j = bcdata(mm)%jcBeg, bcdata(mm)%jcEnd
 
                        ! Determine the two children in this direction.
                        ! Take care of the halo's, as this info is only
                        ! available for owned cells.
 
                        if (j < 2) then
                            jf1 = 1; jf2 = 1
                        else if (j > jjmax) then
                            jf1 = jfine(jjmax, 2) + 1; jf2 = jf1
                        else
                            jf1 = jfine(j, 1); jf2 = jfine(j, 2)
                        end if
 
                        ! Loop in the generalized i-direction.
 
                        do i = bcdata(mm)%icBeg, bcdata(mm)%icEnd
 
                            ! Determine the two children in this direction.
                            ! Same story as in j-direction.
 
                            if (i < 2) then
                                if1 = 1; if2 = 1
                            else if (i > iimax) then
                                if1 = ifine(iimax, 2) + 1; if2 = if1
                            else
                                if1 = ifine(i, 1); if2 = ifine(i, 2)
                            end if
 
                            ! Average the fine grid velocities to the
                            ! coarse grid velocities.
 
                            uslip(i, j, 1) = fourth * (uslipfine(if1, jf1, 1) &
                                                       + uslipfine(if2, jf1, 1) &
                                                       + uslipfine(if1, jf2, 1) &
                                                       + uslipfine(if2, jf2, 1))
 
                            uslip(i, j, 2) = fourth * (uslipfine(if1, jf1, 2) &
                                                       + uslipfine(if2, jf1, 2) &
                                                       + uslipfine(if1, jf2, 2) &
                                                       + uslipfine(if2, jf2, 2))
 
                            uslip(i, j, 3) = fourth * (uslipfine(if1, jf1, 3) &
                                                       + uslipfine(if2, jf1, 3) &
                                                       + uslipfine(if1, jf2, 3) &
                                                       + uslipfine(if2, jf2, 3))
                        end do
                    end do
 
                end do bocoloop
            end do domains
        end do levelloop
 
    end subroutine slipvelocitiescoarselevels
#endif
 
    subroutine gridvelocitiescoarselevels(sps)
        !
        !       gridVelocitiesCoarseLevels computes the grid velocities for
        !       the cell centers and the normal grid velocities for the faces
        !       of moving blocks on the coarser grid levels. GroundLevel is
        !       considered the fine grid level.
        !
        use constants
        use blockpointers
        use iteration
        use utils, only: setpointers
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: sps
        !
        !      Local variables.
        !
        integer(kind=intType) :: nLevels, level, levm1, nn, mm
        integer(kind=intType) :: i, j, k, iie, jje, kke
        integer(kind=intType) :: ii, ii1, jj, jj1, kk, kk1
 
        integer(kind=intType), dimension(:, :), pointer :: iFine, jFine
        integer(kind=intType), dimension(:, :), pointer :: kFine
 
        real(kind=realtype) :: jjweight, kkweight, weight
 
        real(kind=realtype), dimension(:), pointer :: jweight, kweight
 
        real(kind=realtype), dimension(:, :), pointer :: sfine, sface
        real(kind=realtype), dimension(:, :, :, :), pointer :: sf
 
        ! Loop over the number of coarse grid levels, starting at
        ! groundLevel+1,
 
        nlevels = ubound(flowdoms, 2)
        levelloop: do level = groundlevel + 1, nlevels
 
            ! Loop over the number of local blocks.
 
            domains: do nn = 1, ndom
 
                ! Set the pointers for this block.
 
                call setpointers(nn, level, sps)
 
                ! Check for a moving block. As it is possible that in a
                ! multidisicplinary environment additional grid velocities
                ! are set, the test should be done on addGridVelocities
                ! and not on blockIsMoving.
 
                testmoving: if (addgridvelocities) then
                    !
                    !             Grid velocities of the cell centers, including the 1st
                    !             level halo cells. These are determined by accumulating
                    !             the fine grid values. At the end the internal halo's are
                    !             communicated to obtain the correct values.
                    !
                    levm1 = level - 1
 
                    ! Set the pointer sf to the cell velocities on the fine mesh.
 
                    sf => flowdoms(nn, levm1, sps)%s
 
                    ! Loop over the cells, including the 1st level halo's.
                    ! The indices kk, kk1 contain the corresponding fine
                    ! grid indices in k-direction. Idem for jj, jj1, ii
                    ! and ii1.
 
                    do k = 1, ke
                        if (k == 1) then
                            kk = 1; kk1 = 1
                        else if (k == ke) then
                            kk = flowdoms(nn, levm1, sps)%ke; kk1 = kk
                        else
                            kk = mgkfine(k, 1); kk1 = mgkfine(k, 2)
                        end if
 
                        do j = 1, je
                            if (j == 1) then
                                jj = 1; jj1 = 1
                            else if (j == je) then
                                jj = flowdoms(nn, levm1, sps)%je; jj1 = jj
                            else
                                jj = mgjfine(j, 1); jj1 = mgjfine(j, 2)
                            end if
 
                            do i = 1, ie
                                if (i == 1) then
                                    ii = 1; ii1 = 1
                                else if (i == ie) then
                                    ii = flowdoms(nn, levm1, sps)%ie; ii1 = ii
                                else
                                    ii = mgifine(i, 1); ii1 = mgifine(i, 2)
                                end if
 
                                ! Determine the coarse grid velocity by
                                ! averaging the fine grid values.
 
                                s(i, j, k, 1) = (sf(ii1, jj1, kk1, 1) + sf(ii, jj1, kk1, 1) &
                                                 + sf(ii1, jj, kk1, 1) + sf(ii, jj, kk1, 1) &
                                                 + sf(ii1, jj1, kk, 1) + sf(ii, jj1, kk, 1) &
                                                 + sf(ii1, jj, kk, 1) + sf(ii, jj, kk, 1)) &
                                                * eighth
                                s(i, j, k, 2) = (sf(ii1, jj1, kk1, 2) + sf(ii, jj1, kk1, 2) &
                                                 + sf(ii1, jj, kk1, 2) + sf(ii, jj, kk1, 2) &
                                                 + sf(ii1, jj1, kk, 2) + sf(ii, jj1, kk, 2) &
                                                 + sf(ii1, jj, kk, 2) + sf(ii, jj, kk, 2)) &
                                                * eighth
                                s(i, j, k, 3) = (sf(ii1, jj1, kk1, 3) + sf(ii, jj1, kk1, 3) &
                                                 + sf(ii1, jj, kk1, 3) + sf(ii, jj, kk1, 3) &
                                                 + sf(ii1, jj1, kk, 3) + sf(ii, jj1, kk, 3) &
                                                 + sf(ii1, jj, kk, 3) + sf(ii, jj, kk, 3)) &
                                                * eighth
                            end do
                        end do
                    end do
                    !
                    !             Normal grid velocities of the faces.
                    !
                    ! Loop over the three directions.
 
                    loopcoarsedir: do mm = 1, 3
 
                        ! Set some values depending on the situation.
 
                        select case (mm)
 
                        case (1_inttype)       ! Normals in i-direction
                            iie = ie; jje = je; kke = ke
                            ifine => mgifine; jfine => mgjfine; kfine => mgkfine
                            jweight => mgjweight; kweight => mgkweight
 
                        case (2_inttype)       ! Normals in j-direction
                            iie = je; jje = ie; kke = ke
                            ifine => mgjfine; jfine => mgifine; kfine => mgkfine
                            jweight => mgiweight; kweight => mgkweight
 
                        case (3_inttype)       ! Normals in k-direction
                            iie = ke; jje = ie; kke = je
                            ifine => mgkfine; jfine => mgifine; kfine => mgjfine
                            jweight => mgiweight; kweight => mgjweight
 
                        end select
                        !
                        !               Normal grid velocities in generalized i-direction.
                        !               mm == 1: i-direction
                        !               mm == 2: j-direction
                        !               mm == 3: k-direction
                        !
                        do i = 0, iie
 
                            ! Determine the i-index of the corresponding plane on
                            ! the fine grid. Note that halo planes are not entirely
                            ! correct. This is not really a problem.
 
                            if (i < 2) then
                                ii = i
                            else if (i < iie) then
                                ii = ifine(i, 2)
                            else
                                ii = ifine(iie - 1, 2) + 1
                            end if
 
                            ! Set the pointers for sFine and sFace, which will
                            ! contain the mesh velocities for this particular
                            ! plane. The index depends on the value of mm.
 
                            select case (mm)
                            case (1_inttype)
                                sfine => flowdoms(nn, levm1, sps)%sFaceI(ii, :, :)
                                sface => sfacei(i, :, :)
                            case (2_inttype)
                                sfine => flowdoms(nn, levm1, sps)%sFaceJ(:, ii, :)
                                sface => sfacej(:, i, :)
                            case (3_inttype)
                                sfine => flowdoms(nn, levm1, sps)%sFaceK(:, :, ii)
                                sface => sfacek(:, :, i)
                            end select
 
                            ! Loop over the k and j faces for this general i-plane.
                            ! Again the halo's are not entirely correct. Kk, kk1,
                            ! jj and jj1 are the children in k and j-direction
                            ! respectively.
 
                            do k = 1, kke
                                if (k == 1) then
                                    kk = 1; kk1 = 1
                                    kkweight = kweight(2)
                                else if (k == kke) then
                                    kk = kfine(kke - 1, 2) + 1; kk1 = kk
                                    kkweight = kweight(kke - 1)
                                else
                                    kk = kfine(k, 1); kk1 = kfine(k, 2)
                                    kweight = kweight(k)
                                end if
 
                                do j = 1, jje
                                    if (j == 1) then
                                        jj = 1; jj1 = 1
                                        jjweight = jweight(2)
                                    else if (j == jje) then
                                        jj = jfine(jje - 1, 2) + 1; jj1 = jj
                                        jweight = jweight(jje - 1)
                                    else
                                        jj = jfine(j, 1); jj1 = jfine(j, 2)
                                        jjweight = jweight(j)
                                    end if
 
                                    ! Determine the coarse grid normal velocity.
                                    ! Take the averaging weight into account; for
                                    ! a normal coarsening this weight is 1.0.
 
                                    weight = kkweight * jjweight
                                    sface(j, k) = weight * (sfine(jj1, kk1) &
                                                            + sfine(jj, kk1) &
                                                            + sfine(jj1, kk) &
                                                            + sfine(jj, kk))
                                end do
                            end do
                        end do
 
                    end do loopcoarsedir
                end if testmoving
            end do domains
 
            ! Exchange the cell centered velocities.
 
            call exchangecellgridvelocities(level, sps)
 
        end do levelloop
 
    end subroutine gridvelocitiescoarselevels
 
    !      ==================================================================
 
    subroutine exchangecellgridvelocities(level, sps)
        !
        !       exchangeCellGridVelocities exchanges the grid velocities in
        !       the cell centers for the given grid level and spectral
        !       solution.
        !
        use constants
        use block
        use communication
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: level, sps
        !
        !      Local variables.
        !
        integer :: size, procID, ierr, index
        integer, dimension(mpi_status_size) :: mpiStatus
 
        integer(kind=intType) :: i, j, ii, jj
        integer(kind=intType) :: d1, i1, j1, k1, d2, i2, j2, k2
 
        real(kind=realtype) :: alp
        real(kind=realtype), dimension(3) :: vv
 
        !       The 1 to 1 communication.
        !
        ! Send the variables. The data is first copied into
        ! the send buffer after which the buffer is sent asap.
 
        ii = 1
        sends: do i = 1, commpatterncell_1st(level)%nProcSend
 
            ! Store the processor id and the size of the message
            ! a bit easier.
 
            procid = commpatterncell_1st(level)%sendProc(i)
            size = 3 * commpatterncell_1st(level)%nSend(i)
 
            ! Copy the data in the correct part of the send buffer.
 
            jj = ii
            do j = 1, commpatterncell_1st(level)%nSend(i)
 
                ! Store the block id and the indices of the donor
                ! a bit easier.
 
                d1 = commpatterncell_1st(level)%sendList(i)%block(j)
                i1 = commpatterncell_1st(level)%sendList(i)%indices(j, 1)
                j1 = commpatterncell_1st(level)%sendList(i)%indices(j, 2)
                k1 = commpatterncell_1st(level)%sendList(i)%indices(j, 3)
 
                ! Store the grid velocities in sendBuffer, if they
                ! are allocated. Otherwise they are simply zero.
 
                if (flowdoms(d1, level, sps)%addGridVelocities) then
                    sendbuffer(jj) = flowdoms(d1, level, sps)%s(i1, j1, k1, 1)
                    sendbuffer(jj + 1) = flowdoms(d1, level, sps)%s(i1, j1, k1, 2)
                    sendbuffer(jj + 2) = flowdoms(d1, level, sps)%s(i1, j1, k1, 3)
                else
                    sendbuffer(jj) = zero
                    sendbuffer(jj + 1) = zero
                    sendbuffer(jj + 2) = zero
                end if
 
                jj = jj + 3
            end do
 
            ! Send the data.
 
            call mpi_isend(sendbuffer(ii), size, adflow_real, procid, &
                           procid, adflow_comm_world, sendrequests(i), &
                           ierr)
 
            ! Set ii to jj for the next processor.
 
            ii = jj
        end do sends
 
        ! Post the nonblocking receives.
 
        ii = 1
        receives: do i = 1, commpatterncell_1st(level)%nProcRecv
 
            ! Store the processor id and the size of the message
            ! a bit easier.
 
            procid = commpatterncell_1st(level)%recvProc(i)
            size = 3 * commpatterncell_1st(level)%nRecv(i)
 
            ! Post the receive.
 
            call mpi_irecv(recvbuffer(ii), size, adflow_real, procid, &
                           myid, adflow_comm_world, recvrequests(i), ierr)
 
            ! And update ii.
 
            ii = ii + size
        end do receives
 
        ! Copy the local data.
 
        localcopy: do i = 1, internalcell_1st(level)%nCopy
 
            ! Store the block and the indices of the donor a bit easier.
 
            d1 = internalcell_1st(level)%donorBlock(i)
            i1 = internalcell_1st(level)%donorIndices(i, 1)
            j1 = internalcell_1st(level)%donorIndices(i, 2)
            k1 = internalcell_1st(level)%donorIndices(i, 3)
 
            ! Idem for the halo's.
 
            d2 = internalcell_1st(level)%haloBlock(i)
            i2 = internalcell_1st(level)%haloIndices(i, 1)
            j2 = internalcell_1st(level)%haloIndices(i, 2)
            k2 = internalcell_1st(level)%haloIndices(i, 3)
 
            ! Copy the grid velocities, if they are both allocated.
            ! Otherwise they are either set to zero or nothing is done.
 
            if (flowdoms(d2, level, sps)%addGridVelocities) then
                if (flowdoms(d1, level, sps)%addGridVelocities) then
                    flowdoms(d2, level, sps)%s(i2, j2, k2, 1) = &
                        flowdoms(d1, level, sps)%s(i1, j1, k1, 1)
                    flowdoms(d2, level, sps)%s(i2, j2, k2, 2) = &
                        flowdoms(d1, level, sps)%s(i1, j1, k1, 2)
                    flowdoms(d2, level, sps)%s(i2, j2, k2, 3) = &
                        flowdoms(d1, level, sps)%s(i1, j1, k1, 3)
                else
                    flowdoms(d2, level, sps)%s(i2, j2, k2, 1) = zero
                    flowdoms(d2, level, sps)%s(i2, j2, k2, 2) = zero
                    flowdoms(d2, level, sps)%s(i2, j2, k2, 3) = zero
                end if
            end if
 
        end do localcopy
 
        ! Correct the periodic halo's of the internal communication
        ! pattern, if present.
 
        if (internalcell_1st(level)%nPeriodic > 0) &
            call correctperiodicgridvel(level, sps, &
                                        internalcell_1st(level)%nPeriodic, &
                                        internalcell_1st(level)%periodicData)
 
        ! Complete the nonblocking receives in an arbitrary sequence and
        ! copy the variables from the buffer into the halo's.
 
        size = commpatterncell_1st(level)%nProcRecv
        completerecvs: do i = 1, commpatterncell_1st(level)%nProcRecv
 
            ! Complete any of the requests.
 
            call mpi_waitany(size, recvrequests, index, mpistatus, ierr)
 
            ! Copy the data just arrived in the halo's.
 
            ii = index
            jj = 3 * commpatterncell_1st(level)%nRecvCum(ii - 1)
            do j = 1, commpatterncell_1st(level)%nRecv(ii)
 
                ! Store the block and the indices of the halo a bit easier.
 
                d2 = commpatterncell_1st(level)%recvList(ii)%block(j)
                i2 = commpatterncell_1st(level)%recvList(ii)%indices(j, 1)
                j2 = commpatterncell_1st(level)%recvList(ii)%indices(j, 2)
                k2 = commpatterncell_1st(level)%recvList(ii)%indices(j, 3)
 
                ! Copy the grid velocities from recvBuffer if they are
                ! both allocated.
 
                if (flowdoms(d2, level, sps)%addGridVelocities) then
                    flowdoms(d2, level, sps)%s(i2, j2, k2, 1) = recvbuffer(jj + 1)
                    flowdoms(d2, level, sps)%s(i2, j2, k2, 2) = recvbuffer(jj + 2)
                    flowdoms(d2, level, sps)%s(i2, j2, k2, 3) = recvbuffer(jj + 3)
                end if
 
                jj = jj + 3
            end do
 
        end do completerecvs
 
        ! Correct the periodic halo's of the external communication
        ! pattern, if present.
 
        if (commpatterncell_1st(level)%nPeriodic > 0) &
            call correctperiodicgridvel(level, sps, &
                                        commpatterncell_1st(level)%nPeriodic, &
                                        commpatterncell_1st(level)%periodicData)
 
        ! Complete the nonblocking sends.
 
        size = commpatterncell_1st(level)%nProcSend
        do i = 1, commpatterncell_1st(level)%nProcSend
            call mpi_waitany(size, sendrequests, index, mpistatus, ierr)
        end do
 
    end subroutine exchangecellgridvelocities
 
    !      ==================================================================
 
    subroutine correctperiodicgridvel(level, sps, nPeriodic, &
                                      periodicData)
        !
        !       correctPeriodicGridVel applies the periodic transformation
        !       to the grid velocities of the cell halo's in periodicData.
        !
        use block
        use communication
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: level, sps, nPeriodic
        type(periodicdatatype), dimension(:), pointer :: periodicData
        !
        !      Local variables.
        !
        integer(kind=intType) :: nn, mm, ii, i, j, k
        real(kind=realtype) :: vx, vy, vz
 
        real(kind=realtype), dimension(3, 3) :: rotmatrix
 
        ! Loop over the number of periodic transformations.
 
        do nn = 1, nperiodic
 
            ! Store the rotation matrix a bit easier.
 
            rotmatrix = periodicdata(nn)%rotMatrix
 
            ! Loop over the number of halo cells for this transformation.
 
            do ii = 1, periodicdata(nn)%nHalos
 
                ! Store the block and the indices a bit easier.
 
                mm = periodicdata(nn)%block(ii)
                i = periodicdata(nn)%indices(ii, 1)
                j = periodicdata(nn)%indices(ii, 2)
                k = periodicdata(nn)%indices(ii, 3)
 
                ! Check if the grid velocities have been allocated.
 
                if (flowdoms(mm, level, sps)%addGridVelocities) then
 
                    ! Store the original velocities in vx, vy, vz.
 
                    vx = flowdoms(mm, level, sps)%s(i, j, k, 1)
                    vy = flowdoms(mm, level, sps)%s(i, j, k, 2)
                    vz = flowdoms(mm, level, sps)%s(i, j, k, 3)
 
                    ! Compute the new velocity vector.
 
                    flowdoms(mm, level, sps)%s(i, j, k, 1) = rotmatrix(1, 1) * vx &
                                                             + rotmatrix(1, 2) * vy &
                                                             + rotmatrix(1, 3) * vz
                    flowdoms(mm, level, sps)%s(i, j, k, 2) = rotmatrix(2, 1) * vx &
                                                             + rotmatrix(2, 2) * vy &
                                                             + rotmatrix(2, 3) * vz
                    flowdoms(mm, level, sps)%s(i, j, k, 3) = rotmatrix(3, 1) * vx &
                                                             + rotmatrix(3, 2) * vy &
                                                             + rotmatrix(3, 3) * vz
 
                end if
            end do
 
        end do
 
    end subroutine correctperiodicgridvel
 
    !      ==================================================================
 
#endif
end module solverutils