utils.F90

Go to the documentation of this file.

module utils
    implicit none
 
contains
 
    function char2str(charArray, n)
        use constants
        !
        ! some gymnastics to cast a char array to string
        !
        implicit none
        !
        !      Function arguments.
        !
        character, dimension(maxCGNSNameLen), intent(in) :: chararray
        integer(kind=intType), intent(in) :: n
        !
        !      Function type
        !
        character(len=n) :: char2str
        !
        !      Local variables.
        !
        integer(kind=intType) :: i
        do i = 1, n
            char2str(i:i) = chararray(i)
        end do
 
    end function char2str
 
    function tsbeta(degreePolBeta, coefPolBeta, &
                    degreeFourBeta, omegaFourBeta, &
                    cosCoefFourBeta, sinCoefFourBeta, t)
        !
        !       TSbeta computes the angle of attack for a given Time interval
        !       in a time spectral solution.
        !
        use constants
        use inputphysics, only: equationmode
        implicit none
        !
        !      Function type
        !
        real(kind=realtype) :: tsbeta
        !
        !      Function arguments.
        !
        integer(kind=intType), intent(in) :: degreepolbeta
        integer(kind=intType), intent(in) :: degreefourbeta
 
        real(kind=realtype), intent(in) :: omegafourbeta, t
 
        real(kind=realtype), dimension(0:*), intent(in) :: coefpolbeta
        real(kind=realtype), dimension(0:*), intent(in) :: coscoeffourbeta
        real(kind=realtype), dimension(*), intent(in) :: sincoeffourbeta
        !
        !      Local variables.
        !
        integer(kind=intType) :: nn
 
        real(kind=realtype) :: beta, val
 
        ! Return immediately if this is a steady computation.
 
        if (equationmode == steady) then
            tsbeta = zero
            return
        end if
 
        ! Compute the polynomial contribution. If no polynomial was
        ! specified, the value of index 0 is set to zero automatically.
 
        beta = coefpolbeta(0)
        do nn = 1, degreepolbeta
            beta = beta + coefpolbeta(nn) * (t**nn)
        end do
 
        ! Compute the fourier contribution. Again the cosine coefficient
        ! of index 0 is defaulted to zero if not specified.
 
        beta = beta + coscoeffourbeta(0)
        do nn = 1, degreefourbeta
            val = nn * omegafourbeta * t
            beta = beta + coscoeffourbeta(nn) * cos(val) &
                   + sincoeffourbeta(nn) * sin(val)
        end do
 
        ! Set TSBeta to phi.
 
        tsbeta = beta
 
    end function tsbeta
 
    function tsbetadot(degreePolBeta, coefPolBeta, &
                       degreeFourBeta, omegaFourBeta, &
                       cosCoefFourBeta, sinCoefFourBeta, t)
        !
        !       TSbeta computes the angle of attack for a given Time interval
        !       in a time spectral solution.
        !
        use constants
        use inputphysics, only: equationmode
        implicit none
        !
        !      Function type
        !
        real(kind=realtype) :: tsbetadot
        !
        !      Function arguments.
        !
        integer(kind=intType), intent(in) :: degreepolbeta
        integer(kind=intType), intent(in) :: degreefourbeta
 
        real(kind=realtype), intent(in) :: omegafourbeta, t
 
        real(kind=realtype), dimension(0:*), intent(in) :: coefpolbeta
        real(kind=realtype), dimension(0:*), intent(in) :: coscoeffourbeta
        real(kind=realtype), dimension(*), intent(in) :: sincoeffourbeta
        !
        !      Local variables.
        !
        integer(kind=intType) :: nn
 
        real(kind=realtype) :: betadot, val
 
        ! Return immediately if this is a steady computation.
 
        if (equationmode == steady) then
            tsbetadot = zero
            return
        end if
 
        ! Compute the polynomial contribution. If no polynomial was
        ! specified, the value of index 0 is set to zero automatically.
 
        betadot = zero
        do nn = 1, degreepolbeta
            betadot = betadot + nn * coefpolbeta(nn) * (t**(nn - 1))
        end do
 
        ! Compute the fourier contribution. Again the cosine coefficient
        ! of index 0 is defaulted to zero if not specified.
 
        do nn = 1, degreefourbeta
            val = nn * omegafourbeta
            betadot = betadot - val * coscoeffourbeta(nn) * sin(val * t) &
                      + val * sincoeffourbeta(nn) * cos(val * t)
        end do
 
        ! Set TSBeta to phi.
 
        tsbetadot = betadot
 
    end function tsbetadot
 
    function tsmach(degreePolMach, coefPolMach, &
                    degreeFourMach, omegaFourMach, &
                    cosCoefFourMach, sinCoefFourMach, t)
        !
        !       TSMach computes the Mach Number for a given time interval
        !       in a time spectral solution.
        !
        use constants
        use inputphysics, only: equationmode
        implicit none
        !
        !      Function type
        !
        real(kind=realtype) :: tsmach
        !
        !      Function arguments.
        !
        integer(kind=intType), intent(in) :: degreepolmach
        integer(kind=intType), intent(in) :: degreefourmach
 
        real(kind=realtype), intent(in) :: omegafourmach, t
 
        real(kind=realtype), dimension(0:*), intent(in) :: coefpolmach
        real(kind=realtype), dimension(0:*), intent(in) :: coscoeffourmach
        real(kind=realtype), dimension(*), intent(in) :: sincoeffourmach
        !
        !      Local variables.
        !
        integer(kind=intType) :: nn
 
        real(kind=realtype) :: intervalmach, val
 
        ! Return immediately if this is a steady computation.
 
        if (equationmode == steady) then
            tsmach = zero
            return
        end if
 
        ! Compute the polynomial contribution. If no polynomial was
        ! specified, the value of index 0 is set to zero automatically.
 
        intervalmach = coefpolmach(0)
        do nn = 1, degreepolmach
            intervalmach = intervalmach + coefpolmach(nn) * (t**nn)
        end do
 
        ! Compute the fourier contribution. Again the cosine coefficient
        ! of index 0 is defaulted to zero if not specified.
 
        intervalmach = intervalmach + coscoeffourmach(0)
        do nn = 1, degreefourmach
            val = nn * omegafourmach * t
            intervalmach = intervalmach + coscoeffourmach(nn) * cos(val) &
                           + sincoeffourmach(nn) * sin(val)
        end do
        print *, 'inTSMach', intervalmach, nn, val, t
        ! Set TSMach to phi.
 
        tsmach = intervalmach
 
    end function tsmach
 
    function tsmachdot(degreePolMach, coefPolMach, &
                       degreeFourMach, omegaFourMach, &
                       cosCoefFourMach, sinCoefFourMach, t)
        !
        !       TSmach computes the angle of attack for a given Time interval
        !       in a time spectral solution.
        !
        use constants
        use inputphysics, only: equationmode
        implicit none
        !
        !      Function type
        !
        real(kind=realtype) :: tsmachdot
        !
        !      Function arguments.
        !
        integer(kind=intType), intent(in) :: degreepolmach
        integer(kind=intType), intent(in) :: degreefourmach
 
        real(kind=realtype), intent(in) :: omegafourmach, t
 
        real(kind=realtype), dimension(0:*), intent(in) :: coefpolmach
        real(kind=realtype), dimension(0:*), intent(in) :: coscoeffourmach
        real(kind=realtype), dimension(*), intent(in) :: sincoeffourmach
        !
        !      Local variables.
        !
        integer(kind=intType) :: nn
 
        real(kind=realtype) :: machdot, val
 
        ! Return immediately if this is a steady computation.
 
        if (equationmode == steady) then
            tsmachdot = zero
            return
        end if
 
        ! Compute the polynomial contribution. If no polynomial was
        ! specified, the value of index 0 is set to zero automatically.
 
        machdot = zero
        do nn = 1, degreepolmach
            machdot = machdot + nn * coefpolmach(nn) * (t**(nn - 1))
        end do
 
        ! Compute the fourier contribution. Again the cosine coefficient
        ! of index 0 is defaulted to zero if not specified.
 
        do nn = 1, degreefourmach
            val = nn * omegafourmach
            machdot = machdot - val * coscoeffourmach(nn) * sin(val * t) &
                      + val * sincoeffourmach(nn) * cos(val * t)
        end do
 
        ! Set TSMach to phi.
 
        tsmachdot = machdot
 
    end function tsmachdot
 
    function tsalpha(degreePolAlpha, coefPolAlpha, &
                     degreeFourAlpha, omegaFourAlpha, &
                     cosCoefFourAlpha, sinCoefFourAlpha, t)
        !
        !       TSalpha computes the angle of attack for a given Time interval
        !       in a time spectral solution.
        !
        use constants
        use inputphysics, only: equationmode
        implicit none
        !
        !      Function type
        !
        real(kind=realtype) :: tsalpha
        !
        !      Function arguments.
        !
        integer(kind=intType), intent(in) :: degreepolalpha
        integer(kind=intType), intent(in) :: degreefouralpha
 
        real(kind=realtype), intent(in) :: omegafouralpha, t
 
        real(kind=realtype), dimension(0:*), intent(in) :: coefpolalpha
        real(kind=realtype), dimension(0:*), intent(in) :: coscoeffouralpha
        real(kind=realtype), dimension(*), intent(in) :: sincoeffouralpha
        !
        !      Local variables.
        !
        integer(kind=intType) :: nn
 
        real(kind=realtype) :: alpha, val
 
        ! Return immediately if this is a steady computation.
 
        if (equationmode == steady) then
            tsalpha = zero
            return
        end if
 
        ! Compute the polynomial contribution. If no polynomial was
        ! specified, the value of index 0 is set to zero automatically.
        alpha = coefpolalpha(0)
        do nn = 1, degreepolalpha
            alpha = alpha + coefpolalpha(nn) * (t**nn)
        end do
 
        ! Compute the fourier contribution. Again the cosine coefficient
        ! of index 0 is defaulted to zero if not specified.
 
        alpha = alpha + coscoeffouralpha(0)
        do nn = 1, degreefouralpha
            val = nn * omegafouralpha * t
            alpha = alpha + coscoeffouralpha(nn) * cos(val) &
                    + sincoeffouralpha(nn) * sin(val)
        end do
        !print *,'inTSalpha',alpha,nn,val,t
        ! Set TSAlpha to phi.
 
        tsalpha = alpha
 
    end function tsalpha
 
    function tsalphadot(degreePolAlpha, coefPolAlpha, &
                        degreeFourAlpha, omegaFourAlpha, &
                        cosCoefFourAlpha, sinCoefFourAlpha, t)
        !
        !       TSalpha computes the angle of attack for a given Time interval
        !       in a time spectral solution.
        !
        use constants
        use inputphysics, only: equationmode
        implicit none
        !
        !      Function type
        !
        real(kind=realtype) :: tsalphadot
        !
        !      Function arguments.
        !
        integer(kind=intType), intent(in) :: degreepolalpha
        integer(kind=intType), intent(in) :: degreefouralpha
 
        real(kind=realtype), intent(in) :: omegafouralpha, t
 
        real(kind=realtype), dimension(0:*), intent(in) :: coefpolalpha
        real(kind=realtype), dimension(0:*), intent(in) :: coscoeffouralpha
        real(kind=realtype), dimension(*), intent(in) :: sincoeffouralpha
        !
        !      Local variables.
        !
        integer(kind=intType) :: nn
 
        real(kind=realtype) :: alphadot, val
 
        ! Return immediately if this is a steady computation.
 
        if (equationmode == steady) then
            tsalphadot = zero
            return
        end if
 
        ! Compute the polynomial contribution. If no polynomial was
        ! specified, the value of index 0 is set to zero automatically.
 
        alphadot = zero
        do nn = 1, degreepolalpha
            alphadot = alphadot + nn * coefpolalpha(nn) * (t**(nn - 1))
        end do
 
        ! Compute the fourier contribution. Again the cosine coefficient
        ! of index 0 is defaulted to zero if not specified.
 
        do nn = 1, degreefouralpha
            val = nn * omegafouralpha
            alphadot = alphadot - val * coscoeffouralpha(nn) * sin(val * t) &
                       + val * sincoeffouralpha(nn) * cos(val * t)
        end do
 
        ! Set TSAlpha to phi.
 
        tsalphadot = alphadot
 
    end function tsalphadot
 
    function derivativerigidrotangle(degreePolRot, &
                                     coefPolRot, &
                                     degreeFourRot, &
                                     omegaFourRot, &
                                     cosCoefFourRot, &
                                     sinCoefFourRot, t)
        !
        !       derivativeRigidRotAngle computes the time derivative of the
        !       rigid body rotation angle at the given time for the given
        !       arguments. The angle is described by a combination of a
        !       polynomial and fourier series.
        !
        use constants
        use inputphysics, only: equationmode
        use flowvarrefstate, only: timeref
        implicit none
        !
        !      Function type
        !
        real(kind=realtype) :: derivativerigidrotangle
        !
        !      Function arguments.
        !
        integer(kind=intType), intent(in) :: degreepolrot
        integer(kind=intType), intent(in) :: degreefourrot
 
        real(kind=realtype), intent(in) :: omegafourrot, t
 
        real(kind=realtype), dimension(0:*), intent(in) :: coefpolrot
        real(kind=realtype), dimension(0:*), intent(in) :: coscoeffourrot
        real(kind=realtype), dimension(*), intent(in) :: sincoeffourrot
        !
        !      Local variables.
        !
        integer(kind=intType) :: nn
 
        real(kind=realtype) :: dphi, val
 
        ! Return immediately if this is a steady computation.
 
        if (equationmode == steady) then
            derivativerigidrotangle = zero
            return
        end if
 
        ! Compute the polynomial contribution.
 
        dphi = zero
        do nn = 1, degreepolrot
            dphi = dphi + nn * coefpolrot(nn) * (t**(nn - 1))
        end do
 
        ! Compute the fourier contribution.
 
        do nn = 1, degreefourrot
            val = nn * omegafourrot
            dphi = dphi - val * coscoeffourrot(nn) * sin(val * t)
            dphi = dphi + val * sincoeffourrot(nn) * cos(val * t)
        end do
 
        ! Set derivativeRigidRotAngle to dPhi. Multiply by timeRef
        ! to obtain the correct non-dimensional value.
 
        derivativerigidrotangle = timeref * dphi
 
    end function derivativerigidrotangle
 
    function mydim(x, y)
 
        use constants
 
        real(kind=realtype) x, y
        real(kind=realtype) :: mydim
 
        mydim = x - y
        if (mydim < 0.0) then
            mydim = 0.0
        end if
 
    end function mydim
 
    function getcorrectfork()
 
        use constants
        use flowvarrefstate, only: kpresent
        use iteration, only: currentlevel, groundlevel
        implicit none
 
        logical :: getcorrectfork
 
        if (kpresent .and. currentlevel <= groundlevel) then
            getcorrectfork = .true.
        else
            getcorrectfork = .false.
        end if
    end function getcorrectfork
    subroutine terminate(routineName, errorMessage)
        !
        !       terminate writes an error message to standard output and
        !       terminates the execution of the program.
        !
        use constants
        use communication, only: adflow_comm_world, myid
        implicit none
        !
        !      Subroutine arguments
        !
        character(len=*), intent(in) :: routineName
        character(len=*), intent(in) :: errorMessage
#ifndef USE_TAPENADE
 
        !
        !      Local parameter
        !
        integer, parameter :: maxCharLine = 55
        !
        !      Local variables
        !
        integer :: ierr, len, i2
        logical :: firstTime
 
        character(len=len_trim(errorMessage)) :: message
        character(len=8) :: integerString
 
        !
        ! Copy the errorMessage into message. It is not possible to work
        ! with errorMessage directly, because it is modified in this
        ! routine. Sometimes a constant string is passed to this routine
        ! and some compilers simply fail then.
 
        message = errormessage
 
        ! Print a nice error message. In case of a parallel executable
        ! also the processor id is printed.
 
        print "(a)", "#"
        print "(a)", "#--------------------------- !!! Error !!! &
             &----------------------------"
 
        write (integerstring, "(i8)") myid
        integerstring = adjustl(integerstring)
 
        print "(2a)", "#* Terminate called by processor ", &
            trim(integerstring)
 
        ! Write the header of the error message.
 
        print "(2a)", "#* Run-time error in procedure ", &
            trim(routinename)
 
        ! Loop to write the error message. If the message is too long it
        ! is split over several lines.
 
        firsttime = .true.
        do
            ! Determine the remaining error message to be written.
            ! If longer than the maximum number of characters allowed
            ! on a line, it is attempted to split the message.
 
            message = adjustl(message)
            len = len_trim(message)
            i2 = min(maxcharline, len)
 
            if (i2 < len) i2 = index(message(:i2), " ", .true.) - 1
            if (i2 < 0) i2 = index(message, " ") - 1
            if (i2 < 0) i2 = len
 
            ! Write this part of the error message. If it is the first
            ! line of the message some additional stuff is printed.
 
            if (firsttime) then
                print "(2a)", "#* Error message: ", &
                    trim(message(:i2))
                firsttime = .false.
            else
                print "(2a)", "#*                ", &
                    trim(message(:i2))
            end if
 
            ! Exit the loop if the entire message has been written.
 
            if (i2 == len) exit
 
            ! Adapt the string for the next part to be written.
 
            message = message(i2 + 1:)
 
        end do
 
        ! Write the trailing message.
 
        print "(a)", "#*"
        print "(a)", "#* Now exiting"
        print "(a)", "#------------------------------------------&
             &----------------------------"
        print "(a)", "#"
 
        ! Call abort and stop the program. This stop should be done in
        ! abort, but just to be sure.
 
        call mpi_abort(adflow_comm_world, 1, ierr)
        stop
 
#endif
 
    end subroutine terminate
 
    subroutine rotmatrixrigidbody(tNew, tOld, rotationMatrix, &
                                  rotationPoint)
        !
        !       rotMatrixRigidBody determines the rotation matrix and the
        !       rotation point to determine the coordinates of the new time
        !       level starting from the coordinates of the old time level.
        !
        use constants
        use inputmotion
        use flowvarrefstate, only: lref
        implicit none
        !
        !      Subroutine arguments.
        !
        real(kind=realtype), intent(in) :: tnew, told
 
        real(kind=realtype), dimension(3), intent(out) :: rotationpoint
        real(kind=realtype), dimension(3, 3), intent(out) :: rotationmatrix
        !
        !      Local variables.
        !
        integer(kind=intType) :: i, j
 
        real(kind=realtype) :: phi
        real(kind=realtype) :: cosx, cosy, cosz, sinx, siny, sinz
 
        real(kind=realtype), dimension(3, 3) :: mnew, mold
 
        ! Determine the rotation angle around the x-axis for the new
        ! time level and the corresponding values of the sine and cosine.
 
        phi = rigidrotangle(degreepolxrot, coefpolxrot, &
                            degreefourxrot, omegafourxrot, &
                            coscoeffourxrot, sincoeffourxrot, tnew)
        sinx = sin(phi)
        cosx = cos(phi)
 
        ! Idem for the y-axis.
 
        phi = rigidrotangle(degreepolyrot, coefpolyrot, &
                            degreefouryrot, omegafouryrot, &
                            coscoeffouryrot, sincoeffouryrot, tnew)
        siny = sin(phi)
        cosy = cos(phi)
 
        ! Idem for the z-axis.
 
        phi = rigidrotangle(degreepolzrot, coefpolzrot, &
                            degreefourzrot, omegafourzrot, &
                            coscoeffourzrot, sincoeffourzrot, tnew)
        sinz = sin(phi)
        cosz = cos(phi)
 
        ! Construct the transformation matrix at the new time level.
        ! It is assumed that the sequence of rotation is first around the
        ! x-axis then around the y-axis and finally around the z-axis.
 
        mnew(1, 1) = cosy * cosz
        mnew(2, 1) = cosy * sinz
        mnew(3, 1) = -siny
 
        mnew(1, 2) = sinx * siny * cosz - cosx * sinz
        mnew(2, 2) = sinx * siny * sinz + cosx * cosz
        mnew(3, 2) = sinx * cosy
 
        mnew(1, 3) = cosx * siny * cosz + sinx * sinz
        mnew(2, 3) = cosx * siny * sinz - sinx * cosz
        mnew(3, 3) = cosx * cosy
 
        ! Determine the rotation angle around the x-axis for the old
        ! time level and the corresponding values of the sine and cosine.
 
        phi = rigidrotangle(degreepolxrot, coefpolxrot, &
                            degreefourxrot, omegafourxrot, &
                            coscoeffourxrot, sincoeffourxrot, told)
        sinx = sin(phi)
        cosx = cos(phi)
 
        ! Idem for the y-axis.
 
        phi = rigidrotangle(degreepolyrot, coefpolyrot, &
                            degreefouryrot, omegafouryrot, &
                            coscoeffouryrot, sincoeffouryrot, told)
        siny = sin(phi)
        cosy = cos(phi)
 
        ! Idem for the z-axis.
 
        phi = rigidrotangle(degreepolzrot, coefpolzrot, &
                            degreefourzrot, omegafourzrot, &
                            coscoeffourzrot, sincoeffourzrot, told)
        sinz = sin(phi)
        cosz = cos(phi)
 
        ! Construct the transformation matrix at the old time level.
 
        mold(1, 1) = cosy * cosz
        mold(2, 1) = cosy * sinz
        mold(3, 1) = -siny
 
        mold(1, 2) = sinx * siny * cosz - cosx * sinz
        mold(2, 2) = sinx * siny * sinz + cosx * cosz
        mold(3, 2) = sinx * cosy
 
        mold(1, 3) = cosx * siny * cosz + sinx * sinz
        mold(2, 3) = cosx * siny * sinz - sinx * cosz
        mold(3, 3) = cosx * cosy
 
        ! Construct the transformation matrix between the new and the
        ! old time level. This is mNew*inverse(mOld). However the
        ! inverse of mOld is the transpose.
 
        do j = 1, 3
            do i = 1, 3
                rotationmatrix(i, j) = mnew(i, 1) * mold(j, 1) &
                                       + mnew(i, 2) * mold(j, 2) &
                                       + mnew(i, 3) * mold(j, 3)
            end do
        end do
 
        ! Determine the rotation point at the old time level; it is
        ! possible that this value changes due to translation of the grid.
 
        !  aInf = sqrt(gammaInf*pInf/rhoInf)
 
        !  rotationPoint(1) = LRef*rotPoint(1) &
        !                   + MachGrid(1)*aInf*tOld/timeRef
        !  rotationPoint(2) = LRef*rotPoint(2) &
        !                   + MachGrid(2)*aInf*tOld/timeRef
        !  rotationPoint(3) = LRef*rotPoint(3) &
        !                   + MachGrid(3)*aInf*tOld/timeRef
 
        rotationpoint(1) = lref * rotpoint(1)
        rotationpoint(2) = lref * rotpoint(2)
        rotationpoint(3) = lref * rotpoint(3)
 
    end subroutine rotmatrixrigidbody
 
    function secondderivativerigidrotangle(degreePolRot, &
                                           coefPolRot, &
                                           degreeFourRot, &
                                           omegaFourRot, &
                                           cosCoefFourRot, &
                                           sinCoefFourRot, t)
        !
        !       2ndderivativeRigidRotAngle computes the 2nd time derivative of
        !       the rigid body rotation angle at the given time for the given
        !       arguments. The angle is described by a combination of a
        !       polynomial and fourier series.
        !
        use constants
        use flowvarrefstate, only: timeref
        use inputphysics, only: equationmode
        implicit none
        !
        !      Function type
        !
        real(kind=realtype) :: secondderivativerigidrotangle
        !
        !      Function arguments.
        !
        integer(kind=intType), intent(in) :: degreepolrot
        integer(kind=intType), intent(in) :: degreefourrot
 
        real(kind=realtype), intent(in) :: omegafourrot, t
 
        real(kind=realtype), dimension(0:*), intent(in) :: coefpolrot
        real(kind=realtype), dimension(0:*), intent(in) :: coscoeffourrot
        real(kind=realtype), dimension(*), intent(in) :: sincoeffourrot
        !
        !      Local variables.
        !
        integer(kind=intType) :: nn
 
        real(kind=realtype) :: dphi, val
 
        ! Return immediately if this is a steady computation.
 
        if (equationmode == steady) then
            secondderivativerigidrotangle = zero
            return
        end if
 
        ! Compute the polynomial contribution.
 
        dphi = zero
        do nn = 2, degreepolrot
            dphi = dphi + (nn - 1) * nn * coefpolrot(nn) * (t**(nn - 2))
        end do
 
        ! Compute the fourier contribution.
 
        do nn = 1, degreefourrot
            val = nn * omegafourrot
            dphi = dphi - val**2 * sincoeffourrot(nn) * sin(val * t)
            dphi = dphi - val**2 * coscoeffourrot(nn) * cos(val * t)
        end do
 
        ! Set derivativeRigidRotAngle to dPhi. Multiply by timeRef
        ! to obtain the correct non-dimensional value.
 
        secondderivativerigidrotangle = timeref**2 * dphi
 
    end function secondderivativerigidrotangle
 
    function rigidrotangle(degreePolRot, coefPolRot, &
                           degreeFourRot, omegaFourRot, &
                           cosCoefFourRot, sinCoefFourRot, t)
        !
        !       rigidRotAngle computes the rigid body rotation angle at the
        !       given time for the given arguments. The angle is described by
        !       a combination of a polynomial and fourier series.
        !
        use constants
        use inputphysics, only: equationmode
        implicit none
        !
        !      Function type
        !
        real(kind=realtype) :: rigidrotangle
        !
        !      Function arguments.
        !
        integer(kind=intType), intent(in) :: degreepolrot
        integer(kind=intType), intent(in) :: degreefourrot
 
        real(kind=realtype), intent(in) :: omegafourrot, t
 
        real(kind=realtype), dimension(0:*), intent(in) :: coefpolrot
        real(kind=realtype), dimension(0:*), intent(in) :: coscoeffourrot
        real(kind=realtype), dimension(*), intent(in) :: sincoeffourrot
        !
        !      Local variables.
        !
        integer(kind=intType) :: nn
 
        real(kind=realtype) :: phi, val
 
        ! Return immediately if this is a steady computation.
 
        if (equationmode == steady) then
            rigidrotangle = zero
            return
        end if
 
        ! Compute the polynomial contribution. If no polynomial was
        ! specified, the value of index 0 is set to zero automatically.
 
        phi = coefpolrot(0)
        do nn = 1, degreepolrot
            phi = phi + coefpolrot(nn) * (t**nn)
        end do
 
        ! Compute the fourier contribution. Again the cosine coefficient
        ! of index 0 is defaulted to zero if not specified.
 
        phi = phi + coscoeffourrot(0)
        do nn = 1, degreefourrot
            val = nn * omegafourrot * t
            phi = phi + coscoeffourrot(nn) * cos(val) &
                  + sincoeffourrot(nn) * sin(val)
        end do
 
        ! Set rigidRotAngle to phi.
 
        rigidrotangle = phi
 
    end function rigidrotangle
 
    subroutine setbcpointers(nn, spatialPointers)
        !
        !       setBCPointers sets the pointers needed for the boundary
        !       condition treatment on a general face, such that the boundary
        !       routines are only implemented once instead of 6 times.
        !
        use constants
        use blockpointers, only: w, p, rlv, rev, gamma, x, d2wall, &
                                 si, sj, sk, s, globalcell, bcdata, nx, il, ie, ib, &
                                 ny, jl, je, jb, nz, kl, ke, kb, bcfaceid, &
                                 addgridvelocities, sfacei, sfacej, sfacek, addgridvelocities
        use bcpointers, only: ww0, ww1, ww2, ww3, pp0, pp1, pp2, pp3, &
                              rlv0, rlv1, rlv2, rlv3, rev0, rev1, rev2, rev3, &
                              gamma0, gamma1, gamma2, gamma3, gcp, xx, ss, ssi, ssj, ssk, dd2wall, &
                              sface, istart, iend, jstart, jend, isize, jsize
        use inputphysics, only: cpmodel, equations
        implicit none
 
        ! Subroutine arguments.
        integer(kind=intType), intent(in) :: nn
        logical, intent(in) :: spatialPointers
 
        ! Determine the sizes of each face and point to just the range we
        ! need on each face.
        istart = bcdata(nn)%icBeg
        iend = bcdata(nn)%icEnd
        jstart = bcdata(nn)%jcBeg
        jend = bcdata(nn)%jcEnd
 
        ! Set the size of the subface
        isize = iend - istart + 1
        jsize = jend - jstart + 1
 
        ! Determine the face id on which the subface is located and set
        ! the pointers accordinly.
 
        select case (bcfaceid(nn))
 
            !---------------------------------------------------------------------------
        case (imin)
 
            ww3 => w(3, 1:, 1:, :)
            ww2 => w(2, 1:, 1:, :)
            ww1 => w(1, 1:, 1:, :)
            ww0 => w(0, 1:, 1:, :)
 
            pp3 => p(3, 1:, 1:)
            pp2 => p(2, 1:, 1:)
            pp1 => p(1, 1:, 1:)
            pp0 => p(0, 1:, 1:)
 
            rlv3 => rlv(3, 1:, 1:)
            rlv2 => rlv(2, 1:, 1:)
            rlv1 => rlv(1, 1:, 1:)
            rlv0 => rlv(0, 1:, 1:)
 
            rev3 => rev(3, 1:, 1:)
            rev2 => rev(2, 1:, 1:)
            rev1 => rev(1, 1:, 1:)
            rev0 => rev(0, 1:, 1:)
 
            gamma3 => gamma(3, 1:, 1:)
            gamma2 => gamma(2, 1:, 1:)
            gamma1 => gamma(1, 1:, 1:)
            gamma0 => gamma(0, 1:, 1:)
 
            gcp => globalcell(2, 1:, 1:)
            !---------------------------------------------------------------------------
 
        case (imax)
 
            ww3 => w(nx, 1:, 1:, :)
            ww2 => w(il, 1:, 1:, :)
            ww1 => w(ie, 1:, 1:, :)
            ww0 => w(ib, 1:, 1:, :)
 
            pp3 => p(nx, 1:, 1:)
            pp2 => p(il, 1:, 1:)
            pp1 => p(ie, 1:, 1:)
            pp0 => p(ib, 1:, 1:)
 
            rlv3 => rlv(nx, 1:, 1:)
            rlv2 => rlv(il, 1:, 1:)
            rlv1 => rlv(ie, 1:, 1:)
            rlv0 => rlv(ib, 1:, 1:)
 
            rev3 => rev(nx, 1:, 1:)
            rev2 => rev(il, 1:, 1:)
            rev1 => rev(ie, 1:, 1:)
            rev0 => rev(ib, 1:, 1:)
 
            gamma3 => gamma(nx, 1:, 1:)
            gamma2 => gamma(il, 1:, 1:)
            gamma1 => gamma(ie, 1:, 1:)
            gamma0 => gamma(ib, 1:, 1:)
 
            gcp => globalcell(il, 1:, 1:)
            !---------------------------------------------------------------------------
 
        case (jmin)
 
            ww3 => w(1:, 3, 1:, :)
            ww2 => w(1:, 2, 1:, :)
            ww1 => w(1:, 1, 1:, :)
            ww0 => w(1:, 0, 1:, :)
 
            pp3 => p(1:, 3, 1:)
            pp2 => p(1:, 2, 1:)
            pp1 => p(1:, 1, 1:)
            pp0 => p(1:, 0, 1:)
 
            rlv3 => rlv(1:, 3, 1:)
            rlv2 => rlv(1:, 2, 1:)
            rlv1 => rlv(1:, 1, 1:)
            rlv0 => rlv(1:, 0, 1:)
 
            rev3 => rev(1:, 3, 1:)
            rev2 => rev(1:, 2, 1:)
            rev1 => rev(1:, 1, 1:)
            rev0 => rev(1:, 0, 1:)
 
            gamma3 => gamma(1:, 3, 1:)
            gamma2 => gamma(1:, 2, 1:)
            gamma1 => gamma(1:, 1, 1:)
            gamma0 => gamma(1:, 0, 1:)
 
            gcp => globalcell(1:, 2, 1:)
            !---------------------------------------------------------------------------
 
        case (jmax)
 
            ww3 => w(1:, ny, 1:, :)
            ww2 => w(1:, jl, 1:, :)
            ww1 => w(1:, je, 1:, :)
            ww0 => w(1:, jb, 1:, :)
 
            pp3 => p(1:, ny, 1:)
            pp2 => p(1:, jl, 1:)
            pp1 => p(1:, je, 1:)
            pp0 => p(1:, jb, 1:)
 
            rlv3 => rlv(1:, ny, 1:)
            rlv2 => rlv(1:, jl, 1:)
            rlv1 => rlv(1:, je, 1:)
            rlv0 => rlv(1:, jb, 1:)
 
            rev3 => rev(1:, ny, 1:)
            rev2 => rev(1:, jl, 1:)
            rev1 => rev(1:, je, 1:)
            rev0 => rev(1:, jb, 1:)
 
            gamma3 => gamma(1:, ny, 1:)
            gamma2 => gamma(1:, jl, 1:)
            gamma1 => gamma(1:, je, 1:)
            gamma0 => gamma(1:, jb, 1:)
 
            gcp => globalcell(1:, jl, 1:)
            !---------------------------------------------------------------------------
 
        case (kmin)
 
            ww3 => w(1:, 1:, 3, :)
            ww2 => w(1:, 1:, 2, :)
            ww1 => w(1:, 1:, 1, :)
            ww0 => w(1:, 1:, 0, :)
 
            pp3 => p(1:, 1:, 3)
            pp2 => p(1:, 1:, 2)
            pp1 => p(1:, 1:, 1)
            pp0 => p(1:, 1:, 0)
 
            rlv3 => rlv(1:, 1:, 3)
            rlv2 => rlv(1:, 1:, 2)
            rlv1 => rlv(1:, 1:, 1)
            rlv0 => rlv(1:, 1:, 0)
 
            rev3 => rev(1:, 1:, 3)
            rev2 => rev(1:, 1:, 2)
            rev1 => rev(1:, 1:, 1)
            rev0 => rev(1:, 1:, 0)
 
            gamma3 => gamma(1:, 1:, 3)
            gamma2 => gamma(1:, 1:, 2)
            gamma1 => gamma(1:, 1:, 1)
            gamma0 => gamma(1:, 1:, 0)
 
            gcp => globalcell(1:, 1:, 2)
            !---------------------------------------------------------------------------
 
        case (kmax)
 
            ww3 => w(1:, 1:, nz, :)
            ww2 => w(1:, 1:, kl, :)
            ww1 => w(1:, 1:, ke, :)
            ww0 => w(1:, 1:, kb, :)
 
            pp3 => p(1:, 1:, nz)
            pp2 => p(1:, 1:, kl)
            pp1 => p(1:, 1:, ke)
            pp0 => p(1:, 1:, kb)
 
            rlv3 => rlv(1:, 1:, nz)
            rlv2 => rlv(1:, 1:, kl)
            rlv1 => rlv(1:, 1:, ke)
            rlv0 => rlv(1:, 1:, kb)
 
            rev3 => rev(1:, 1:, nz)
            rev2 => rev(1:, 1:, kl)
            rev1 => rev(1:, 1:, ke)
            rev0 => rev(1:, 1:, kb)
 
            gamma3 => gamma(1:, 1:, nz)
            gamma2 => gamma(1:, 1:, kl)
            gamma1 => gamma(1:, 1:, ke)
            gamma0 => gamma(1:, 1:, kb)
 
            gcp => globalcell(1:, 1:, kl)
        end select
 
        if (spatialpointers) then
            select case (bcfaceid(nn))
            case (imin)
                xx => x(1, :, :, :)
                ssi => si(1, :, :, :)
                ssj => sj(2, :, :, :)
                ssk => sk(2, :, :, :)
                ss => s(2, :, :, :)
            case (imax)
                xx => x(il, :, :, :)
                ssi => si(il, :, :, :)
                ssj => sj(il, :, :, :)
                ssk => sk(il, :, :, :)
                ss => s(il, :, :, :)
            case (jmin)
                xx => x(:, 1, :, :)
                ssi => sj(:, 1, :, :)
                ssj => si(:, 2, :, :)
                ssk => sk(:, 2, :, :)
                ss => s(:, 2, :, :)
            case (jmax)
                xx => x(:, jl, :, :)
                ssi => sj(:, jl, :, :)
                ssj => si(:, jl, :, :)
                ssk => sk(:, jl, :, :)
                ss => s(:, jl, :, :)
            case (kmin)
                xx => x(:, :, 1, :)
                ssi => sk(:, :, 1, :)
                ssj => si(:, :, 2, :)
                ssk => sj(:, :, 2, :)
                ss => s(:, :, 2, :)
            case (kmax)
                xx => x(:, :, kl, :)
                ssi => sk(:, :, kl, :)
                ssj => si(:, :, kl, :)
                ssk => sj(:, :, kl, :)
                ss => s(:, :, kl, :)
            end select
 
            if (addgridvelocities) then
                select case (bcfaceid(nn))
                case (imin)
                    sface => sfacei(1, :, :)
                case (imax)
                    sface => sfacei(il, :, :)
                case (jmin)
                    sface => sfacej(:, 1, :)
                case (jmax)
                    sface => sfacej(:, jl, :)
                case (kmin)
                    sface => sfacek(:, :, 1)
                case (kmax)
                    sface => sfacek(:, :, kl)
                end select
            end if
 
            if (equations == ransequations) then
                select case (bcfaceid(nn))
                case (imin)
                    dd2wall => d2wall(2, :, :)
                case (imax)
                    dd2wall => d2wall(il, :, :)
                case (jmin)
                    dd2wall => d2wall(:, 2, :)
                case (jmax)
                    dd2wall => d2wall(:, jl, :)
                case (kmin)
                    dd2wall => d2wall(:, :, 2)
                case (kmax)
                    dd2wall => d2wall(:, :, kl)
                end select
            end if
        end if
    end subroutine setbcpointers
 
    subroutine computerootbendingmoment(cf, cm, bendingMoment)
 
        !                                                      *
        ! Compute a normalized bending moment coefficient from *
        ! the force and moment coefficient. At the moment this *
        ! Routine only works for a half body. Additional logic *
        ! would be needed for a full body.                     *
        !                                                      *
 
        use constants
        use inputphysics, only: lengthref, pointref, pointrefec, liftindex
        implicit none
 
        !input/output variables
        real(kind=realtype), intent(in), dimension(3) :: cf, cm
        real(kind=realtype), intent(out) :: bendingmoment
 
        !Subroutine Variables
        real(kind=realtype) :: elasticmomentx, elasticmomenty, elasticmomentz
        bendingmoment = zero
        if (liftindex == 2) then
            !z out wing sum momentx,momentz
            elasticmomentx = cm(1) + cf(2) * (pointrefec(3) - &
                                              pointref(3)) / lengthref - cf(3) * &
                             (pointrefec(2) - pointref(2)) / lengthref
            elasticmomentz = cm(3) - cf(2) * (pointrefec(1) - &
                                              pointref(1)) / lengthref + cf(1) * &
                             (pointrefec(2) - pointref(2)) / lengthref
            bendingmoment = sqrt(elasticmomentx**2 + elasticmomentz**2)
        elseif (liftindex == 3) then
            !y out wing sum momentx,momenty
            elasticmomentx = cm(1) + cf(3) * (pointrefec(2) - &
                                              pointref(2)) / lengthref + &
                             cf(3) * (pointrefec(3) - pointref(3)) / lengthref
            elasticmomenty = cm(2) + cf(3) * (pointrefec(1) - &
                                              pointref(1)) / lengthref + &
                             cf(1) * (pointrefec(3) - pointref(3)) / lengthref
            bendingmoment = sqrt(elasticmomentx**2 + elasticmomenty**2)
        end if
 
    end subroutine computerootbendingmoment
 
    subroutine computeleastsquaresregression(y, x, npts, m, b)
        !
        !       Computes the slope of best fit for a set of x,y data of length
        !       npts
        !
        use constants
        implicit none
        !Subroutine arguments
        integer(kind=intType) :: npts
        real(kind=realtype), dimension(npts) :: x, y
        real(kind=realtype) :: m, b
 
        !local variables
        real(kind=realtype) :: sumx, sumy, sumx2, sumxy
        integer(kind=intType) :: i
 
        !begin execution
        sumx = 0.0
        sumy = 0.0
        sumx2 = 0.0
        sumxy = 0.0
        do i = 1, npts
 
            sumx = sumx + x(i)
            sumy = sumy + y(i)
            sumx2 = sumx2 + x(i) * x(i)
            sumxy = sumxy + x(i) * y(i)
        end do
 
        m = ((npts * sumxy) - (sumy * sumx)) / ((npts * sumx2) - (sumx)**2)
        b = (sumy * sumx2 - (sumx * sumxy)) / ((npts * sumx2) - (sumx)**2)
 
    end subroutine computeleastsquaresregression
 
    subroutine computetsderivatives(force, moment, coef0, dcdalpha, &
                                    dcdalphadot, dcdq, dcdqdot)
        !
        !      Computes the stability derivatives based on the time spectral
        !      solution of a given mesh. Takes in the force coefficients at
        !      all time instantces and computes the agregate parameters
        !
        use constants
        use communication
        use inputphysics
        use inputtimespectral
        use inputtsstabderiv
        use flowvarrefstate
        use monitor
        use section
        use inputmotion
        implicit none
 
        !
        !     Subroutine arguments.
        !
        real(kind=realtype), dimension(3, nTimeIntervalsSpectral) :: force, moment
        real(kind=realtype), dimension(8) :: dcdq, dcdqdot
        real(kind=realtype), dimension(8) :: dcdalpha, dcdalphadot
        real(kind=realtype), dimension(8) :: coef0
 
        ! Working Variables
        real(kind=realtype), dimension(nTimeIntervalsSpectral, 8) :: basecoef
        real(kind=realtype), dimension(8) :: coef0dot
        real(kind=realtype), dimension(nTimeIntervalsSpectral, 8) :: resbasecoef
        real(kind=realtype), dimension(nTimeIntervalsSpectral) :: intervalalpha, intervalalphadot
        real(kind=realtype), dimension(nTimeIntervalsSpectral) :: intervalmach, intervalmachdot
        real(kind=realtype), dimension(nSections) :: t
        integer(kind=intType) :: i, sps, nn
        !speed of sound: for normalization of q derivatives
        real(kind=realtype) :: a
        real(kind=realtype) :: fact, factmoment
        ! Functions
        real(kind=realtype), dimension(nTimeIntervalsSpectral) :: dphix, dphiy, dphiz
        real(kind=realtype), dimension(nTimeIntervalsSpectral) :: dphixdot, dphiydot, dphizdot
        real(kind=realtype) :: derivativerigidrotangle, secondderivativerigidrotangle
 
        fact = two / (gammainf * pinf * machcoef**2 &
                      * surfaceref * lref**2)
        factmoment = fact / (lengthref * lref)
 
        if (tsqmode) then
 
            print *, 'TS Q Mode code needs to be updated in computeTSDerivatives!'
            stop
 
            ! !q is pitch
            ! do sps =1,nTimeIntervalsSpectral
            !    !compute the time of this intervavc
            !    t = timeUnsteadyRestart
 
            !    if(equationMode == timeSpectral) then
            !       do nn=1,nSections
            !          t(nn) = t(nn) + (sps-1)*sections(nn)%timePeriod &
            !               /         (nTimeIntervalsSpectral*1.0)
            !       enddo
            !    endif
 
            !    ! Compute the time derivative of the rotation angles around the
            !    ! z-axis. i.e. compute q
 
            !    dphiZ(sps) = derivativeRigidRotAngle(degreePolZRot,   &
            !         coefPolZRot,     &
            !         degreeFourZRot,  &
            !         omegaFourZRot,   &
            !         cosCoefFourZRot, &
            !         sinCoefFourZRot, t)
 
            !    ! add in q_dot computation
            !    dphiZdot(sps) = secondDerivativeRigidRotAngle(degreePolZRot,   &
            !         coefPolZRot,     &
            !         degreeFourZRot,  &
            !         omegaFourZRot,   &
            !         cosCoefFourZRot, &
            !         sinCoefFourZRot, t)
            ! end do
 
            ! !now compute dCl/dq
            ! do i =1,8
            !    call computeLeastSquaresRegression(BaseCoef(:,i),dphiz,nTimeIntervalsSpectral,dcdq(i),coef0(i))
            ! end do
 
            ! ! now subtract off estimated cl,cmz and use remainder to compute
            ! ! clqdot and cmzqdot.
            ! do i = 1,8
            !    do sps = 1,nTimeIntervalsSpectral
            !       ResBaseCoef(sps,i) = BaseCoef(sps,i)-(dcdq(i)*dphiz(sps)+Coef0(i))
            !    enddo
            ! enddo
 
            ! !now normalize the results...
            ! a  = sqrt(gammaInf*pInfDim/rhoInfDim)
            ! dcdq = dcdq*timeRef*2*(machGrid*a)/lengthRef
 
            ! !now compute dCl/dpdot
            ! do i = 1,8
            !    call computeLeastSquaresRegression(ResBaseCoef(:,i),dphizdot,nTimeIntervalsSpectral,dcdqdot(i),Coef0dot(i))
            ! enddo
 
        elseif (tsalphamode) then
 
            do sps = 1, ntimeintervalsspectral
 
                !compute the time of this interval
                t = timeunsteadyrestart
 
                if (equationmode == timespectral) then
                    do nn = 1, nsections
                        t(nn) = t(nn) + (sps - 1) * sections(nn)%timePeriod &
                                / (ntimeintervalsspectral * 1.0)
                    end do
                end if
 
                intervalalpha(sps) = tsalpha(degreepolalpha, coefpolalpha, &
                                             degreefouralpha, omegafouralpha, &
                                             coscoeffouralpha, sincoeffouralpha, t(1))
 
                intervalalphadot(sps) = tsalphadot(degreepolalpha, coefpolalpha, &
                                                   degreefouralpha, omegafouralpha, &
                                                   coscoeffouralpha, sincoeffouralpha, t(1))
 
                ! THIS CALL IS WRONG!!!!
                !call getDirAngle(velDirFreestream,liftDirection,liftIndex,alpha+intervalAlpha(sps), beta)
 
                basecoef(sps, 1) = fact * ( &
                                   force(1, sps) * liftdirection(1) + &
                                   force(2, sps) * liftdirection(2) + &
                                   force(3, sps) * liftdirection(3))
                basecoef(sps, 2) = fact * ( &
                                   force(1, sps) * dragdirection(1) + &
                                   force(2, sps) * dragdirection(2) + &
                                   force(3, sps) * dragdirection(3))
                basecoef(sps, 3) = force(1, sps) * fact
                basecoef(sps, 4) = force(2, sps) * fact
                basecoef(sps, 5) = force(3, sps) * fact
                basecoef(sps, 6) = moment(1, sps) * factmoment
                basecoef(sps, 7) = moment(2, sps) * factmoment
                basecoef(sps, 8) = moment(3, sps) * factmoment
            end do
 
            !now compute dCl/dalpha
            do i = 1, 8
                call computeleastsquaresregression(basecoef(:, i), &
                                                   intervalalpha, ntimeintervalsspectral, dcdalpha(i), coef0(i))
            end do
 
            ! now subtract off estimated cl,cmz and use remainder to compute
            ! clalphadot and cmzalphadot.
            do i = 1, 8
                do sps = 1, ntimeintervalsspectral
                    resbasecoef(sps, i) = basecoef(sps, i) - (dcdalpha(i) * intervalalpha(sps) + coef0(i))
                end do
            end do
 
            !now compute dCi/dalphadot
            do i = 1, 8
                call computeleastsquaresregression(resbasecoef(:, i), &
                                                   intervalalphadot, ntimeintervalsspectral, &
                                                   dcdalphadot(i), coef0dot(i))
            end do
 
            a = sqrt(gammainf * pinfdim / rhoinfdim)
            dcdalphadot = dcdalphadot * 2 * (machgrid * a) / lengthref
 
        else
            call terminate('computeTSDerivatives', 'Not a valid stability motion')
        end if
 
    end subroutine computetsderivatives
 
    subroutine getdirangle(freeStreamAxis, liftAxis, liftIndex, alpha, beta)
        !
        !      Convert the wind axes to angle of attack and side slip angle.
        !      The direction angles alpha and beta are computed given the
        !      components of the wind direction vector (freeStreamAxis), the
        !      lift direction vector (liftAxis) and assuming that the
        !      body direction (xb,yb,zb) is in the default ijk coordinate
        !      system. The rotations are determined by first determining
        !      whether the lift is primarily in the j or k direction and then
        !      determining the angles accordingly.
        !      direction vector:
        !        1) Rotation about the zb or yb -axis: alpha clockwise (CW)
        !           (xb,yb,zb) -> (x1,y1,z1)
        !        2) Rotation about the yl or z1 -axis: beta counter-clockwise
        !           (CCW) (x1,y1,z1) -> (xw,yw,zw)
        !         input arguments:
        !            freeStreamAxis = wind vector in body axes
        !            liftAxis       = lift direction vector in body axis
        !         output arguments:
        !            alpha    = angle of attack in radians
        !            beta     = side slip angle in radians
        !
        use constants
 
        implicit none
        !
        !     Subroutine arguments.
        !
        !      real(kind=realType), intent(in)  :: xw, yw, zw
        real(kind=realtype), dimension(3), intent(in) :: freestreamaxis
        real(kind=realtype), dimension(3), intent(in) :: liftaxis
        real(kind=realtype), intent(out) :: alpha, beta
        integer(kind=intType), intent(out) :: liftIndex
        !
        !     Local variables.
        !
        real(kind=realtype) :: rnorm
        integer(kind=intType) :: flowIndex, i
        real(kind=realtype), dimension(3) :: freestreamaxisnorm
        integer(kind=intType) :: temp
 
        ! Assume domoniate flow is x
 
        flowindex = 1
 
        ! Determine the dominant lift direction
        if (abs(liftaxis(1)) > abs(liftaxis(2)) .and. &
            abs(liftaxis(1)) > abs(liftaxis(3))) then
            temp = 1
        else if (abs(liftaxis(2)) > abs(liftaxis(1)) .and. &
                 abs(liftaxis(2)) > abs(liftaxis(3))) then
            temp = 2
        else
            temp = 3
        end if
 
        liftindex = temp
 
        ! Normalize the freeStreamDirection vector.
        rnorm = sqrt(freestreamaxis(1)**2 + freestreamaxis(2)**2 + freestreamaxis(3)**2)
        do i = 1, 3
            freestreamaxisnorm(i) = freestreamaxis(i) / rnorm
        end do
 
        if (liftindex == 2) then
            ! different coordinate system for aerosurf
            ! Wing is in z- direction
            ! Compute angle of attack alpha.
 
            alpha = asin(freestreamaxisnorm(2))
 
            ! Compute side-slip angle beta.
 
            beta = -atan2(freestreamaxisnorm(3), freestreamaxisnorm(1))
 
        elseif (liftindex == 3) then
            ! Wing is in y- direction
 
            ! Compute angle of attack alpha.
 
            alpha = asin(freestreamaxisnorm(3))
 
            ! Compute side-slip angle beta.
 
            beta = atan2(freestreamaxisnorm(2), freestreamaxisnorm(1))
        else
            call terminate('getDirAngle', 'Invalid Lift Direction')
        end if
    end subroutine getdirangle
 
    subroutine stabilityderivativedriver
        !
        !      Runs the Time spectral stability derivative routines from the
        !      main program file
        !
        use precision
        implicit none
        !
        !     Local variables.
        !
        real(kind=realtype), dimension(8) :: dcdalpha, dcdalphadot, dcdbeta, &
                                             dcdbetadot, dcdmach, dcdmachdot
        real(kind=realtype), dimension(8) :: dcdp, dcdpdot, dcdq, dcdqdot, dcdr, dcdrdot
        real(kind=realtype), dimension(8) :: coef0, coef0dot
 
        !call computeTSDerivatives(coef0,dcdalpha,dcdalphadot,dcdq,dcdqdot)
 
    end subroutine stabilityderivativedriver
    subroutine setcoeftimeintegrator
        !
        !       setCoefTimeIntegrator determines the coefficients of the
        !       time integration scheme in unsteady mode. Normally these are
        !       equal to the coefficients corresponding to the specified
        !       accuracy. However during the initial phase there are not
        !       enough states in the past and the accuracy is reduced.
        !
        use constants
        use inputunsteady
        use inputphysics
        use iteration
        use monitor
        implicit none
        !
        !      Local variables.
        !
        integer(kind=intType) :: nn, nLevelsSet
 
        ! Determine which time integrator must be used.
 
        ! Modified by HDN
        select case (timeaccuracy)
        case (firstorder)
 
            ! 1st order. No need to check the number of available
            ! states in the past. Set the two coefficients and
            ! nLevelsSet to 2.
 
            coeftime(0) = 1.0_realtype
            coeftime(1) = -1.0_realtype
 
            if (useale .and. equationmode .eq. unsteady) then
                coeftimeale(1) = 1.0_realtype
                coefmeshale(1, 1) = half
                coefmeshale(1, 2) = half
            end if
 
            nlevelsset = 2
 
            !--------------------------------------------------
 
        case (secondorder)
 
            ! Second order time integrator. Determine the amount of
            ! available states and set the coefficients accordingly.
            select case (noldsolavail)
 
            case (1_inttype)
                coeftime(0) = 1.0_realtype
                coeftime(1) = -1.0_realtype
 
                if (useale .and. equationmode .eq. unsteady) then
                    coeftimeale(1) = half
                    coeftimeale(2) = half
                    coeftimeale(3) = zero
                    coeftimeale(4) = zero
 
                    coefmeshale(1, 1) = half
                    coefmeshale(1, 2) = half
                    coefmeshale(2, 1) = half
                    coefmeshale(2, 2) = half
                end if
 
                nlevelsset = 2
 
            case default   ! 2 or bigger.
                coeftime(0) = 1.5_realtype
                coeftime(1) = -2.0_realtype
                coeftime(2) = 0.5_realtype
 
                if (useale .and. equationmode .eq. unsteady) then
                    coeftimeale(1) = threefourth
                    coeftimeale(2) = threefourth
                    coeftimeale(3) = -fourth
                    coeftimeale(4) = -fourth
 
                    coefmeshale(1, 1) = half * (1.0_realtype + 1.0_realtype / sqrtthree)
                    coefmeshale(1, 2) = half * (1.0_realtype - 1.0_realtype / sqrtthree)
                    coefmeshale(2, 1) = coefmeshale(1, 2)
                    coefmeshale(2, 2) = coefmeshale(1, 1)
                end if
 
                nlevelsset = 3
 
            end select
 
            !--------------------------------------------------
 
        case (thirdorder)
 
            ! Third order time integrator.  Determine the amount of
            ! available states and set the coefficients accordingly.
 
            select case (noldsolavail)
 
            case (1_inttype)
                coeftime(0) = 1.0_realtype
                coeftime(1) = -1.0_realtype
 
                if (useale .and. equationmode .eq. unsteady) then
                    coeftimeale(1) = 1.0_realtype
                    coefmeshale(1, 1) = half
                    coefmeshale(1, 2) = half
                end if
 
                nlevelsset = 2
 
            case (2_inttype)
                coeftime(0) = 1.5_realtype
                coeftime(1) = -2.0_realtype
                coeftime(2) = 0.5_realtype
 
                if (useale .and. equationmode .eq. unsteady) then
                    coeftimeale(1) = threefourth
                    coeftimeale(2) = -fourth
                    coefmeshale(1, 1) = half * (1.0_realtype + 1.0_realtype / sqrtthree)
                    coefmeshale(1, 2) = half * (1.0_realtype - 1.0_realtype / sqrtthree)
                    coefmeshale(2, 1) = coefmeshale(1, 2)
                    coefmeshale(2, 2) = coefmeshale(1, 1)
                end if
 
                nlevelsset = 3
 
            case default   ! 3 or bigger.
                coeftime(0) = 11.0_realtype / 6.0_realtype
                coeftime(1) = -3.0_realtype
                coeftime(2) = 1.5_realtype
                coeftime(3) = -1.0_realtype / 3.0_realtype
 
                ! These numbers are NOT correct
                ! DO NOT use 3rd order ALE for now
                if (useale .and. equationmode .eq. unsteady) then
                    print *, 'Third-order ALE not implemented yet.'
                    coeftimeale(1) = threefourth
                    coeftimeale(2) = threefourth
                    coeftimeale(3) = -fourth
                    coeftimeale(4) = -fourth
                    coefmeshale(1, 1) = half * (1.0_realtype + 1.0_realtype / sqrtthree)
                    coefmeshale(1, 2) = half * (1.0_realtype - 1.0_realtype / sqrtthree)
                    coefmeshale(2, 1) = coefmeshale(1, 2)
                    coefmeshale(2, 2) = coefmeshale(1, 1)
                    coefmeshale(3, 1) = coefmeshale(1, 2)
                    coefmeshale(3, 2) = coefmeshale(1, 1)
                end if
 
                nlevelsset = 4
 
            end select
 
        end select
 
        ! Set the rest of the coefficients to 0 if not enough states
        ! in the past are available.
 
        do nn = nlevelsset, noldlevels
            coeftime(nn) = zero
        end do
 
    end subroutine setcoeftimeintegrator
 
    function mynorm2(x)
        use constants
        implicit none
        real(kind=realtype), dimension(3), intent(in) :: x
        real(kind=realtype) :: mynorm2
        mynorm2 = sqrt(x(1)**2 + x(2)**2 + x(3)**2)
    end function mynorm2
 
    function iswalltype(bType)
 
        use constants
        implicit none
        integer(kind=intType) :: btype
        logical :: iswalltype
 
        iswalltype = .false.
        if (btype == nswalladiabatic .or. &
            btype == nswallisothermal .or. &
            btype == eulerwall) then
            iswalltype = .true.
        end if
 
    end function iswalltype
 
    subroutine cross_prod(a, b, c)
 
        use precision
 
        ! Inputs
        real(kind=realtype), dimension(3), intent(in) :: a, b
 
        ! Outputs
        real(kind=realtype), dimension(3), intent(out) :: c
 
        c(1) = a(2) * b(3) - a(3) * b(2)
        c(2) = a(3) * b(1) - a(1) * b(3)
        c(3) = a(1) * b(2) - a(2) * b(1)
 
    end subroutine cross_prod
 
    subroutine siangle(angle, mult, trans)
 
        use constants
        use su_cgns, only: radian, degree
        implicit none
        !
        !      Subroutine arguments.
        !
        integer, intent(in) :: angle
        real(kind=realtype), intent(out) :: mult, trans
 
        ! Determine the situation we are having here.
 
        if (angle == radian) then
 
            ! Angle is already given in radIans. No need for a conversion.
 
            mult = one
            trans = zero
 
        else if (angle == degree) then
 
            ! Angle is given in degrees. A multiplication must be performed.
 
            mult = pi / 180.0_realtype
            trans = zero
 
        else
 
            call terminate("siAngle", &
                           "No idea how to convert this to SI units")
 
        end if
 
    end subroutine siangle
 
    subroutine sidensity(mass, len, mult, trans)
        !
        !       siDensity computes the conversion from the given density
        !       unit, which can be constructed from mass and length, to the
        !       SI-unit kg/m^3. The conversion will look like:
        !       density in kg/m^3 = mult*(density in NCU) + trans.
        !       NCU means non-christian units, i.e. everything that is not SI.
        !
        use constants
        use su_cgns, only: kilogram, meter
        implicit none
        !
        !      Subroutine arguments.
        !
        integer, intent(in) :: mass, len
        real(kind=realtype), intent(out) :: mult, trans
 
        ! Determine the situation we are having here.
 
        if (mass == kilogram .and. len == meter) then
 
            ! Density is given in kg/m^3, i.e. no need for a conversion.
 
            mult = one
            trans = zero
 
        else
 
            call terminate("siDensity", &
                           "No idea how to convert this to SI units")
 
        end if
 
    end subroutine sidensity
 
    subroutine silen(len, mult, trans)
        !
        !       siLen computes the conversion from the given length unit to
        !       the SI-unit meter. The conversion will look like:
        !       length in meter = mult*(length in NCU) + trans.
        !       NCU means non-christian units, i.e. everything that is not SI.
        !
        use constants
        use su_cgns, only: meter, centimeter, millimeter, foot, inch
        implicit none
        !
        !      Subroutine arguments.
        !
        integer, intent(in) :: len
        real(kind=realtype), intent(out) :: mult, trans
 
        ! Determine the situation we are having here.
 
        select case (len)
 
        case (meter)
            mult = one; trans = zero
 
        case (centimeter)
            mult = 0.01_realtype; trans = zero
 
        case (millimeter)
            mult = 0.001_realtype; trans = zero
 
        case (foot)
            mult = 0.3048_realtype; trans = zero
 
        case (inch)
            mult = 0.0254_realtype; trans = zero
 
        case default
            call terminate("siLen", &
                           "No idea how to convert this to SI units")
 
        end select
 
    end subroutine silen
 
    subroutine sipressure(mass, len, time, mult, trans)
        !
        !       siPressure computes the conversion from the given pressure
        !       unit, which can be constructed from mass, length and time, to
        !       the SI-unit Pa. The conversion will look like:
        !       pressure in Pa = mult*(pressure in NCU) + trans.
        !       NCU means non-christian units, i.e. everything that is not SI.
        !
        use constants
        use su_cgns, only: kilogram, meter, second
        implicit none
        !
        !      Subroutine arguments.
        !
        integer, intent(in) :: mass, len, time
        real(kind=realtype), intent(out) :: mult, trans
 
        ! Determine the situation we are having here.
 
        if (mass == kilogram .and. len == meter .and. time == second) then
 
            ! Pressure is given in Pa, i.e. no need for a conversion.
 
            mult = one
            trans = zero
 
        else
 
            call terminate("siPressure", &
                           "No idea how to convert this to SI units")
 
        end if
 
    end subroutine sipressure
 
    subroutine sitemperature(temp, mult, trans)
        !
        !       siTemperature computes the conversion from the given
        !       temperature unit to the SI-unit kelvin. The conversion will
        !       look like:
        !       temperature in K = mult*(temperature in NCU) + trans.
        !       NCU means non-christian units, i.e. everything that is not SI.
        !
        use constants
        use su_cgns, only: kelvin, celsius, rankine, fahrenheit
        implicit none
        !
        !      Subroutine arguments.
        !
        integer, intent(in) :: temp
        real(kind=realtype), intent(out) :: mult, trans
 
        ! Determine the situation we are having here.
 
        select case (temp)
 
        case (kelvin)
 
            ! Temperature is already given in Kelvin. No need to convert.
 
            mult = one
            trans = zero
 
        case (celsius)      ! is it Celcius or Celsius?
 
            ! Temperature is in Celsius. Only an offset must be applied.
 
            mult = one
            trans = 273.16_realtype
 
        case (rankine)
 
            ! Temperature is in Rankine. Only a multiplication needs to
            ! be performed.
 
            mult = 5.0_realtype / 9.0_realtype
            trans = zero
 
        case (fahrenheit)
 
            ! Temperature is in Fahrenheit. Both a multiplication and an
            ! offset must be applied.
 
            mult = 5.0_realtype / 9.0_realtype
            trans = 255.382
 
        case default
 
            ! Unknown temperature unit.
 
            call terminate("siTemperature", &
                           "No idea how to convert this to SI units")
 
        end select
 
    end subroutine sitemperature
    subroutine siturb(mass, len, time, temp, turbName, mult, trans)
        !
        !       siTurb computes the conversion from the given turbulence
        !       unit, which can be constructed from mass, len, time and temp,
        !       to the SI-unit for the given variable. The conversion will
        !       look like: var in SI = mult*(var in NCU) + trans.
        !       NCU means non-christian units, i.e. everything that is not SI.
        !
        use constants
        use su_cgns, only: kilogram, meter, second, kelvin
        implicit none
        !
        !      Subroutine arguments.
        !
        integer, intent(in) :: mass, len, time, temp
        character(len=*), intent(in) :: turbName
        real(kind=realtype), intent(out) :: mult, trans
 
        ! Determine the situation we are having here.
 
        if (mass == kilogram .and. len == meter .and. &
            time == second .and. temp == kelvin) then
 
            ! Everthing is already in SI units. No conversion needed.
 
            mult = one
            trans = zero
 
        else
 
            call terminate("siTurb", &
                           "No idea how to convert this to SI units")
 
        end if
 
    end subroutine siturb
 
    subroutine sivelocity(length, time, mult, trans)
        !
        !       siVelocity computes the conversion from the given velocity
        !       unit, which can be constructed from length and time, to the
        !       SI-unit m/s. The conversion will look like:
        !       velocity in m/s = mult*(velocity in ncu) + trans.
        !       Ncu means non-christian units, i.e. everything that is not SI.
        !
        use constants
        use su_cgns, only: meter, centimeter, millimeter, foot, inch, second
        implicit none
        !
        !      Subroutine arguments.
        !
        integer, intent(in) :: length, time
        real(kind=realtype), intent(out) :: mult, trans
 
        ! Determine the situation we are having here.
        ! First the length.
 
        select case (length)
 
        case (meter)
            mult = one; trans = zero
 
        case (centimeter)
            mult = 0.01_realtype; trans = zero
 
        case (millimeter)
            mult = 0.001_realtype; trans = zero
 
        case (foot)
            mult = 0.3048_realtype; trans = zero
 
        case (inch)
            mult = 0.0254_realtype; trans = zero
 
        case default
            call terminate("siVelocity", &
                           "No idea how to convert this length to SI units")
 
        end select
 
        ! And the time.
 
        select case (time)
 
        case (second)
            mult = mult
 
        case default
            call terminate("siVelocity", &
                           "No idea how to convert this time to SI units")
 
        end select
 
    end subroutine sivelocity
 
    ! ----------------------------------------------------------------------
    !                                                                      |
    !                    No Tapenade Routine below this line               |
    !                                                                      |
    ! ----------------------------------------------------------------------
 
#ifndef  USE_TAPENADE
    subroutine setbcpointers_d(nn, spatialpointers)
!
!       setbcpointers sets the pointers needed for the boundary
!       condition treatment on a general face, such that the boundary
!       routines are only implemented once instead of 6 times.
!
        use constants
        use blockpointers, only: w, wd, p, pd, rlv, rlvd, rev, revd, &
    &   gamma, x, xd, d2wall, d2walld, si, sid, sj, sjd, sk, skd, s, sd, &
    &   globalcell, bcdata, bcdatad, nx, il, ie, ib, ny, jl, je, jb, nz, kl,&
    &   ke, kb, bcfaceid, addgridvelocities, sfacei, sfaceid, sfacej, &
    &   sfacejd, sfacek, sfacekd, addgridvelocities
        use bcpointers_d, only: ww0, ww0d, ww1, ww1d, ww2, ww2d, ww3, ww3d,&
    &   pp0, pp0d, pp1, pp1d, pp2, pp2d, pp3, pp3d, rlv0, rlv0d, rlv1, rlv1d&
    &   , rlv2, rlv2d, rlv3, rlv3d, rev0, rev0d, rev1, rev1d, rev2, rev2d, &
    &   rev3, rev3d, gamma0, gamma1, gamma2, gamma3, gcp, xx, xxd, ss, ssd, &
    &   ssi, ssid, ssj, ssjd, ssk, sskd, dd2wall, sface, istart, iend, &
    &   jstart, jend, isize, jsize
        use inputphysics, only: cpmodel, equations
        implicit none
! subroutine arguments.
        integer(kind=inttype), intent(in) :: nn
        logical, intent(in) :: spatialpointers
! determine the sizes of each face and point to just the range we
! need on each face.
        istart = bcdata(nn)%icbeg
        iend = bcdata(nn)%icend
        jstart = bcdata(nn)%jcbeg
        jend = bcdata(nn)%jcend
! set the size of the subface
        isize = iend - istart + 1
        jsize = jend - jstart + 1
! determine the face id on which the subface is located and set
! the pointers accordinly.
        select case (bcfaceid(nn))
        case (imin)
!---------------------------------------------------------------------------
            ww3d => wd(3, 1:, 1:, :)
            ww3 => w(3, 1:, 1:, :)
            ww2d => wd(2, 1:, 1:, :)
            ww2 => w(2, 1:, 1:, :)
            ww1d => wd(1, 1:, 1:, :)
            ww1 => w(1, 1:, 1:, :)
            ww0d => wd(0, 1:, 1:, :)
            ww0 => w(0, 1:, 1:, :)
            pp3d => pd(3, 1:, 1:)
            pp3 => p(3, 1:, 1:)
            pp2d => pd(2, 1:, 1:)
            pp2 => p(2, 1:, 1:)
            pp1d => pd(1, 1:, 1:)
            pp1 => p(1, 1:, 1:)
            pp0d => pd(0, 1:, 1:)
            pp0 => p(0, 1:, 1:)
            rlv3d => rlvd(3, 1:, 1:)
            rlv3 => rlv(3, 1:, 1:)
            rlv2d => rlvd(2, 1:, 1:)
            rlv2 => rlv(2, 1:, 1:)
            rlv1d => rlvd(1, 1:, 1:)
            rlv1 => rlv(1, 1:, 1:)
            rlv0d => rlvd(0, 1:, 1:)
            rlv0 => rlv(0, 1:, 1:)
            rev3d => revd(3, 1:, 1:)
            rev3 => rev(3, 1:, 1:)
            rev2d => revd(2, 1:, 1:)
            rev2 => rev(2, 1:, 1:)
            rev1d => revd(1, 1:, 1:)
            rev1 => rev(1, 1:, 1:)
            rev0d => revd(0, 1:, 1:)
            rev0 => rev(0, 1:, 1:)
            gamma3 => gamma(3, 1:, 1:)
            gamma2 => gamma(2, 1:, 1:)
            gamma1 => gamma(1, 1:, 1:)
            gamma0 => gamma(0, 1:, 1:)
            gcp => globalcell(2, 1:, 1:)
        case (imax)
!---------------------------------------------------------------------------
            ww3d => wd(nx, 1:, 1:, :)
            ww3 => w(nx, 1:, 1:, :)
            ww2d => wd(il, 1:, 1:, :)
            ww2 => w(il, 1:, 1:, :)
            ww1d => wd(ie, 1:, 1:, :)
            ww1 => w(ie, 1:, 1:, :)
            ww0d => wd(ib, 1:, 1:, :)
            ww0 => w(ib, 1:, 1:, :)
            pp3d => pd(nx, 1:, 1:)
            pp3 => p(nx, 1:, 1:)
            pp2d => pd(il, 1:, 1:)
            pp2 => p(il, 1:, 1:)
            pp1d => pd(ie, 1:, 1:)
            pp1 => p(ie, 1:, 1:)
            pp0d => pd(ib, 1:, 1:)
            pp0 => p(ib, 1:, 1:)
            rlv3d => rlvd(nx, 1:, 1:)
            rlv3 => rlv(nx, 1:, 1:)
            rlv2d => rlvd(il, 1:, 1:)
            rlv2 => rlv(il, 1:, 1:)
            rlv1d => rlvd(ie, 1:, 1:)
            rlv1 => rlv(ie, 1:, 1:)
            rlv0d => rlvd(ib, 1:, 1:)
            rlv0 => rlv(ib, 1:, 1:)
            rev3d => revd(nx, 1:, 1:)
            rev3 => rev(nx, 1:, 1:)
            rev2d => revd(il, 1:, 1:)
            rev2 => rev(il, 1:, 1:)
            rev1d => revd(ie, 1:, 1:)
            rev1 => rev(ie, 1:, 1:)
            rev0d => revd(ib, 1:, 1:)
            rev0 => rev(ib, 1:, 1:)
            gamma3 => gamma(nx, 1:, 1:)
            gamma2 => gamma(il, 1:, 1:)
            gamma1 => gamma(ie, 1:, 1:)
            gamma0 => gamma(ib, 1:, 1:)
            gcp => globalcell(il, 1:, 1:)
        case (jmin)
!---------------------------------------------------------------------------
            ww3d => wd(1:, 3, 1:, :)
            ww3 => w(1:, 3, 1:, :)
            ww2d => wd(1:, 2, 1:, :)
            ww2 => w(1:, 2, 1:, :)
            ww1d => wd(1:, 1, 1:, :)
            ww1 => w(1:, 1, 1:, :)
            ww0d => wd(1:, 0, 1:, :)
            ww0 => w(1:, 0, 1:, :)
            pp3d => pd(1:, 3, 1:)
            pp3 => p(1:, 3, 1:)
            pp2d => pd(1:, 2, 1:)
            pp2 => p(1:, 2, 1:)
            pp1d => pd(1:, 1, 1:)
            pp1 => p(1:, 1, 1:)
            pp0d => pd(1:, 0, 1:)
            pp0 => p(1:, 0, 1:)
            rlv3d => rlvd(1:, 3, 1:)
            rlv3 => rlv(1:, 3, 1:)
            rlv2d => rlvd(1:, 2, 1:)
            rlv2 => rlv(1:, 2, 1:)
            rlv1d => rlvd(1:, 1, 1:)
            rlv1 => rlv(1:, 1, 1:)
            rlv0d => rlvd(1:, 0, 1:)
            rlv0 => rlv(1:, 0, 1:)
            rev3d => revd(1:, 3, 1:)
            rev3 => rev(1:, 3, 1:)
            rev2d => revd(1:, 2, 1:)
            rev2 => rev(1:, 2, 1:)
            rev1d => revd(1:, 1, 1:)
            rev1 => rev(1:, 1, 1:)
            rev0d => revd(1:, 0, 1:)
            rev0 => rev(1:, 0, 1:)
            gamma3 => gamma(1:, 3, 1:)
            gamma2 => gamma(1:, 2, 1:)
            gamma1 => gamma(1:, 1, 1:)
            gamma0 => gamma(1:, 0, 1:)
            gcp => globalcell(1:, 2, 1:)
        case (jmax)
!---------------------------------------------------------------------------
            ww3d => wd(1:, ny, 1:, :)
            ww3 => w(1:, ny, 1:, :)
            ww2d => wd(1:, jl, 1:, :)
            ww2 => w(1:, jl, 1:, :)
            ww1d => wd(1:, je, 1:, :)
            ww1 => w(1:, je, 1:, :)
            ww0d => wd(1:, jb, 1:, :)
            ww0 => w(1:, jb, 1:, :)
            pp3d => pd(1:, ny, 1:)
            pp3 => p(1:, ny, 1:)
            pp2d => pd(1:, jl, 1:)
            pp2 => p(1:, jl, 1:)
            pp1d => pd(1:, je, 1:)
            pp1 => p(1:, je, 1:)
            pp0d => pd(1:, jb, 1:)
            pp0 => p(1:, jb, 1:)
            rlv3d => rlvd(1:, ny, 1:)
            rlv3 => rlv(1:, ny, 1:)
            rlv2d => rlvd(1:, jl, 1:)
            rlv2 => rlv(1:, jl, 1:)
            rlv1d => rlvd(1:, je, 1:)
            rlv1 => rlv(1:, je, 1:)
            rlv0d => rlvd(1:, jb, 1:)
            rlv0 => rlv(1:, jb, 1:)
            rev3d => revd(1:, ny, 1:)
            rev3 => rev(1:, ny, 1:)
            rev2d => revd(1:, jl, 1:)
            rev2 => rev(1:, jl, 1:)
            rev1d => revd(1:, je, 1:)
            rev1 => rev(1:, je, 1:)
            rev0d => revd(1:, jb, 1:)
            rev0 => rev(1:, jb, 1:)
            gamma3 => gamma(1:, ny, 1:)
            gamma2 => gamma(1:, jl, 1:)
            gamma1 => gamma(1:, je, 1:)
            gamma0 => gamma(1:, jb, 1:)
            gcp => globalcell(1:, jl, 1:)
        case (kmin)
!---------------------------------------------------------------------------
            ww3d => wd(1:, 1:, 3, :)
            ww3 => w(1:, 1:, 3, :)
            ww2d => wd(1:, 1:, 2, :)
            ww2 => w(1:, 1:, 2, :)
            ww1d => wd(1:, 1:, 1, :)
            ww1 => w(1:, 1:, 1, :)
            ww0d => wd(1:, 1:, 0, :)
            ww0 => w(1:, 1:, 0, :)
            pp3d => pd(1:, 1:, 3)
            pp3 => p(1:, 1:, 3)
            pp2d => pd(1:, 1:, 2)
            pp2 => p(1:, 1:, 2)
            pp1d => pd(1:, 1:, 1)
            pp1 => p(1:, 1:, 1)
            pp0d => pd(1:, 1:, 0)
            pp0 => p(1:, 1:, 0)
            rlv3d => rlvd(1:, 1:, 3)
            rlv3 => rlv(1:, 1:, 3)
            rlv2d => rlvd(1:, 1:, 2)
            rlv2 => rlv(1:, 1:, 2)
            rlv1d => rlvd(1:, 1:, 1)
            rlv1 => rlv(1:, 1:, 1)
            rlv0d => rlvd(1:, 1:, 0)
            rlv0 => rlv(1:, 1:, 0)
            rev3d => revd(1:, 1:, 3)
            rev3 => rev(1:, 1:, 3)
            rev2d => revd(1:, 1:, 2)
            rev2 => rev(1:, 1:, 2)
            rev1d => revd(1:, 1:, 1)
            rev1 => rev(1:, 1:, 1)
            rev0d => revd(1:, 1:, 0)
            rev0 => rev(1:, 1:, 0)
            gamma3 => gamma(1:, 1:, 3)
            gamma2 => gamma(1:, 1:, 2)
            gamma1 => gamma(1:, 1:, 1)
            gamma0 => gamma(1:, 1:, 0)
            gcp => globalcell(1:, 1:, 2)
        case (kmax)
!---------------------------------------------------------------------------
            ww3d => wd(1:, 1:, nz, :)
            ww3 => w(1:, 1:, nz, :)
            ww2d => wd(1:, 1:, kl, :)
            ww2 => w(1:, 1:, kl, :)
            ww1d => wd(1:, 1:, ke, :)
            ww1 => w(1:, 1:, ke, :)
            ww0d => wd(1:, 1:, kb, :)
            ww0 => w(1:, 1:, kb, :)
            pp3d => pd(1:, 1:, nz)
            pp3 => p(1:, 1:, nz)
            pp2d => pd(1:, 1:, kl)
            pp2 => p(1:, 1:, kl)
            pp1d => pd(1:, 1:, ke)
            pp1 => p(1:, 1:, ke)
            pp0d => pd(1:, 1:, kb)
            pp0 => p(1:, 1:, kb)
            rlv3d => rlvd(1:, 1:, nz)
            rlv3 => rlv(1:, 1:, nz)
            rlv2d => rlvd(1:, 1:, kl)
            rlv2 => rlv(1:, 1:, kl)
            rlv1d => rlvd(1:, 1:, ke)
            rlv1 => rlv(1:, 1:, ke)
            rlv0d => rlvd(1:, 1:, kb)
            rlv0 => rlv(1:, 1:, kb)
            rev3d => revd(1:, 1:, nz)
            rev3 => rev(1:, 1:, nz)
            rev2d => revd(1:, 1:, kl)
            rev2 => rev(1:, 1:, kl)
            rev1d => revd(1:, 1:, ke)
            rev1 => rev(1:, 1:, ke)
            rev0d => revd(1:, 1:, kb)
            rev0 => rev(1:, 1:, kb)
            gamma3 => gamma(1:, 1:, nz)
            gamma2 => gamma(1:, 1:, kl)
            gamma1 => gamma(1:, 1:, ke)
            gamma0 => gamma(1:, 1:, kb)
            gcp => globalcell(1:, 1:, kl)
        end select
        if (spatialpointers) then
            select case (bcfaceid(nn))
            case (imin)
                xxd => xd(1, :, :, :)
                xx => x(1, :, :, :)
                ssid => sid(1, :, :, :)
                ssi => si(1, :, :, :)
                ssjd => sjd(2, :, :, :)
                ssj => sj(2, :, :, :)
                sskd => skd(2, :, :, :)
                ssk => sk(2, :, :, :)
                ssd => sd(2, :, :, :)
                ss => s(2, :, :, :)
            case (imax)
                xxd => xd(il, :, :, :)
                xx => x(il, :, :, :)
                ssid => sid(il, :, :, :)
                ssi => si(il, :, :, :)
                ssjd => sjd(il, :, :, :)
                ssj => sj(il, :, :, :)
                sskd => skd(il, :, :, :)
                ssk => sk(il, :, :, :)
                ssd => sd(il, :, :, :)
                ss => s(il, :, :, :)
            case (jmin)
                xxd => xd(:, 1, :, :)
                xx => x(:, 1, :, :)
                ssid => sjd(:, 1, :, :)
                ssi => sj(:, 1, :, :)
                ssjd => sid(:, 2, :, :)
                ssj => si(:, 2, :, :)
                sskd => skd(:, 2, :, :)
                ssk => sk(:, 2, :, :)
                ssd => sd(:, 2, :, :)
                ss => s(:, 2, :, :)
            case (jmax)
                xxd => xd(:, jl, :, :)
                xx => x(:, jl, :, :)
                ssid => sjd(:, jl, :, :)
                ssi => sj(:, jl, :, :)
                ssjd => sid(:, jl, :, :)
                ssj => si(:, jl, :, :)
                sskd => skd(:, jl, :, :)
                ssk => sk(:, jl, :, :)
                ssd => sd(:, jl, :, :)
                ss => s(:, jl, :, :)
            case (kmin)
                xxd => xd(:, :, 1, :)
                xx => x(:, :, 1, :)
                ssid => skd(:, :, 1, :)
                ssi => sk(:, :, 1, :)
                ssjd => sid(:, :, 2, :)
                ssj => si(:, :, 2, :)
                sskd => sjd(:, :, 2, :)
                ssk => sj(:, :, 2, :)
                ssd => sd(:, :, 2, :)
                ss => s(:, :, 2, :)
            case (kmax)
                xxd => xd(:, :, kl, :)
                xx => x(:, :, kl, :)
                ssid => skd(:, :, kl, :)
                ssi => sk(:, :, kl, :)
                ssjd => sid(:, :, kl, :)
                ssj => si(:, :, kl, :)
                sskd => sjd(:, :, kl, :)
                ssk => sj(:, :, kl, :)
                ssd => sd(:, :, kl, :)
                ss => s(:, :, kl, :)
            end select
            if (addgridvelocities) then
                select case (bcfaceid(nn))
                case (imin)
                    sface => sfacei(1, :, :)
                case (imax)
                    sface => sfacei(il, :, :)
                case (jmin)
                    sface => sfacej(:, 1, :)
                case (jmax)
                    sface => sfacej(:, jl, :)
                case (kmin)
                    sface => sfacek(:, :, 1)
                case (kmax)
                    sface => sfacek(:, :, kl)
                end select
            end if
            if (equations .eq. ransequations) then
                select case (bcfaceid(nn))
                case (imin)
                    dd2wall => d2wall(2, :, :)
                case (imax)
                    dd2wall => d2wall(il, :, :)
                case (jmin)
                    dd2wall => d2wall(:, 2, :)
                case (jmax)
                    dd2wall => d2wall(:, jl, :)
                case (kmin)
                    dd2wall => d2wall(:, :, 2)
                case (kmax)
                    dd2wall => d2wall(:, :, kl)
                end select
            end if
        end if
    end subroutine setbcpointers_d
 
    subroutine maxeddyv(eddyvisMax)
        !
        !       maxEddyv determines the maximum value of the eddy viscosity
        !       ratio of the block given by the pointers in blockPointes.
        !
        use constants
        use blockpointers, only: il, jl, kl, rlv, rev
        use flowvarrefstate, only: nwf, eddymodel
        implicit none
        !
        !      Subroutine arguments.
        !
        real(kind=realtype), intent(out) :: eddyvismax
        !
        !      Local variables.
        !
        integer(kind=intType) :: i, j, k
 
        real(kind=realtype) :: eddyvis
 
        ! Initialize the maximum value to zero and return immediately if
        ! not an eddy viscosity model is used.
 
        eddyvismax = zero
        if (.not. eddymodel) return
 
        ! Loop over the owned cells of this block.
 
        do k = 2, kl
            do j = 2, jl
                do i = 2, il
 
                    ! Compute the local viscosity ratio and take the maximum
                    ! with the currently stored value.
 
                    eddyvis = rev(i, j, k) / rlv(i, j, k)
                    eddyvismax = max(eddyvismax, eddyvis)
 
                end do
            end do
        end do
 
    end subroutine maxeddyv
 
    subroutine maxhdiffmach(hdiffMax, MachMax)
        !
        !       maxHdiffMach determines the maximum value of the Mach number
        !       and total enthalpy (or better the relative total enthalpy
        !       difference with the freestream).
        !
        use constants
        use blockpointers, only: il, jl, kl, w, p, gamma
        use flowvarrefstate, only: pinfcorr, rhoinf, winf
        use monitor, only: monmachorhmax
        implicit none
        !
        !      Subroutine arguments.
        !
        real(kind=realtype), intent(out) :: hdiffmax, machmax
        !
        !      Local variables.
        !
        integer(kind=intType) :: i, j, k
 
        real(kind=realtype) :: hdiff, hinf, mach2
 
        ! Initialize the maximum values to zero.
 
        hdiffmax = zero
        machmax = zero
 
        ! In case none of the two variables needs to be monitored,
        ! a return is made.
 
        if (.not. monmachorhmax) return
 
        ! Set the free stream value of the total enthalpy.
 
        hinf = (winf(irhoe) + pinfcorr) / rhoinf
 
        ! Loop over the owned cells of this block.
 
        do k = 2, kl
            do j = 2, jl
                do i = 2, il
 
                    ! Compute the local total enthalpy and Mach number squared.
 
                    hdiff = abs((w(i, j, k, irhoe) + p(i, j, k)) / w(i, j, k, irho) - hinf)
                    mach2 = (w(i, j, k, ivx)**2 + w(i, j, k, ivy)**2 &
                             + w(i, j, k, ivz)**2) * w(i, j, k, irho) / (gamma(i, j, k) * p(i, j, k))
 
                    ! Determine the maximum of these values and the
                    ! currently stored maximum values.
 
                    hdiffmax = max(hdiffmax, hdiff)
                    machmax = max(machmax, mach2)
 
                end do
            end do
        end do
 
        ! Currently the maximum Mach number squared is stored in
        ! MachMax. Take the square root. Also create a relative
        ! total enthalpy difference.
 
        machmax = sqrt(machmax)
        hdiffmax = hdiffmax / hinf
 
    end subroutine maxhdiffmach
 
    function delta(val1, val2)
        !
        !       delta is a function used to determine the contents of the full
        !       transformation matrix from the shorthand form. It returns 1
        !       if the absolute value of the two arguments are identical.
        !       Otherwise it returns 0.
        !
        use constants
        implicit none
        !
        !      Function type.
        !
        integer(kind=intType) :: delta
        !
        !      Function arguments.
        !
        integer(kind=intType) :: val1, val2
 
        if (abs(val1) == abs(val2)) then
            delta = 1_inttype
        else
            delta = 0_inttype
        end if
 
    end function delta
 
    subroutine nullifycgnsdompointers(nn)
        !
        !       nullifyCGNSDomPointers nullifies all the pointers of the
        !       given CGNS block.
        !
        use constants
        use cgnsgrid, only: cgnsdoms
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: nn
        !
        nullify (cgnsdoms(nn)%procStored)
        nullify (cgnsdoms(nn)%conn1to1)
        nullify (cgnsdoms(nn)%connNonMatchAbutting)
        nullify (cgnsdoms(nn)%bocoInfo)
 
    end subroutine nullifycgnsdompointers
 
    subroutine nullifyflowdompointers(nn, level, sps)
        !
        !       nullifyFlowDomPointers nullifies all the pointers of the
        !       given block.
        !
        use constants
        use block, only: flowdoms
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: nn, level, sps
 
        nullify (flowdoms(nn, level, sps)%BCType)
        nullify (flowdoms(nn, level, sps)%BCFaceID)
        nullify (flowdoms(nn, level, sps)%cgnsSubface)
 
        nullify (flowdoms(nn, level, sps)%inBeg)
        nullify (flowdoms(nn, level, sps)%jnBeg)
        nullify (flowdoms(nn, level, sps)%knBeg)
        nullify (flowdoms(nn, level, sps)%inEnd)
        nullify (flowdoms(nn, level, sps)%jnEnd)
        nullify (flowdoms(nn, level, sps)%knEnd)
 
        nullify (flowdoms(nn, level, sps)%dinBeg)
        nullify (flowdoms(nn, level, sps)%djnBeg)
        nullify (flowdoms(nn, level, sps)%dknBeg)
        nullify (flowdoms(nn, level, sps)%dinEnd)
        nullify (flowdoms(nn, level, sps)%djnEnd)
        nullify (flowdoms(nn, level, sps)%dknEnd)
 
        nullify (flowdoms(nn, level, sps)%icBeg)
        nullify (flowdoms(nn, level, sps)%jcBeg)
        nullify (flowdoms(nn, level, sps)%kcBeg)
        nullify (flowdoms(nn, level, sps)%icEnd)
        nullify (flowdoms(nn, level, sps)%jcEnd)
        nullify (flowdoms(nn, level, sps)%kcEnd)
 
        nullify (flowdoms(nn, level, sps)%neighBlock)
        nullify (flowdoms(nn, level, sps)%neighProc)
        nullify (flowdoms(nn, level, sps)%l1)
        nullify (flowdoms(nn, level, sps)%l2)
        nullify (flowdoms(nn, level, sps)%l3)
        nullify (flowdoms(nn, level, sps)%groupNum)
 
        nullify (flowdoms(nn, level, sps)%iblank)
        nullify (flowdoms(nn, level, sps)%forcedRecv)
        nullify (flowdoms(nn, level, sps)%status)
        nullify (flowdoms(nn, level, sps)%fringes)
        nullify (flowdoms(nn, level, sps)%orphans)
 
        nullify (flowdoms(nn, level, sps)%BCData)
        nullify (flowdoms(nn, level, sps)%viscSubface)
 
        nullify (flowdoms(nn, level, sps)%viscIminPointer)
        nullify (flowdoms(nn, level, sps)%viscImaxPointer)
        nullify (flowdoms(nn, level, sps)%viscJminPointer)
        nullify (flowdoms(nn, level, sps)%viscJmaxPointer)
        nullify (flowdoms(nn, level, sps)%viscKminPointer)
        nullify (flowdoms(nn, level, sps)%viscKmaxPointer)
 
        nullify (flowdoms(nn, level, sps)%x)
        nullify (flowdoms(nn, level, sps)%xOld)
        nullify (flowdoms(nn, level, sps)%si)
        nullify (flowdoms(nn, level, sps)%sj)
        nullify (flowdoms(nn, level, sps)%sk)
        nullify (flowdoms(nn, level, sps)%vol)
        nullify (flowdoms(nn, level, sps)%volRef)
        nullify (flowdoms(nn, level, sps)%volOld)
 
        nullify (flowdoms(nn, level, sps)%pori)
        nullify (flowdoms(nn, level, sps)%porj)
        nullify (flowdoms(nn, level, sps)%pork)
 
        nullify (flowdoms(nn, level, sps)%indFamilyI)
        nullify (flowdoms(nn, level, sps)%indFamilyJ)
        nullify (flowdoms(nn, level, sps)%indFamilyK)
 
        nullify (flowdoms(nn, level, sps)%factFamilyI)
        nullify (flowdoms(nn, level, sps)%factFamilyJ)
        nullify (flowdoms(nn, level, sps)%factFamilyK)
 
        nullify (flowdoms(nn, level, sps)%rotMatrixI)
        nullify (flowdoms(nn, level, sps)%rotMatrixJ)
        nullify (flowdoms(nn, level, sps)%rotMatrixK)
 
        nullify (flowdoms(nn, level, sps)%sFaceI)
        nullify (flowdoms(nn, level, sps)%sFaceJ)
        nullify (flowdoms(nn, level, sps)%sFaceK)
 
        nullify (flowdoms(nn, level, sps)%w)
        nullify (flowdoms(nn, level, sps)%wOld)
        nullify (flowdoms(nn, level, sps)%p)
        nullify (flowdoms(nn, level, sps)%aa)
        nullify (flowdoms(nn, level, sps)%gamma)
        nullify (flowdoms(nn, level, sps)%rlv)
        nullify (flowdoms(nn, level, sps)%rev)
        nullify (flowdoms(nn, level, sps)%s)
 
        nullify (flowdoms(nn, level, sps)%ux)
        nullify (flowdoms(nn, level, sps)%uy)
        nullify (flowdoms(nn, level, sps)%uz)
 
        nullify (flowdoms(nn, level, sps)%vx)
        nullify (flowdoms(nn, level, sps)%vy)
        nullify (flowdoms(nn, level, sps)%vz)
 
        nullify (flowdoms(nn, level, sps)%wx)
        nullify (flowdoms(nn, level, sps)%wy)
        nullify (flowdoms(nn, level, sps)%wz)
 
        nullify (flowdoms(nn, level, sps)%qx)
        nullify (flowdoms(nn, level, sps)%qy)
        nullify (flowdoms(nn, level, sps)%qz)
 
        nullify (flowdoms(nn, level, sps)%dw)
        nullify (flowdoms(nn, level, sps)%fw)
        nullify (flowdoms(nn, level, sps)%scratch)
        nullify (flowdoms(nn, level, sps)%shockSensor)
 
        nullify (flowdoms(nn, level, sps)%dwOldRK)
 
        nullify (flowdoms(nn, level, sps)%p1)
        nullify (flowdoms(nn, level, sps)%w1)
        nullify (flowdoms(nn, level, sps)%wr)
 
        nullify (flowdoms(nn, level, sps)%mgIFine)
        nullify (flowdoms(nn, level, sps)%mgJFine)
        nullify (flowdoms(nn, level, sps)%mgKFine)
 
        nullify (flowdoms(nn, level, sps)%mgIWeight)
        nullify (flowdoms(nn, level, sps)%mgJWeight)
        nullify (flowdoms(nn, level, sps)%mgKWeight)
 
        nullify (flowdoms(nn, level, sps)%mgICoarse)
        nullify (flowdoms(nn, level, sps)%mgJCoarse)
        nullify (flowdoms(nn, level, sps)%mgKCoarse)
 
        nullify (flowdoms(nn, level, sps)%ico)
        nullify (flowdoms(nn, level, sps)%jco)
        nullify (flowdoms(nn, level, sps)%kco)
 
        nullify (flowdoms(nn, level, sps)%wn)
        nullify (flowdoms(nn, level, sps)%pn)
        nullify (flowdoms(nn, level, sps)%dtl)
        nullify (flowdoms(nn, level, sps)%radI)
        nullify (flowdoms(nn, level, sps)%radJ)
        nullify (flowdoms(nn, level, sps)%radK)
 
        nullify (flowdoms(nn, level, sps)%d2Wall)
 
        nullify (flowdoms(nn, level, sps)%bmti1)
        nullify (flowdoms(nn, level, sps)%bmti2)
        nullify (flowdoms(nn, level, sps)%bmtj1)
        nullify (flowdoms(nn, level, sps)%bmtj2)
        nullify (flowdoms(nn, level, sps)%bmtk1)
        nullify (flowdoms(nn, level, sps)%bmtk2)
 
        nullify (flowdoms(nn, level, sps)%bvti1)
        nullify (flowdoms(nn, level, sps)%bvti2)
        nullify (flowdoms(nn, level, sps)%bvtj1)
        nullify (flowdoms(nn, level, sps)%bvtj2)
        nullify (flowdoms(nn, level, sps)%bvtk1)
        nullify (flowdoms(nn, level, sps)%bvtk2)
 
        nullify (flowdoms(nn, level, sps)%globalCell)
        nullify (flowdoms(nn, level, sps)%globalNode)
        nullify (flowdoms(nn, level, sps)%surfNodeIndices)
        nullify (flowdoms(nn, level, sps)%uv)
        nullify (flowdoms(nn, level, sps)%wallInd)
        nullify (flowdoms(nn, level, sps)%xSeed)
 
        ! Added by HDN
        nullify (flowdoms(nn, level, sps)%xALE)
        nullify (flowdoms(nn, level, sps)%sIALE)
        nullify (flowdoms(nn, level, sps)%sJALE)
        nullify (flowdoms(nn, level, sps)%sKALE)
        nullify (flowdoms(nn, level, sps)%sFaceIALE)
        nullify (flowdoms(nn, level, sps)%sFaceJALE)
        nullify (flowdoms(nn, level, sps)%sFaceKALE)
        nullify (flowdoms(nn, level, sps)%dwALE)
        nullify (flowdoms(nn, level, sps)%fwALE)
#ifndef USE_TAPENADE
        nullify (flowdoms(nn, level, sps)%PCMat)
        nullify (flowdoms(nn, level, sps)%i_D_Fact)
        nullify (flowdoms(nn, level, sps)%i_L_Fact)
        nullify (flowdoms(nn, level, sps)%i_U_Fact)
        nullify (flowdoms(nn, level, sps)%i_U2_Fact)
 
        nullify (flowdoms(nn, level, sps)%j_D_Fact)
        nullify (flowdoms(nn, level, sps)%j_L_Fact)
        nullify (flowdoms(nn, level, sps)%j_U_Fact)
        nullify (flowdoms(nn, level, sps)%j_U2_Fact)
 
        nullify (flowdoms(nn, level, sps)%k_D_Fact)
        nullify (flowdoms(nn, level, sps)%k_L_Fact)
        nullify (flowdoms(nn, level, sps)%k_U_Fact)
        nullify (flowdoms(nn, level, sps)%k_U2_Fact)
#endif
 
    end subroutine nullifyflowdompointers
 
    subroutine reallocateinteger(intArray, newSize, oldSize, &
                                 alwaysFreeMem)
        !
        !       reallocateInteger reallocates the given integer array to the
        !       given new size. The old values of the array are copied. Note
        !       that newSize can be both smaller and larger than oldSize.
        !
        use constants
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), dimension(:), pointer :: intArray
        integer(kind=intType), intent(in) :: newSize, oldSize
        logical, intent(in) :: alwaysFreeMem
        !
        !      Local variables.
        !
        integer(kind=intType), dimension(:), pointer :: tmp
 
        integer(kind=intType) :: i, nn, ll
 
        integer :: ierr
 
        ! Determine the minimum of newSize and oldSize.
 
        nn = min(newsize, oldsize)
 
        ! Set the pointer for tmp to intArray.
 
        tmp => intarray
 
        ! Allocate the memory for intArray in case newSize is larger
        ! than 0 or if alwaysFreeMem is .true. And copy the old data
        ! into it. Preserve the lower bound.
 
        if (newsize > 0 .or. alwaysfreemem) then
 
            ll = 1
            if (associated(intarray)) ll = lbound(intarray, 1)
 
            allocate (intarray(ll:newsize + ll - 1), stat=ierr)
            if (ierr /= 0) &
                call terminate("reallocateInteger", &
                               "Memory allocation failure for intArray")
            do i = ll, ll + nn - 1
                intarray(i) = tmp(i)
            end do
        end if
 
        ! Release the memory for tmp in case oldSize is larger than 0 or
        ! if alwaysFreeMem is .true.
 
        if (oldsize > 0 .or. alwaysfreemem) then
            deallocate (tmp, stat=ierr)
            if (ierr /= 0) &
                call terminate("reallocateInteger", &
                               "Deallocation error for tmp")
        end if
 
    end subroutine reallocateinteger
 
    !---------------------------------------------------------------------------=
 
    subroutine reallocatempioffsetkindinteger(intArray, newSize, &
                                              oldSize, alwaysFreeMem)
        !
        !       reallocateMpiOffsetKindInteger reallocates the given
        !       mpi_offset_kind integer array to the given new size. The old
        !       values of the array are copied. Note that newSize can be both
        !       smaller and larger than oldSize.
        !
        use constants
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=mpi_offset_kind), dimension(:), pointer :: intArray
        integer(kind=intType), intent(in) :: newSize, oldSize
        logical, intent(in) :: alwaysFreeMem
        !
        !      Local variables.
        !
        integer(kind=mpi_offset_kind), dimension(:), pointer :: tmp
 
        integer(kind=intType) :: i, nn, ll
 
        integer :: ierr
 
        ! Determine the minimum of newSize and oldSize.
 
        nn = min(newsize, oldsize)
 
        ! Set the pointer for tmp to intArray.
 
        tmp => intarray
 
        ! Allocate the memory for intArray in case newSize is larger
        ! than 0 or if alwaysFreeMem is .true. And copy the old data
        ! into it. Preserve the lower bound.
 
        if (newsize > 0 .or. alwaysfreemem) then
 
            ll = 1
            if (associated(intarray)) ll = lbound(intarray, 1)
 
            allocate (intarray(ll:newsize + ll - 1), stat=ierr)
            if (ierr /= 0) &
                call terminate("reallocateMpiOffsetKindInteger", &
                               "Memory allocation failure for intArray")
            do i = ll, ll + nn - 1
                intarray(i) = tmp(i)
            end do
        end if
 
        ! Release the memory for tmp in case oldSize is larger than 0 or
        ! if alwaysFreeMem is .true.
 
        if (oldsize > 0 .or. alwaysfreemem) then
            deallocate (tmp, stat=ierr)
            if (ierr /= 0) &
                call terminate("reallocateMpiOffsetKindInteger", &
                               "Deallocation error for tmp")
        end if
 
    end subroutine reallocatempioffsetkindinteger
 
    !---------------------------------------------------------------------------=
 
    subroutine reallocateinteger2(intArray, newSize1, newSize2, &
                                  oldSize1, oldSize2, &
                                  alwaysFreeMem)
        !
        !       reallocateInteger2 reallocates the given 2D integer array to
        !       the given new sizes. The old values of the array are copied.
        !       Note that the newSizes can be both smaller and larger than
        !       the oldSizes.
        !
        use constants
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), dimension(:, :), pointer :: intArray
        integer(kind=intType), intent(in) :: newSize1, newSize2, &
                                             oldSize1, oldSize2
        logical, intent(in) :: alwaysFreeMem
        !
        !      Local variables.
        !
        integer(kind=intType), dimension(:, :), pointer :: tmp
 
        integer(kind=intType) :: newSize, oldSize
        integer(kind=intType) :: nn1, nn2, nn
 
        integer(kind=intType) :: i, j
 
        integer :: ierr
 
        ! Determine the total new and old size.
 
        newsize = newsize1 * newsize2
        oldsize = oldsize1 * oldsize2
 
        ! Determine for each of the 2 components the minimum of the new
        ! and the old size. Multiply these values to obtain the total
        ! amount of data that must be copied.
 
        nn1 = min(newsize1, oldsize1)
        nn2 = min(newsize2, oldsize2)
 
        nn = nn1 * nn2
 
        ! Set the pointer for tmp.
 
        tmp => intarray
 
        ! Allocate the memory for intArray in case newSize is larger
        ! than 0 or if alwaysFreeMem is .true. and copy the old data
        ! into it.
 
        if (newsize > 0 .or. alwaysfreemem) then
            allocate (intarray(newsize1, newsize2), stat=ierr)
            if (ierr /= 0) &
                call terminate("reallocateInteger2", &
                               "Memory allocation failure for intArray")
            do j = 1, nn2
                do i = 1, nn1
                    intarray(i, j) = tmp(i, j)
                end do
            end do
        end if
 
        ! Release the memory of tmp in case oldSize is larger than 0
        ! or if alwaysFreeMem is .true..
 
        if (oldsize > 0 .or. alwaysfreemem) then
            deallocate (tmp, stat=ierr)
            if (ierr /= 0) &
                call terminate("reallocateInteger2", &
                               "Deallocation error for tmp")
        end if
 
    end subroutine reallocateinteger2
 
    subroutine reallocatereal(realArray, newSize, oldSize, &
                              alwaysFreeMem)
        !
        !       ReallocateReal reallocates the given real array to the given
        !       new size. The old values of the array are copied. Note that
        !       newSize can be both smaller and larger than oldSize.
        !
        use constants
        implicit none
        !
        !      Subroutine arguments.
        !
        real(kind=realtype), dimension(:), pointer :: realarray
        integer(kind=intType), intent(in) :: newSize, oldSize
        logical, intent(in) :: alwaysFreeMem
        !
        !      Local variables.
        !
        real(kind=realtype), dimension(:), pointer :: tmp
 
        integer(kind=intType) :: i, nn
 
        integer :: ierr
 
        ! Determine the minimum of newSize and oldSize.
 
        nn = min(newsize, oldsize)
 
        ! Set the pointer for tmp to realArray.
 
        tmp => realarray
 
        ! Allocate the memory for realArray in case newSize is larger
        ! than 0 or if alwaysFreeMem is .True. And copy the old data
        ! into it.
 
        if (newsize > 0 .or. alwaysfreemem) then
            allocate (realarray(newsize), stat=ierr)
            if (ierr /= 0) &
                call terminate("reallocateReal", &
                               "Memory allocation failure for realArray")
            do i = 1, nn
                realarray(i) = tmp(i)
            end do
        end if
 
        ! Release the memory for tmp in case oldSize is larger than 0 or
        ! if alwaysFreeMem is .True.
 
        if (oldsize > 0 .or. alwaysfreemem) then
            deallocate (tmp, stat=ierr)
            if (ierr /= 0) &
                call terminate("reallocateReal", &
                               "Deallocation error for tmp")
        end if
 
    end subroutine reallocatereal
 
    !---------------------------------------------------------------------------=
 
    subroutine reallocatereal2(realArray, newSize1, newSize2, &
                               oldSize1, oldSize2, &
                               alwaysFreeMem)
        !
        !       ReallocateReal2 reallocates the given 2d integer array to
        !       the given new sizes. The old values of the array are copied.
        !       Note that the newSizes can be both smaller and larger than
        !       the oldSizes.
        !
        use constants
        implicit none
        !
        !      Subroutine arguments.
        !
        real(kind=realtype), dimension(:, :), pointer :: realarray
        integer(kind=intType), intent(in) :: newSize1, newSize2, &
                                             oldSize1, oldSize2
        logical, intent(in) :: alwaysFreeMem
        !
        !      Local variables.
        !
        real(kind=realtype), dimension(:, :), pointer :: tmp
 
        integer(kind=intType) :: newSize, oldSize
        integer(kind=intType) :: nn1, nn2, nn
 
        integer(kind=intType) :: i, j
 
        integer :: ierr
 
        ! Determine the total new and old size.
 
        newsize = newsize1 * newsize2
        oldsize = oldsize1 * oldsize2
 
        ! Determine for each of the 2 components the minimum of the new
        ! and the old size. Multiply these values to obtain the total
        ! amount of data that must be copied.
 
        nn1 = min(newsize1, oldsize1)
        nn2 = min(newsize2, oldsize2)
 
        nn = nn1 * nn2
 
        ! Set the pointer for tmp.
 
        tmp => realarray
 
        ! Allocate the memory for realArray in case newSize is larger
        ! than 0 or if alwaysFreeMem is .True. And copy the old data
        ! into it.
 
        if (newsize > 0 .or. alwaysfreemem) then
            allocate (realarray(newsize1, newsize2), stat=ierr)
            if (ierr /= 0) &
                call terminate("reallocateReal2", &
                               "Memory allocation failure for realArray")
            do j = 1, nn2
                do i = 1, nn1
                    realarray(i, j) = tmp(i, j)
                end do
            end do
        end if
 
        ! Release the memory of tmp in case oldSize is larger than 0
        ! or if alwaysFreeMem is .True..
 
        if (oldsize > 0 .or. alwaysfreemem) then
            deallocate (tmp, stat=ierr)
            if (ierr /= 0) &
                call terminate("reallocateReal2", &
                               "Deallocation error for tmp")
        end if
 
    end subroutine reallocatereal2
 
    subroutine setbuffersizes(level, sps, determine1to1Buf, determineOversetBuf)
        !
        !       setBufferSizes determines the size of the send and receive
        !       buffers for this grid level. After that the maximum value of
        !       these sizes and the currently stored value is taken, such that
        !       for all mg levels the same buffer can be used. Normally the
        !       size on the finest grid should be enough, but it is just as
        !       safe to check on all mg levels.
        !
        use constants
        use communication, only: commpatternnode_1st, commpatterncell_2nd, &
                                 commpatternoverset, recvbuffersize, recvbuffersize_1to1, &
                                 recvbuffersizeover, recvbuffersizeover, sendbuffersize, recvbuffersize, &
                                 sendbuffersizeover, sendbuffersize_1to1
        use flowvarrefstate, only: nw, eddymodel, viscous
        use inputphysics, only: cpmodel
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: level, sps
        logical, intent(in) :: determine1to1Buf
        logical, intent(in) :: determineOversetBuf
        !
        !      Local variables.
        !
        integer(kind=intType) :: i
        integer(kind=intType) :: sendSize, recvSize, nVarComm
 
        ! Determine the maximum number of variables to be communicated.
 
        nvarcomm = nw + 1
        if (cpmodel == cptempcurvefits) nvarcomm = nvarcomm + 1
        if (viscous) nvarcomm = nvarcomm + 1
        if (eddymodel) nvarcomm = nvarcomm + 1
 
        ! Check if the 1 to 1 communication must be considered.
 
        if (determine1to1buf) then
 
            ! Store the send and receive buffer sizes needed for the nodal
            ! exchange. Determine the maximum for the number of send and
            ! receive processors.
 
            i = commpatternnode_1st(level)%nProcSend
            sendsize = commpatternnode_1st(level)%nsendCum(i)
 
            i = commpatternnode_1st(level)%nProcRecv
            recvsize = commpatternnode_1st(level)%nrecvCum(i)
 
            ! Determine the buffer sizes for the 2nd level cell exchange and
            ! set the size for this processor to the maximum needed. Note
            ! that it is not needed to test the 1st level cell halo, because
            ! it is entirely incorporated in the 2nd level.
            ! Determine the maximum for the number of send and receive
            ! processors as well.
 
            i = commpatterncell_2nd(level)%nProcSend
            sendsize = max(sendsize, &
                           commpatterncell_2nd(level)%nsendCum(i))
 
            i = commpatterncell_2nd(level)%nProcRecv
            recvsize = max(recvsize, &
                           commpatterncell_2nd(level)%nrecvCum(i))
 
            ! Multiply sendSize and recvSize with the number of variables to
            ! be communicated.
 
            sendsize = sendsize * nvarcomm
            recvsize = recvsize * nvarcomm
 
            ! Store the maximum of the current values and the old values
            ! in sendBufferSize1to1 and recvBufferSize1to1.
 
            sendbuffersize_1to1 = max(sendbuffersize_1to1, sendsize)
            recvbuffersize_1to1 = max(recvbuffersize_1to1, recvsize)
 
        end if
 
        ! Check if the overset communication must be considered.
 
        if (determineoversetbuf) then
 
            ! Same deal for the overset communication.
 
            i = commpatternoverset(level, sps)%nProcSend
            sendsize = commpatternoverset(level, sps)%nsendCum(i)
 
            i = commpatternoverset(level, sps)%nProcRecv
            recvsize = commpatternoverset(level, sps)%nrecvCum(i)
 
            ! Multiply sendSize and recvSize with the number of variables to
            ! be communicated.
 
            sendsize = sendsize * nvarcomm
            recvsize = recvsize * nvarcomm
 
            ! Store the maximum of the current values and the old values.
 
            sendbuffersizeover = max(sendbuffersizeover, sendsize)
            recvbuffersizeover = max(recvbuffersizeover, recvsize)
 
        end if
 
        ! Take the maximum for of all the buffers to
        ! obtain the actual size to be allocated.
 
        sendbuffersize = max(sendbuffersize_1to1, &
                             sendbuffersizeover)
        recvbuffersize = max(recvbuffersize_1to1, &
                             recvbuffersizeover)
 
    end subroutine setbuffersizes
 
    subroutine setpointers(nn, mm, ll)
        !
        !       setPointers makes the variables in blockPointers point to
        !       block nn for grid level mm and spectral solution ll.
        !
        ! Make an exception to use..only. We literally need everything
        ! from blockPointers so use a bare use.
        use constants
        use blockpointers
        implicit none
        !
        !      Subroutine arguments
        !
        integer(kind=intType), intent(in) :: nn, mm, ll
 
        ! Store the info of the current block, such that inside the
        ! module blockPointers it is known to which block the data
        ! belongs.
 
        sectionid = 1 ! We currently are only ever allowed 1 section
        nbklocal = nn
        nbkglobal = flowdoms(nn, mm, ll)%cgnsBlockID
        mglevel = mm
        spectralsol = ll
 
        ! Block dimensions.
 
        nx = flowdoms(nn, mm, ll)%nx
        ny = flowdoms(nn, mm, ll)%ny
        nz = flowdoms(nn, mm, ll)%nz
 
        il = flowdoms(nn, mm, ll)%il
        jl = flowdoms(nn, mm, ll)%jl
        kl = flowdoms(nn, mm, ll)%kl
 
        ie = flowdoms(nn, mm, ll)%ie
        je = flowdoms(nn, mm, ll)%je
        ke = flowdoms(nn, mm, ll)%ke
 
        ib = flowdoms(nn, mm, ll)%ib
        jb = flowdoms(nn, mm, ll)%jb
        kb = flowdoms(nn, mm, ll)%kb
 
        imaxdim = max(ie, je)
        jmaxdim = max(je, ke)
 
        righthanded = flowdoms(nn, mm, ll)%righthanded
 
        ! Point range in the corresponding cgns block
 
        ibegor = flowdoms(nn, mm, ll)%iBegor
        iendor = flowdoms(nn, mm, ll)%iEndor
        jbegor = flowdoms(nn, mm, ll)%jBegor
        jendor = flowdoms(nn, mm, ll)%jEndor
        kbegor = flowdoms(nn, mm, ll)%kBegor
        kendor = flowdoms(nn, mm, ll)%kEndor
 
        ! Subface info. Note that the pointers point to the 1st spectral
        ! mode, because this is the only one allocated. The info is the
        ! same for all modes.
 
        nsubface = flowdoms(nn, mm, ll)%nSubface
        n1to1 = flowdoms(nn, mm, ll)%n1to1
        nbocos = flowdoms(nn, mm, ll)%nBocos
        nviscbocos = flowdoms(nn, mm, ll)%nViscBocos
 
        bctype => flowdoms(nn, mm, 1)%BCType
        bcfaceid => flowdoms(nn, mm, 1)%BCFaceID
        cgnssubface => flowdoms(nn, mm, 1)%cgnsSubface
 
        inbeg => flowdoms(nn, mm, 1)%inBeg
        jnbeg => flowdoms(nn, mm, 1)%jnBeg
        knbeg => flowdoms(nn, mm, 1)%knBeg
        inend => flowdoms(nn, mm, 1)%inEnd
        jnend => flowdoms(nn, mm, 1)%jnEnd
        knend => flowdoms(nn, mm, 1)%knEnd
 
        dinbeg => flowdoms(nn, mm, 1)%dinBeg
        djnbeg => flowdoms(nn, mm, 1)%djnBeg
        dknbeg => flowdoms(nn, mm, 1)%dknBeg
        dinend => flowdoms(nn, mm, 1)%dinEnd
        djnend => flowdoms(nn, mm, 1)%djnEnd
        dknend => flowdoms(nn, mm, 1)%dknEnd
 
        icbeg => flowdoms(nn, mm, 1)%icBeg
        jcbeg => flowdoms(nn, mm, 1)%jcBeg
        kcbeg => flowdoms(nn, mm, 1)%kcBeg
        icend => flowdoms(nn, mm, 1)%icEnd
        jcend => flowdoms(nn, mm, 1)%jcEnd
        kcend => flowdoms(nn, mm, 1)%kcEnd
 
        neighblock => flowdoms(nn, mm, 1)%neighBlock
        neighproc => flowdoms(nn, mm, 1)%neighProc
        l1 => flowdoms(nn, mm, 1)%l1
        l2 => flowdoms(nn, mm, 1)%l2
        l3 => flowdoms(nn, mm, 1)%l3
        groupnum => flowdoms(nn, mm, 1)%groupNum
 
        ! Overset boundary and hole info.
        iblank => flowdoms(nn, mm, ll)%iblank
        status => flowdoms(nn, mm, ll)%status
        forcedrecv => flowdoms(nn, mm, ll)%forcedRecv
 
        fringes => flowdoms(nn, mm, ll)%fringes
        fringeptr => flowdoms(nn, mm, ll)%fringePtr
        gind => flowdoms(nn, mm, ll)%gInd
        ndonors => flowdoms(nn, mm, ll)%nDonors
 
        orphans => flowdoms(nn, mm, ll)%orphans
        norphans = flowdoms(nn, mm, ll)%nOrphans
 
        ! The data for boundary subfaces.
 
        bcdata => flowdoms(nn, mm, ll)%BCData
 
        ! The stress tensor and heat flux vector at viscous wall faces
        ! as well as the face pointers to these viscous wall faces.
        ! The latter point to the 1st spectral mode, because they are
        ! the only ones allocated. The info is the same for all modes.
 
        viscsubface => flowdoms(nn, mm, ll)%viscSubface
 
        visciminpointer => flowdoms(nn, mm, 1)%viscIminPointer
        viscimaxpointer => flowdoms(nn, mm, 1)%viscImaxPointer
        viscjminpointer => flowdoms(nn, mm, 1)%viscJminPointer
        viscjmaxpointer => flowdoms(nn, mm, 1)%viscJmaxPointer
        visckminpointer => flowdoms(nn, mm, 1)%viscKminPointer
        visckmaxpointer => flowdoms(nn, mm, 1)%viscKmaxPointer
 
        ! Mesh related variables. The porosities point to the 1st
        ! spectral mode, because they are the only ones allocated.
        ! The info is the same for all modes.
        ! Note that xOld and volOld always point to the finest
        ! grid level.
 
        x => flowdoms(nn, mm, ll)%x
        xold => flowdoms(nn, 1, ll)%xOld
 
        si => flowdoms(nn, mm, ll)%si
        sj => flowdoms(nn, mm, ll)%sj
        sk => flowdoms(nn, mm, ll)%sk
 
        vol => flowdoms(nn, mm, ll)%vol
        volref => flowdoms(nn, mm, ll)%volRef
        volold => flowdoms(nn, 1, ll)%volOld
 
        skew => flowdoms(nn, mm, ll)%skew
 
        pori => flowdoms(nn, mm, 1)%porI
        porj => flowdoms(nn, mm, 1)%porJ
        pork => flowdoms(nn, mm, 1)%porK
 
        indfamilyi => flowdoms(nn, mm, 1)%indFamilyI
        indfamilyj => flowdoms(nn, mm, 1)%indFamilyJ
        indfamilyk => flowdoms(nn, mm, 1)%indFamilyK
 
        factfamilyi => flowdoms(nn, mm, 1)%factFamilyI
        factfamilyj => flowdoms(nn, mm, 1)%factFamilyJ
        factfamilyk => flowdoms(nn, mm, 1)%factFamilyK
 
        rotmatrixi => flowdoms(nn, mm, ll)%rotMatrixI
        rotmatrixj => flowdoms(nn, mm, ll)%rotMatrixJ
        rotmatrixk => flowdoms(nn, mm, ll)%rotMatrixK
 
        blockismoving = flowdoms(nn, mm, ll)%blockIsMoving
        addgridvelocities = flowdoms(nn, mm, ll)%addGridVelocities
 
        sfacei => flowdoms(nn, mm, ll)%sFaceI
        sfacej => flowdoms(nn, mm, ll)%sFaceJ
        sfacek => flowdoms(nn, mm, ll)%sFaceK
 
        ! Flow variables. Note that wOld, gamma and the laminar viscosity
        ! point to the entries on the finest mesh. The reason is that
        ! they are computed from the other variables. For the eddy
        ! viscosity this is not the case because in a decoupled solver
        ! its values are obtained from the fine grid level.
 
        w => flowdoms(nn, mm, ll)%w
        wold => flowdoms(nn, 1, ll)%wOld
        p => flowdoms(nn, mm, ll)%p
        aa => flowdoms(nn, mm, ll)%aa
        shocksensor => flowdoms(nn, mm, ll)%shockSensor
 
        gamma => flowdoms(nn, 1, ll)%gamma
        rlv => flowdoms(nn, 1, ll)%rlv
        rev => flowdoms(nn, mm, ll)%rev
        s => flowdoms(nn, mm, ll)%s
 
        ux => flowdoms(nn, mm, ll)%ux
        uy => flowdoms(nn, mm, ll)%uy
        uz => flowdoms(nn, mm, ll)%uz
 
        vx => flowdoms(nn, mm, ll)%vx
        vy => flowdoms(nn, mm, ll)%vy
        vz => flowdoms(nn, mm, ll)%vz
 
        wx => flowdoms(nn, mm, ll)%wx
        wy => flowdoms(nn, mm, ll)%wy
        wz => flowdoms(nn, mm, ll)%wz
 
        qx => flowdoms(nn, mm, ll)%qx
        qy => flowdoms(nn, mm, ll)%qy
        qz => flowdoms(nn, mm, ll)%qz
 
        ! Residual and multigrid variables. The residual point to the
        ! finest grid entry, the multigrid variables to their own level.
 
        dw => flowdoms(nn, 1, ll)%dw
        fw => flowdoms(nn, 1, ll)%fw
        dwoldrk => flowdoms(nn, 1, ll)%dwOldRK
        scratch => flowdoms(nn, 1, ll)%scratch
 
        p1 => flowdoms(nn, mm, ll)%p1
        w1 => flowdoms(nn, mm, ll)%w1
        wr => flowdoms(nn, mm, ll)%wr
 
        ! Variables, which allow a more flexible multigrid treatment.
        ! They are the same for all spectral modes and therefore they
        ! point to the 1st mode.
 
        mgifine => flowdoms(nn, mm, 1)%mgIFine
        mgjfine => flowdoms(nn, mm, 1)%mgJFine
        mgkfine => flowdoms(nn, mm, 1)%mgKFine
 
        mgiweight => flowdoms(nn, mm, 1)%mgIWeight
        mgjweight => flowdoms(nn, mm, 1)%mgJWeight
        mgkweight => flowdoms(nn, mm, 1)%mgKWeight
 
        mgicoarse => flowdoms(nn, mm, 1)%mgICoarse
        mgjcoarse => flowdoms(nn, mm, 1)%mgJCoarse
        mgkcoarse => flowdoms(nn, mm, 1)%mgKCoarse
 
        ! Time-stepping variables and spectral radIi.
        ! They all point to the fine mesh entry.
 
        wn => flowdoms(nn, 1, ll)%wn
        pn => flowdoms(nn, 1, ll)%pn
        dtl => flowdoms(nn, 1, ll)%dtl
 
        radi => flowdoms(nn, 1, ll)%radI
        radj => flowdoms(nn, 1, ll)%radJ
        radk => flowdoms(nn, 1, ll)%radK
 
        ! Wall distance for the turbulence models.
 
        d2wall => flowdoms(nn, mm, ll)%d2Wall
        filterdes => flowdoms(nn, mm, ll)%filterDES  ! eran-des
 
        ! Arrays used for the implicit treatment of the turbulent wall
        ! boundary conditions. As these variables are only allocated for
        ! the 1st spectral solution of the fine mesh, the pointers point
        ! to those arrays.
 
        bmti1 => flowdoms(nn, 1, 1)%bmti1
        bmti2 => flowdoms(nn, 1, 1)%bmti2
        bmtj1 => flowdoms(nn, 1, 1)%bmtj1
        bmtj2 => flowdoms(nn, 1, 1)%bmtj2
        bmtk1 => flowdoms(nn, 1, 1)%bmtk1
        bmtk2 => flowdoms(nn, 1, 1)%bmtk2
 
        bvti1 => flowdoms(nn, 1, 1)%bvti1
        bvti2 => flowdoms(nn, 1, 1)%bvti2
        bvtj1 => flowdoms(nn, 1, 1)%bvtj1
        bvtj2 => flowdoms(nn, 1, 1)%bvtj2
        bvtk1 => flowdoms(nn, 1, 1)%bvtk1
        bvtk2 => flowdoms(nn, 1, 1)%bvtk2
 
        ! Pointers for globalCell/Node
        globalcell => flowdoms(nn, mm, ll)%globalCell
        globalnode => flowdoms(nn, mm, ll)%globalNode
 
        xseed => flowdoms(nn, mm, ll)%xSeed
        wallind => flowdoms(nn, mm, ll)%wallInd
 
        ! Added by HDN
        ! Kept the same dim as their counterparts
        xale => flowdoms(nn, mm, ll)%xALE
        sveloiale => flowdoms(nn, mm, ll)%sVeloIALE
        svelojale => flowdoms(nn, mm, ll)%sVeloJALE
        svelokale => flowdoms(nn, mm, ll)%sVeloKALE
        siale => flowdoms(nn, mm, ll)%sIALE
        sjale => flowdoms(nn, mm, ll)%sJALE
        skale => flowdoms(nn, mm, ll)%sKALE
        sfaceiale => flowdoms(nn, mm, ll)%sFaceIALE
        sfacejale => flowdoms(nn, mm, ll)%sFaceJALE
        sfacekale => flowdoms(nn, mm, ll)%sFaceKALE
        dwale => flowdoms(nn, 1, ll)%dwALE
        fwale => flowdoms(nn, 1, ll)%fwALE
 
        ! Pointers for PC
        pcmat => flowdoms(nn, mm, ll)%pcMat
 
        i_d_fact => flowdoms(nn, mm, ll)%i_D_fact
        i_l_fact => flowdoms(nn, mm, ll)%i_L_fact
        i_u_fact => flowdoms(nn, mm, ll)%i_U_fact
        i_u2_fact => flowdoms(nn, mm, ll)%i_U2_fact
 
        j_d_fact => flowdoms(nn, mm, ll)%j_D_fact
        j_l_fact => flowdoms(nn, mm, ll)%j_L_fact
        j_u_fact => flowdoms(nn, mm, ll)%j_U_fact
        j_u2_fact => flowdoms(nn, mm, ll)%j_U2_fact
 
        k_d_fact => flowdoms(nn, mm, ll)%k_D_fact
        k_l_fact => flowdoms(nn, mm, ll)%k_L_fact
        k_u_fact => flowdoms(nn, mm, ll)%k_U_fact
        k_u2_fact => flowdoms(nn, mm, ll)%k_U2_fact
 
        pcvec1 => flowdoms(nn, mm, ll)%PCVec1
        pcvec2 => flowdoms(nn, mm, ll)%PCVec2
 
        i_ipiv => flowdoms(nn, mm, ll)%i_ipiv
        j_ipiv => flowdoms(nn, mm, ll)%j_ipiv
        k_ipiv => flowdoms(nn, mm, ll)%k_ipiv
 
    end subroutine setpointers
 
    subroutine setpointers_b(nn, level, sps)
        use constants
        implicit none
        integer(kind=intType), intent(in) :: nn, level, sps
 
        ! Alias for setPonters_d
        call setpointers_d(nn, level, sps)
 
    end subroutine setpointers_b
 
    ! Set the pointers for the derivative values AND the normal pointers
    subroutine setpointers_d(nn, level, sps)
 
        use block, only: flowdomsd
        use blockpointers
        implicit none
        !
        !      Subroutine arguments
        !
        integer(kind=intType), intent(in) :: nn, level, sps
 
        ! Set normal pointers
        call setpointers(nn, level, sps)
 
        viscsubfaced => flowdomsd(nn, 1, sps)%viscSubface
 
        xd => flowdomsd(nn, 1, sps)%x
 
        sid => flowdomsd(nn, 1, sps)%si
        sjd => flowdomsd(nn, 1, sps)%sj
        skd => flowdomsd(nn, 1, sps)%sk
 
        vold => flowdomsd(nn, 1, sps)%vol
 
        rotmatrixid => flowdomsd(nn, 1, sps)%rotMatrixI
        rotmatrixjd => flowdomsd(nn, 1, sps)%rotMatrixJ
        rotmatrixkd => flowdomsd(nn, 1, sps)%rotMatrixK
 
        sfaceid => flowdomsd(nn, 1, sps)%sFaceI
        sfacejd => flowdomsd(nn, 1, sps)%sFaceJ
        sfacekd => flowdomsd(nn, 1, sps)%sFaceK
 
        ! Flow variables. Note that wOld, gamma and the laminar viscosity
        ! point to the entries on the finest mesh. The reason is that
        ! they are computed from the other variables. For the eddy
        ! viscosity this is not the case because in a decoupled solver
        ! its values are obtained from the fine grid level.
 
        wd => flowdomsd(nn, 1, sps)%w
        pd => flowdomsd(nn, 1, sps)%p
 
        gammad => flowdomsd(nn, 1, sps)%gamma
        aad => flowdomsd(nn, 1, sps)%aa
        rlvd => flowdomsd(nn, 1, sps)%rlv
        revd => flowdomsd(nn, 1, sps)%rev
        sd => flowdomsd(nn, 1, sps)%s
 
        uxd => flowdomsd(nn, 1, sps)%ux
        uyd => flowdomsd(nn, 1, sps)%uy
        uzd => flowdomsd(nn, 1, sps)%uz
 
        vxd => flowdomsd(nn, 1, sps)%vx
        vyd => flowdomsd(nn, 1, sps)%vy
        vzd => flowdomsd(nn, 1, sps)%vz
 
        wxd => flowdomsd(nn, 1, sps)%wx
        wyd => flowdomsd(nn, 1, sps)%wy
        wzd => flowdomsd(nn, 1, sps)%wz
 
        qxd => flowdomsd(nn, 1, sps)%qx
        qyd => flowdomsd(nn, 1, sps)%qy
        qzd => flowdomsd(nn, 1, sps)%qz
 
        ! Residual and multigrid variables. The residual point to the
        ! finest grid entry, the multigrid variables to their own level.
 
        dwd => flowdomsd(nn, 1, sps)%dw
        fwd => flowdomsd(nn, 1, sps)%fw
        scratchd => flowdomsd(nn, 1, sps)%scratch
 
        dtld => flowdomsd(nn, 1, sps)%dtl
 
        ! Time-stepping variables and spectral radIi.
        ! They asps point to the fine mesh entry.
 
        radid => flowdomsd(nn, 1, sps)%radI
        radjd => flowdomsd(nn, 1, sps)%radJ
        radkd => flowdomsd(nn, 1, sps)%radK
 
        d2walld => flowdomsd(nn, 1, sps)%d2Wall
 
        ! Arrays used for the implicit treatment of the turbulent wasps
        ! boundary conditions. As these variables are only aspocated for
        ! the 1st spectral solution of the fine mesh, the pointers point
        ! to those arrays.
 
        bmti1d => flowdomsd(nn, 1, 1)%bmti1
        bmti2d => flowdomsd(nn, 1, 1)%bmti2
        bmtj1d => flowdomsd(nn, 1, 1)%bmtj1
        bmtj2d => flowdomsd(nn, 1, 1)%bmtj2
        bmtk1d => flowdomsd(nn, 1, 1)%bmtk1
        bmtk2d => flowdomsd(nn, 1, 1)%bmtk2
 
        bvti1d => flowdomsd(nn, 1, 1)%bvti1
        bvti2d => flowdomsd(nn, 1, 1)%bvti2
        bvtj1d => flowdomsd(nn, 1, 1)%bvtj1
        bvtj2d => flowdomsd(nn, 1, 1)%bvtj2
        bvtk1d => flowdomsd(nn, 1, 1)%bvtk1
        bvtk2d => flowdomsd(nn, 1, 1)%bvtk2
 
        !BCData Array
        bcdatad => flowdomsd(nn, 1, sps)%BCdata
 
    end subroutine setpointers_d
 
    subroutine spectralinterpolcoef(nsps, t, alpScal, alpMat)
        !
        !       spectralInterpolCoef determines the scalar and matrix
        !       spectral interpolation coefficients for the given number of
        !       spectral solutions for the given t, where t is the ratio of
        !       the time and the periodic interval time. Note that the index
        !       of the spectral solutions of both alpScal and alpMat start
        !       at 0. In this way these coefficients are easier to determine.
        !
        use constants
        use inputtimespectral, only: ntimeintervalsspectral, rotmatrixspectral
        use section, only: nsections, sections
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: nsps
        real(kind=realtype), intent(in) :: t
 
        real(kind=realtype), dimension(0:nsps - 1), intent(out) :: alpscal
        real(kind=realtype), dimension(nSections, 0:nsps - 1, 3, 3), &
            intent(out) :: alpmat
        !
        !      Local variables.
        !
        integer(kind=intType) :: jj, nn, j, p, r, nhalfM1, m, mhalfM1
 
        real(kind=realtype) :: nspsinv, minv, tm, alp
 
        real(kind=realtype), dimension(3, 3) :: rp, tmp
 
        !       Scalar coefficients.
        !
        ! Loop over the number of spectral solutions to compute the
        ! coefficients. Note that the loop starts at 0.
 
        if (mod(nsps, 2) .eq. 0) then
            nhalfm1 = nsps / 2 - 1
        else
            nhalfm1 = (nsps - 1) / 2
        end if
 
        nspsinv = one / real(nsps, realtype)
 
        do j = 0, (nsps - 1)
            if (mod(nsps, 2) .eq. 0) then
                alpscal(j) = one + cos(j * pi) * cos(nsps * pi * t)
            else
                alpscal(j) = one + cos(j * pi * (nsps + 1) / nsps) * cos((nsps + 1) * pi * t)
            end if
 
            do r = 1, nhalfm1
                alpscal(j) = alpscal(j) &
                             + two * cos(r * j * two * pi * nspsinv) * cos(r * two * pi * t) &
                             + two * sin(r * j * two * pi * nspsinv) * sin(r * two * pi * t)
            end do
 
            alpscal(j) = alpscal(j) * nspsinv
 
        end do
        !
        !       Matrix coefficients. These are (can be) different for every
        !       section and they must therefore be determined for every
        !       section.
        !
        ! Loop over the number of sections in the grid.
 
        sectionloop: do nn = 1, nsections
 
            ! Compute the numbers for the entire wheel for this section.
            ! Note that also t must be adapted, because t is a ratio between
            ! the actual time and the periodic time.
 
            m = nsps * sections(nn)%nSlices
            if (mod(m, 2) .eq. 0) then
                mhalfm1 = m / 2 - 1
            else
                mhalfm1 = (m - 1) / 2
            end if
            minv = one / real(m, realtype)
            tm = t / real(sections(nn)%nSlices, realtype)
 
            ! Loop over the number of spectral solutions.
 
            spectralloop: do jj = 0, (nsps - 1)
 
                ! Initialize the matrix coefficients to zero and the matrix
                ! rp to the identity matrix. Rp is the rotation matrix of this
                ! section to the power p, which starts at 0, i.e. rp = i.
 
                alpmat(nn, jj, 1, 1) = zero
                alpmat(nn, jj, 1, 2) = zero
                alpmat(nn, jj, 1, 3) = zero
 
                alpmat(nn, jj, 2, 1) = zero
                alpmat(nn, jj, 2, 2) = zero
                alpmat(nn, jj, 2, 3) = zero
 
                alpmat(nn, jj, 3, 1) = zero
                alpmat(nn, jj, 3, 2) = zero
                alpmat(nn, jj, 3, 3) = zero
 
                rp(1, 1) = one
                rp(1, 2) = zero
                rp(1, 3) = zero
 
                rp(2, 1) = zero
                rp(2, 2) = one
                rp(2, 3) = zero
 
                rp(3, 1) = zero
                rp(3, 2) = zero
                rp(3, 3) = one
 
                ! Loop over the number of slices of this section. Note that
                ! this loop starts at zero, which simplifies the formulas.
 
                slicesloop: do p = 0, (sections(nn)%nSlices - 1)
 
                    ! Determine the index j, the index of alp in the entire
                    ! wheel.
 
                    j = jj + p * nsps
 
                    ! Compute the scalar coefficient alp of the index j in
                    ! the entire wheel.
 
                    if (mod(m, 2) .eq. 0) then
                        alp = one + cos(j * pi) * cos(m * pi * tm)
                    else
                        alp = one + cos(j * pi * (m + 1) / m) * cos((m + 1) * pi * tm)
                    end if
                    do r = 1, mhalfm1
                        alp = alp + two * cos(r * j * two * pi * minv) * cos(r * two * pi * tm) &
                              + two * sin(r * j * two * pi * minv) * sin(r * two * pi * tm)
                    end do
 
                    alp = alp * minv
 
                    ! Update the matrix coefficient.
 
                    do r = 1, 3
                        do j = 1, 3
                            alpmat(nn, jj, r, j) = alpmat(nn, jj, r, j) + alp * rp(r, j)
                        end do
                    end do
 
                    ! Multiply rp by the rotation matrix to obtain the correct
                    ! matrix for the next slice. Use tmp as temporary storage.
 
                    do r = 1, 3
                        do j = 1, 3
                            tmp(r, j) = rp(r, 1) * rotmatrixspectral(nn, 1, j) &
                                        + rp(r, 2) * rotmatrixspectral(nn, 2, j) &
                                        + rp(r, 3) * rotmatrixspectral(nn, 3, j)
                        end do
                    end do
 
                    rp = tmp
 
                end do slicesloop
            end do spectralloop
        end do sectionloop
 
    end subroutine spectralinterpolcoef
 
    subroutine deallocatetempmemory(resNeeded)
        !
        !       deallocateTempMemory deallocates memory used in the solver,
        !       but which is not needed to store the actual solution. In this
        !       way the memory can be used differently, e.g. when writing the
        !       solution or computing the wall distances.
        !
        use constants
        use block, only: flowdoms, ndom
        use communication, only: sendbuffer, recvbuffer
        use inputiteration, only: smoother
        use inputtimespectral, only: ntimeintervalsspectral
        implicit none
        !
        !      Subroutine arguments.
        !
        logical, intent(in) :: resNeeded
        !
        !      Local variables.
        !
        integer :: ierr
 
        integer(kind=intType) :: nn, mm
 
        ! Deallocate the communication buffers
 
        deallocate (sendbuffer, recvbuffer, stat=ierr)
        if (ierr /= 0) &
            call terminate("deallocateTempMemory", &
                           "Deallocation error for communication buffers")
 
        ! Loop over the spectral modes and domains. Note that only memory
        ! on the finest grid is released, because a) most of these
        ! variables are only allocated on the fine grid and b) the coarser
        ! grids do not contribute that much in the memory usage anyway.
 
        spectralmodes: do mm = 1, ntimeintervalsspectral
            domains: do nn = 1, ndom
 
                ! Check if the residual, time step, etc. Is needed.
 
                if (.not. resneeded) then
 
                    ! Residual, etc. Not needed.
                    ! Deallocate residual, the time step and the spectral radii
                    ! of the fine level.
 
                    deallocate (flowdoms(nn, 1, mm)%dw, flowdoms(nn, 1, mm)%fw, &
                                flowdoms(nn, 1, mm)%dtl, flowdoms(nn, 1, mm)%radI, &
                                flowdoms(nn, 1, mm)%radJ, flowdoms(nn, 1, mm)%radK, &
                                stat=ierr)
                    if (ierr /= 0) &
                         call terminate("deallocateTempMemory", &
                         "Deallocation error for dw, fw, dtl and &
                         &spectral radii.")
                end if
 
                ! The memory for the zeroth Runge Kutta stage
                ! if a Runge Kutta scheme is used.
 
                if (smoother == rungekutta) then
 
                    deallocate (flowdoms(nn, 1, mm)%wn, flowdoms(nn, 1, mm)%pn, &
                                stat=ierr)
                    if (ierr /= 0) &
                        call terminate("deallocateTempMemory", &
                                       "Deallocation error for wn and pn")
                end if
 
            end do domains
        end do spectralmodes
 
    end subroutine deallocatetempmemory
 
    subroutine allocatetempmemory(resNeeded)
        !
        !       AllocateTempMemory allocates the memory again that was
        !       temporarily deallocted by deallocateTempMemory.
        !
        use constants
        use block, only: flowdoms, ndom
        use communication, only: sendbuffer, recvbuffer, sendbuffersize, recvbuffersize
        use inputiteration, only: smoother
        use inputtimespectral, only: ntimeintervalsspectral
        use flowvarrefstate, only: nw, nwf
 
        implicit none
        !
        !      Subroutine arguments.
        !
        logical, intent(in) :: resNeeded
        !
        !      Local variables.
        !
        integer :: ierr
 
        integer(kind=intType) :: nn, mm
        integer(kind=intType) :: il, jl, kl, ie, je, ke, ib, jb, kb
 
        ! The memory for the receive buffers.
 
        allocate (sendbuffer(sendbuffersize), &
                  recvbuffer(recvbuffersize), stat=ierr)
        if (ierr /= 0) &
            call terminate("allocateTempMemory", &
                           "Memory allocation failure for comm buffers")
 
        ! Loop over the spectral modes and domains. Note that only memory
        ! on the finest mesh level needs to be reallocated, because the
        ! memory on the coarser levels has not been released or is not
        ! needed .
 
        spectralmodes: do mm = 1, ntimeintervalsspectral
            domains: do nn = 1, ndom
 
                ! Store some dimensions a bit easier.
 
                il = flowdoms(nn, 1, mm)%il
                jl = flowdoms(nn, 1, mm)%jl
                kl = flowdoms(nn, 1, mm)%kl
 
                ie = flowdoms(nn, 1, mm)%ie
                je = flowdoms(nn, 1, mm)%je
                ke = flowdoms(nn, 1, mm)%ke
 
                ib = flowdoms(nn, 1, mm)%ib
                jb = flowdoms(nn, 1, mm)%jb
                kb = flowdoms(nn, 1, mm)%kb
 
                ! Check if the residual, time step, etc. was deallocated.
 
                if (.not. resneeded) then
 
                    ! Allocate the residual, the time step and
                    ! the spectral radii.
 
                    allocate (flowdoms(nn, 1, mm)%dw(0:ib, 0:jb, 0:kb, 1:nw), &
                              flowdoms(nn, 1, mm)%fw(0:ib, 0:jb, 0:kb, 1:nwf), &
                              flowdoms(nn, 1, mm)%dtl(1:ie, 1:je, 1:ke), &
                              flowdoms(nn, 1, mm)%radI(1:ie, 1:je, 1:ke), &
                              flowdoms(nn, 1, mm)%radJ(1:ie, 1:je, 1:ke), &
                              flowdoms(nn, 1, mm)%radK(1:ie, 1:je, 1:ke), stat=ierr)
                    if (ierr /= 0) &
                         call terminate("allocateTempMemory", &
                         "Memory allocation failure for dw, fw, &
                         &dtl and the spectral radii.")
 
                    ! Initialize dw and fw to zero to avoid possible overflows
                    ! of the halo's.
 
                    flowdoms(nn, 1, mm)%dw = zero
                    flowdoms(nn, 1, mm)%fw = zero
 
                end if
 
                ! The memory for the zeroth runge kutta stage
                ! if a runge kutta scheme is used.
 
                if (smoother == rungekutta) then
 
                    allocate (flowdoms(nn, 1, mm)%wn(2:il, 2:jl, 2:kl, 1:nwf), &
                              flowdoms(nn, 1, mm)%pn(2:il, 2:jl, 2:kl), stat=ierr)
                    if (ierr /= 0) &
                        call terminate("allocateTempMemory", &
                                       "Memory allocation failure for wn and pn")
                end if
 
            end do domains
        end do spectralmodes
 
    end subroutine allocatetempmemory
 
    subroutine getliftdirfromsymmetry(liftDir)
 
        ! The purpose of this function is to determine what coordinate
        ! direction the mirror plane is in. It does NOT handle multiple mirror
        ! planes. It is used just to determine what the lift direction is.
 
        use constants
        use blockpointers, only: x, il, jl, kl, bctype, ndom, bcdata, bcfaceid, nbocos
        use communication, only: adflow_comm_world
        implicit none
 
        ! Output
        integer(kind=intType), intent(out) :: liftDir
        integer(kind=intType), dimension(3) :: sym_local, sym
 
        ! Working
        integer(kind=intType) :: nn, i_index(1), mm, ierr
        real(kind=realtype), dimension(:, :, :), pointer :: xx
        real(kind=realtype) :: cp(3), v1(3), v2(3)
        ! Loop over each block and each subFace
 
        sym_local = 0_inttype
        sym = 0_inttype
        liftdir = 0_inttype
        do nn = 1, ndom
            call setpointers(nn, 1, 1)
            do mm = 1, nbocos
                if (bctype(mm) == symm) then
 
                    select case (bcfaceid(mm))
                    case (imin)
                        xx => x(1, :, :, :)
                    case (imax)
                        xx => x(il, :, :, :)
                    case (jmin)
                        xx => x(:, 1, :, :)
                    case (jmax)
                        xx => x(:, jl, :, :)
                    case (kmin)
                        xx => x(:, :, 1, :)
                    case (kmax)
                        xx => x(:, :, kl, :)
                    end select
 
                    ! Take the cross product
                    v1(:) = xx(bcdata(mm)%inEnd, bcdata(mm)%jnEnd, :) - &
                            xx(bcdata(mm)%inBeg, bcdata(mm)%jnBeg, :)
                    v2(:) = xx(bcdata(mm)%inBeg, bcdata(mm)%jnEnd, :) - &
                            xx(bcdata(mm)%inEnd, bcdata(mm)%jnBeg, :)
 
                    ! Cross Product
                    cp(1) = (v1(2) * v2(3) - v1(3) * v2(2))
                    cp(2) = (v1(3) * v2(1) - v1(1) * v2(3))
                    cp(3) = (v1(1) * v2(2) - v1(2) * v2(1))
 
                    ! Only interesed in abs values
                    cp = abs(cp)
 
                    ! Location, ie coordiante direction of dominate direction
                    i_index = maxloc(real(cp))
 
                    sym_local(i_index(1)) = 1_inttype
                end if
            end do
        end do
 
        ! Now we have a bunch of sym_locals, mpi_allreduce them and SUM
 
        call mpi_allreduce(sym_local, sym, 3, adflow_integer, &
                           mpi_sum, adflow_comm_world, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Now we should make sure that only ONE of the values is
        ! non-zero. If more than one value is zero, it means we have
        ! multiple symmetry planes which we can't support.
        if (sym(1) == 0 .and. sym(2) == 0 .and. sym(3) == 0) then
            ! Pass - no sym, can't determine lift dir:
        else if (sym(1) .ne. 0 .and. sym(2) == 0 .and. sym(3) == 0) then
            ! Pass - x dir can't be symmetry
        else if (sym(1) == 0 .and. sym(2) .ne. 0 .and. sym(3) == 0) then
            liftdir = 3
        else if (sym(1) == 0 .and. sym(2) == 0 .and. sym(3) .ne. 0) then
            liftdir = 2
        else
            ! Multiple orientations...can't do anything
        end if
 
    end subroutine getliftdirfromsymmetry
 
    subroutine writeintromessage
        !
        !       writeIntroMessage writes a message to stdout with
        !       information how the executable was built, e.g. whether single
        !       or double precision is used for the integers and reals, etc.
        !       To avoid a messy output only processor 0 prints this info.
        !
        use constants
        use communication, only: myid, nproc
        implicit none
        !
        !      Local variables
        !
        character(len=7) :: integerString
 
        ! Return if this is not processor 0.
 
        if (myid > 0) return
 
        ! I'm processor 0. Write the info to stdout.
 
        print "(a)", "#"
        print "(a)", "# ADflow, multiblock structured flow solver"
        print "(a)", "#"
        print "(a)", "# This code solves the 3D RANS, laminar NS or &
             &Euler equations"
        print "(a)", "# on multiblock structured hexahedral grids."
 
        write (integerstring, "(i7)") nproc
        integerstring = adjustl(integerstring)
        print "(3a)", "# This is a parallel executable running on ", &
            trim(integerstring), " processors."
        print "(a)", "# It has been compiled with the &
             &following options:"
 
        if (debug) then
            print "(a)", "# - Debug mode."
        else
            print "(a)", "# - Optimized mode."
        end if
 
#ifdef USE_LONG_INT
        print "(a)", "# - Size of standard integers: 8 bytes."
#else
        print "(a)", "# - Size of standard integers: 4 bytes."
#endif
 
#ifdef USE_SINGLE_PRECISION
        print "(a)", "# - Size of standard floating point types: &
             &4 bytes."
 
#elif  USE_QUADRUPLE_PRECISION
        print "(a)", "# - Size of standard floating point types: &
             &16 bytes."
#else
        print "(a)", "# - Size of standard floating point types: &
             &8 bytes."
#endif
 
#ifdef USE_NO_CGNS
        print "(a)", "# - Without cgns support"
#else
        print "(a)", "# - With cgns support"
#endif
 
#ifdef USE_NO_SIGNALS
        print "(a)", "# - Without support for signals."
#else
        print "(a)", "# - With support for signals."
#endif
 
        print "(a)", "#"
 
    end subroutine writeintromessage
 
    subroutine pointreduce(pts, N, tol, uniquePts, link, nUnique)
 
        ! Given a list of N points (pts) in three space, with possible
        ! duplicates, (to within tol) return a list of the nUnique
        ! uniquePoints of points and a link array of length N, that points
        ! into the unique list
 
        use constants
        use kdtree2_module
        implicit none
 
        ! Input Parameters
        integer(kind=intType), intent(in) :: N
        real(kind=realtype), intent(in), dimension(3, N) :: pts
        real(kind=realtype), intent(in) :: tol
 
        ! Output Parametres
        real(kind=realtype), intent(out), dimension(3, N) :: uniquepts
        integer(kind=intType), intent(out), dimension(N) :: link
        integer(kind=intType), intent(out) :: nUnique
 
        ! Working paramters
        type(kdtree2), pointer :: mytree
        real(kind=realtype) :: tol2, timeb, timea
        integer(kind=intType) :: nFound, i, j, nAlloc
        type(kdtree2_result), allocatable, dimension(:) :: results
 
        if (n == 0) then
            nunique = 0
            return
        end if
 
        ! We will use the KD_tree to do most of the heavy lifting here:
 
        mytree => kdtree2_create(pts, sort=.true.)
 
        ! KD tree works with the square of the tolerance
        tol2 = tol**2
 
        ! Unlikely we'll have more than 20 points same, but there is a
        ! safetly check anwyay.
        nalloc = 20
        allocate (results(nalloc))
 
        link = 0
        nunique = 0
 
        ! Loop over all nodes
        do i = 1, n
            if (link(i) == 0) then
                call kdtree2_r_nearest(mytree, pts(:, i), tol2, nfound, nalloc, results)
 
                ! Expand if necesary and re-run
                if (nfound > nalloc) then
                    deallocate (results)
                    nalloc = nfound
                    allocate (results(nalloc))
                    call kdtree2_r_nearest(mytree, pts(:, i), tol2, nfound, nalloc, results)
                end if
 
                if (nfound == 1) then
                    ! This one is easy, it is already a unique node
                    nunique = nunique + 1
                    link(i) = nunique
                    uniquepts(:, nunique) = pts(:, i)
                else
                    if (link(i) == 0) then
                        ! This node hasn't been assigned yet:
                        nunique = nunique + 1
                        uniquepts(:, nunique) = pts(:, i)
 
                        do j = 1, nfound
                            link(results(j)%idx) = nunique
                        end do
                    end if
                end if
            end if
        end do
 
        ! Done with the tree and the result vector
        call kdtree2destroy(mytree)
        deallocate (results)
 
    end subroutine pointreduce
    subroutine releasememorypart1
        !
        !       releaseMemoryPart1 releases all the memory on the coarser
        !       grids of flowDoms and the fine grid memory which is not needed
        !       for the possible interpolation of the spectral solution.
        !
 
        ! This is a free-for-all on the imports. Oh well.
        use block
        use inputiteration
        use inputtimespectral
        use inputphysics
        use inputunsteady
        use monitor
        use cgnsgrid
        use communication
        use iteration
        use cgnsgrid
        use section
        use walldistancedata
        use adjointvars
        use adjointpetsc
        use surfacefamilies
        implicit none
        !
        !      Local variables
        !
        integer :: ierr
 
        integer(kind=intType) :: sps, nLevels, level, nn, l, i, j
 
        ! Determine the number of grid levels present in flowDoms.
 
        nlevels = ubound(flowdoms, 2)
 
        ! Loop over the number of spectral solutions.
 
        spectralloop: do sps = 1, ntimeintervalsspectral
 
            ! Loop over the coarser grid levels and local blocks and
            ! deallocate all the memory.
 
            do level = 2, nlevels
                do nn = 1, ndom
                    call deallocateblock(nn, level, sps)
                end do
            end do
 
            ! Release some memory of the fine grid, which is not needed
            ! anymore.
 
            do nn = 1, ndom
                ! Modified by HDN
                ! Added dwALE, fwALE
                deallocate ( &
                    flowdoms(nn, 1, sps)%dw, flowdoms(nn, 1, sps)%fw, &
                    flowdoms(nn, 1, sps)%dtl, flowdoms(nn, 1, sps)%radI, &
                    flowdoms(nn, 1, sps)%radJ, flowdoms(nn, 1, sps)%radK, &
                    flowdoms(nn, 1, sps)%shockSensor, &
                    stat=ierr)
                if (ierr /= 0) &
                     call terminate("releaseMemoryPart1", &
                     "Deallocation error for dw, fw, dwALE, fwALE, dtl and &
                     &spectral radii.")
 
                ! Extra variables for ALE
                if (equationmode == unsteady .and. useale) then
                    deallocate ( &
                        flowdoms(nn, 1, sps)%dwALE, &
                        flowdoms(nn, 1, sps)%fwALE, &
                        stat=ierr)
                    if (ierr /= 0) &
                        call terminate("releaseMemoryPart1", &
                                       "Deallocation error for dwALE, fwALE.")
                end if
 
                ! Nullify the pointers, such that no attempt is made to
                ! release the memory again.
 
                nullify (flowdoms(nn, 1, sps)%dw)
                nullify (flowdoms(nn, 1, sps)%fw)
                nullify (flowdoms(nn, 1, sps)%dwALE) ! Added by HDN
                nullify (flowdoms(nn, 1, sps)%fwALE) ! Added by HDN
                nullify (flowdoms(nn, 1, sps)%dtl)
                nullify (flowdoms(nn, 1, sps)%radI)
                nullify (flowdoms(nn, 1, sps)%radJ)
                nullify (flowdoms(nn, 1, sps)%radK)
                nullify (flowdoms(nn, 1, sps)%scratch)
                nullify (flowdoms(nn, 1, sps)%shockSensor)
                ! Check if the zeroth stage runge kutta memory has been
                ! allocated. If so deallocate it and nullify the pointers.
 
                if (smoother == rungekutta) then
 
                    deallocate (flowdoms(nn, 1, sps)%wn, flowdoms(nn, 1, sps)%pn, &
                                stat=ierr)
                    if (ierr /= 0) &
                        call terminate("releaseMemoryPart1", &
                                       "Deallocation error for wn and pn")
 
                    nullify (flowdoms(nn, 1, sps)%wn)
                    nullify (flowdoms(nn, 1, sps)%pn)
 
                end if
 
                ! Release the memory of the old residuals for the time
                ! accurate Runge-Kutta schemes.
 
                if (equationmode == unsteady .and. &
                    timeintegrationscheme == explicitrk) then
 
                    deallocate (flowdoms(nn, 1, sps)%dwOldRK, stat=ierr)
                    if (ierr /= 0) &
                        call terminate("releaseMemoryPart1", &
                                       "Deallocation error for dwOldRK,")
 
                    nullify (flowdoms(nn, 1, sps)%dwOldRK)
                end if
 
            end do
 
        end do spectralloop
 
        ! derivative values
        if (derivvarsallocated) then
            call deallocderivativevalues(1)
        end if
 
        ! Bunch of extra sutff that hasn't been deallocated
        if (allocated(cyclestrategy)) then
            deallocate (cyclestrategy)
        end if
 
        if (allocated(monnames)) then
            deallocate (monnames)
        end if
 
        if (allocated(monloc)) then
            deallocate (monloc)
        end if
 
        if (allocated(monglob)) then
            deallocate (monglob)
        end if
 
        if (allocated(monref)) then
            deallocate (monref)
        end if
 
        if (allocated(cgnsfamilies)) then
            deallocate (cgnsfamilies)
        end if
 
        if (allocated(cgnsdomsd)) then
            deallocate (cgnsdomsd)
        end if
        ! deallocate(famIDsDomainInterfaces, &
        !      bcIDsDomainInterfaces,  &
        !      famIDsSliding)
        if (allocated(sections)) then
            deallocate (sections)
        end if
 
        if (allocated(bcfamexchange)) then
            do j = 1, size(bcfamexchange, 2)
                do i = 1, size(bcfamexchange, 1)
                    call destroyfamilyexchange(bcfamexchange(i, j))
                end do
            end do
            deallocate (bcfamexchange)
        end if
 
        if (allocated(ncellglobal)) then
            deallocate (ncellglobal)
        end if
 
        ! Now deallocate the containers and communication objects.
        if (allocated(commpatterncell_1st)) then
            do l = 1, nlevels
                call deallocatecommtype(commpatterncell_1st(l))
            end do
            deallocate (commpatterncell_1st)
        end if
        if (allocated(commpatterncell_2nd)) then
            do l = 1, nlevels
                call deallocatecommtype(commpatterncell_2nd(l))
            end do
            deallocate (commpatterncell_2nd)
        end if
        if (allocated(commpatternnode_1st)) then
            do l = 1, nlevels
                call deallocatecommtype(commpatternnode_1st(l))
            end do
            deallocate (commpatternnode_1st)
        end if
        if (allocated(internalcell_1st)) then
            do l = 1, nlevels
                call deallocateinternalcommtype(internalcell_1st(l))
            end do
            deallocate (internalcell_1st)
        end if
        if (allocated(internalcell_2nd)) then
            do l = 1, nlevels
                call deallocateinternalcommtype(internalcell_2nd(l))
            end do
            deallocate (internalcell_2nd)
        end if
        if (allocated(internalnode_1st)) then
            do l = 1, nlevels
                call deallocateinternalcommtype(internalnode_1st(l))
            end do
            deallocate (internalnode_1st)
        end if
 
        ! Send/recv buffer
        if (allocated(sendbuffer)) then
            deallocate (sendbuffer)
        end if
 
        if (allocated(recvbuffer)) then
            deallocate (recvbuffer)
        end if
 
        ! massFlow stuff from setFamilyInfoFaces.f90
        if (allocated(massflowfamilyinv)) then
            deallocate (massflowfamilyinv)
        end if
        if (allocated(massflowfamilydiss)) then
            deallocate (massflowfamilydiss)
        end if
 
    end subroutine releasememorypart1
 
    subroutine deallocatecommtype(comm)
        use communication
        implicit none
        integer(kind=intType) :: ierr, i
 
        type(commtype) :: comm
        ! Deallocate memory in comm
 
        ! Deallocate the sendLists
        do i = 1, comm%nProcSend
            deallocate (comm%sendList(i)%block, stat=ierr)
            call echk(ierr, __file__, __line__)
 
            deallocate (comm%sendList(i)%indices, stat=ierr)
            call echk(ierr, __file__, __line__)
 
            deallocate (comm%sendList(i)%interp, stat=ierr)
            call echk(ierr, __file__, __line__)
        end do
 
        ! Deallocate the recvLists
        do i = 1, comm%nProcRecv
            deallocate (comm%recvList(i)%block, stat=ierr)
            call echk(ierr, __file__, __line__)
 
            deallocate (comm%recvList(i)%indices, stat=ierr)
            call echk(ierr, __file__, __line__)
        end do
 
        deallocate (comm%sendProc, stat=ierr)
        call echk(ierr, __file__, __line__)
 
        deallocate (comm%nsend, stat=ierr)
        call echk(ierr, __file__, __line__)
 
        deallocate (comm%nsendcum, stat=ierr)
        call echk(ierr, __file__, __line__)
 
        deallocate (comm%sendlist, stat=ierr)
        call echk(ierr, __file__, __line__)
 
        deallocate (comm%recvProc, stat=ierr)
        call echk(ierr, __file__, __line__)
 
        deallocate (comm%nrecv, stat=ierr)
        call echk(ierr, __file__, __line__)
 
        deallocate (comm%nrecvcum, stat=ierr)
        call echk(ierr, __file__, __line__)
 
        deallocate (comm%recvlist, stat=ierr)
        call echk(ierr, __file__, __line__)
 
        deallocate (comm%indexsendproc, stat=ierr)
        call echk(ierr, __file__, __line__)
 
        deallocate (comm%indexrecvproc, stat=ierr)
        call echk(ierr, __file__, __line__)
 
        if (comm%nPeriodic > 0) then
            do i = 1, comm%nPeriodic
                deallocate (comm%periodicData(i)%block, stat=ierr)
                call echk(ierr, __file__, __line__)
 
                deallocate (comm%periodicData(i)%indices)
                call echk(ierr, __file__, __line__)
            end do
 
            deallocate (comm%periodicData, stat=ierr)
            call echk(ierr, __file__, __line__)
        end if
 
    end subroutine deallocatecommtype
 
    subroutine deallocateinternalcommtype(comm)
        use communication
        implicit none
        integer(kind=intType) :: ierr, i
 
        type(internalcommtype) :: comm
        ! Deallocate memory in comm
        deallocate (comm%donorBlock, stat=ierr)
        call echk(ierr, __file__, __line__)
 
        deallocate (comm%donorIndices, stat=ierr)
        call echk(ierr, __file__, __line__)
 
        deallocate (comm%donorInterp, stat=ierr)
        call echk(ierr, __file__, __line__)
 
        deallocate (comm%haloBlock, stat=ierr)
        call echk(ierr, __file__, __line__)
 
        deallocate (comm%haloIndices, stat=ierr)
        call echk(ierr, __file__, __line__)
 
        if (comm%nPeriodic > 0) then
            do i = 1, comm%nPeriodic
                deallocate (comm%periodicData(i)%block, stat=ierr)
                call echk(ierr, __file__, __line__)
 
                deallocate (comm%periodicData(i)%indices)
                call echk(ierr, __file__, __line__)
            end do
            deallocate (comm%periodicData, stat=ierr)
            call echk(ierr, __file__, __line__)
        end if
 
    end subroutine deallocateinternalcommtype
 
    subroutine deallocderivativevalues(level)
 
        use constants
        use block, only: flowdomsd, flowdoms, ndom
        use inputtimespectral, only: ntimeintervalsspectral
        use walldistancedata, only: xsurfvec, xsurfvecd
        use flowvarrefstate, only: winfd
        use inputphysics, only: walldistanceneeded
        use adjointvars, only: derivvarsallocated
        use bcpointers_b
 
        implicit none
 
        ! Input Parameters
        integer(kind=intType) :: level
 
        ! Local variables
        integer(kind=intType) :: nn, sps, stat, mm, ierr
 
        do nn = 1, ndom
            do sps = 1, ntimeintervalsspectral
 
                deallocate ( &
                    flowdomsd(nn, level, sps)%x, &
                    flowdomsd(nn, level, sps)%si, &
                    flowdomsd(nn, level, sps)%sj, &
                    flowdomsd(nn, level, sps)%sk, &
                    flowdomsd(nn, level, sps)%vol, &
                    flowdomsd(nn, level, sps)%rotMatrixI, &
                    flowdomsd(nn, level, sps)%rotMatrixJ, &
                    flowdomsd(nn, level, sps)%rotMatrixK, &
                    flowdomsd(nn, level, sps)%s, &
                    flowdomsd(nn, level, sps)%sFaceI, &
                    flowdomsd(nn, level, sps)%sFaceJ, &
                    flowdomsd(nn, level, sps)%sFaceK, &
                    flowdomsd(nn, level, sps)%w, &
                    flowdomsd(nn, level, sps)%dw, &
                    flowdomsd(nn, level, sps)%fw, &
                    flowdomsd(nn, level, sps)%scratch, &
                    flowdomsd(nn, level, sps)%p, &
                    flowdomsd(nn, level, sps)%gamma, &
                    flowdomsd(nn, level, sps)%aa, &
                    flowdomsd(nn, level, sps)%rlv, &
                    flowdomsd(nn, level, sps)%rev, &
                    flowdomsd(nn, level, sps)%dtl, &
                    flowdomsd(nn, level, sps)%radI, &
                    flowdomsd(nn, level, sps)%radJ, &
                    flowdomsd(nn, level, sps)%radK, &
                    flowdomsd(nn, level, sps)%ux, &
                    flowdomsd(nn, level, sps)%uy, &
                    flowdomsd(nn, level, sps)%uz, &
                    flowdomsd(nn, level, sps)%vx, &
                    flowdomsd(nn, level, sps)%vy, &
                    flowdomsd(nn, level, sps)%vz, &
                    flowdomsd(nn, level, sps)%wx, &
                    flowdomsd(nn, level, sps)%wy, &
                    flowdomsd(nn, level, sps)%wz, &
                    flowdomsd(nn, level, sps)%qx, &
                    flowdomsd(nn, level, sps)%qy, &
                    flowdomsd(nn, level, sps)%qz, &
                    flowdomsd(nn, level, sps)%bmti1, &
                    flowdomsd(nn, level, sps)%bmti2, &
                    flowdomsd(nn, level, sps)%bmtj1, &
                    flowdomsd(nn, level, sps)%bmtj2, &
                    flowdomsd(nn, level, sps)%bmtk1, &
                    flowdomsd(nn, level, sps)%bmtk2, &
                    flowdomsd(nn, level, sps)%bvti1, &
                    flowdomsd(nn, level, sps)%bvti2, &
                    flowdomsd(nn, level, sps)%bvtj1, &
                    flowdomsd(nn, level, sps)%bvtj2, &
                    flowdomsd(nn, level, sps)%bvtk1, &
                    flowdomsd(nn, level, sps)%bvtk2, &
                    flowdomsd(nn, level, sps)%d2Wall, &
                    stat=ierr)
                call echk(ierr, __file__, __line__)
 
                ! Deallocate allocated boundary data
                do mm = 1, flowdoms(nn, level, sps)%nBocos
                    deallocate ( &
                        flowdomsd(nn, level, sps)%BCData(mm)%norm, &
                        flowdomsd(nn, level, sps)%BCData(mm)%rface, &
                        flowdomsd(nn, level, sps)%BCData(mm)%Fp, &
                        flowdomsd(nn, level, sps)%BCData(mm)%Fv, &
                        flowdomsd(nn, level, sps)%BCData(mm)%Tp, &
                        flowdomsd(nn, level, sps)%BCData(mm)%Tv, &
                        flowdomsd(nn, level, sps)%BCData(mm)%F, &
                        flowdomsd(nn, level, sps)%BCData(mm)%T, &
                        flowdomsd(nn, level, sps)%BCData(mm)%area, &
                        flowdomsd(nn, level, sps)%BCData(mm)%uSlip, &
                        flowdomsd(nn, level, sps)%BCData(mm)%TNS_Wall, &
                        stat=ierr)
                    call echk(ierr, __file__, __line__)
                end do
 
                deallocate (flowdomsd(nn, level, sps)%BCData, stat=ierr)
                call echk(ierr, __file__, __line__)
 
                viscbocoloop: do mm = 1, flowdoms(nn, level, sps)%nViscBocos
                    deallocate ( &
                        flowdomsd(nn, level, sps)%viscSubface(mm)%tau, &
                        flowdomsd(nn, level, sps)%viscSubface(mm)%q, &
                        stat=ierr)
                    call echk(ierr, __file__, __line__)
                end do viscbocoloop
 
                deallocate (flowdomsd(nn, level, sps)%viscSubFace, stat=ierr)
                call echk(ierr, __file__, __line__)
 
            end do
        end do
 
        ! Also dealloc winfd
        deallocate (winfd)
 
        ! Finally deallocate flowdomsd
        deallocate (flowdomsd, stat=ierr)
        call echk(ierr, __file__, __line__)
 
        ! And the petsc vector(s)
        if (.not. walldistanceneeded) then
            do sps = 1, ntimeintervalsspectral
                call vecdestroy(xsurfvec(1, sps), ierr)
            end do
        end if
 
        do sps = 1, ntimeintervalsspectral
            call vecdestroy(xsurfvecd(sps), ierr)
            call echk(ierr, __file__, __line__)
        end do
        deallocate (xsurfvecd)
 
        derivvarsallocated = .false.
    end subroutine deallocderivativevalues
 
    !      ---------------------------------------------------------------------------
 
    subroutine releasememorypart2
        !
        !       releaseMemoryPart2 releases all the memory of flowDoms on the
        !       finest grid as well as the memory allocated in the other
        !       modules.
        !
        use block
        use inputtimespectral
        use inputphysics, only: cpmin_family, sepsenmaxfamily
        use adjointpetsc
        use cgnsgrid
        implicit none
        !
        !      Local variables
        !
        integer :: ierr
 
        integer(kind=intType) :: nn, sps
 
        ! Release the memory of flowDoms of the finest grid and of the
        ! array flowDoms afterwards.
        if (allocated(flowdoms)) then
            do sps = 1, ntimeintervalsspectral
                do nn = 1, ndom
                    call deallocateblock(nn, 1_inttype, sps)
                end do
            end do
            deallocate (flowdoms, stat=ierr)
            if (ierr /= 0) &
                call terminate("releaseMemoryPart2", &
                               "Deallocation failure for flowDoms")
        end if
 
        ! Some more memory should be deallocated if this code is to
        ! be used in combination with adaptation.
 
        ! deallocate the cpmin_family array allocated in inputParamRoutines
        if (allocated(cpmin_family)) &
            deallocate (cpmin_family)
 
        ! deallocate the sepSenMaxFamily array allocated in inputParamRoutines
        if (allocated(sepsenmaxfamily)) &
            deallocate (sepsenmaxfamily)
 
        ! Destroy variables allocated in preprocessingAdjoint
        if (adjointpetscpreprocvarsallocated) then
            call vecdestroy(w_like1, petscierr)
            call echk(petscierr, __file__, __line__)
 
            call vecdestroy(w_like2, petscierr)
            call echk(petscierr, __file__, __line__)
 
            call vecdestroy(psi_like1, petscierr)
            call echk(petscierr, __file__, __line__)
 
            call vecdestroy(psi_like2, petscierr)
            call echk(petscierr, __file__, __line__)
 
            call vecdestroy(psi_like3, petscierr)
            call echk(petscierr, __file__, __line__)
 
            call vecdestroy(x_like, petscierr)
            call echk(petscierr, __file__, __line__)
        end if
 
        ! Finally delete cgnsDoms...but there is still more
        ! pointers that need to be deallocated...
        if (allocated(cgnsdoms)) then
            do nn = 1, cgnsndom
                if (associated(cgnsdoms(nn)%procStored)) &
                    deallocate (cgnsdoms(nn)%procStored)
 
                if (associated(cgnsdoms(nn)%conn1to1)) &
                    deallocate (cgnsdoms(nn)%conn1to1)
 
                if (associated(cgnsdoms(nn)%connNonMatchAbutting)) &
                    deallocate (cgnsdoms(nn)%connNonMatchAbutting)
 
                if (associated(cgnsdoms(nn)%bocoInfo)) &
                    deallocate (cgnsdoms(nn)%bocoInfo)
 
                deallocate ( &
                    cgnsdoms(nn)%iBegOr, cgnsdoms(nn)%iEndOr, &
                    cgnsdoms(nn)%jBegOr, cgnsdoms(nn)%jEndOr, &
                    cgnsdoms(nn)%kBegOr, cgnsdoms(nn)%kEndOr, &
                    cgnsdoms(nn)%localBlockID)
            end do
        end if
 
    end subroutine releasememorypart2
 
    subroutine deallocateblock(nn, level, sps)
        !
        !       deallocateBlock deallocates all the allocated memory of the
        !       given block.
        !
        use constants
        use inputunsteady
        use inputphysics
        use iteration
        use inputiteration, only: useskewnesscheck
        use block, only: viscsubfacetype, bcdatatype, flowdoms
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: nn, level, sps
        !
        !      Local variables.
        !
        integer :: ierr
 
        integer(kind=intType) :: i
 
        type(viscsubfacetype), dimension(:), pointer :: viscSubface
        type(bcdatatype), dimension(:), pointer :: BCData
 
        logical :: deallocationFailure
 
        ! Initialize deallocationFailure to .false.
 
        deallocationfailure = .false.
 
        ! Set the pointer for viscSubface and deallocate the memory
        ! stored in there. Initialize ierr to 0, such that the terminate
        ! routine is only called at the end if a memory deallocation
        ! failure occurs.
        ierr = 0
        viscsubface => flowdoms(nn, level, sps)%viscSubface
        do i = 1, flowdoms(nn, level, sps)%nViscBocos
            deallocate (viscsubface(i)%tau, viscsubface(i)%q, &
                        viscsubface(i)%utau, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            nullify (viscsubface(i)%tau)
            nullify (viscsubface(i)%q)
            nullify (viscsubface(i)%utau)
        end do
 
        ! Set the pointer for BCData and deallocate the memory
        ! stored in there.
        bcdata => flowdoms(nn, level, sps)%BCData
        do i = 1, flowdoms(nn, level, sps)%nBocos
 
            if (associated(bcdata(i)%norm)) &
                deallocate (bcdata(i)%norm, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(bcdata(i)%area)) &
                deallocate (bcdata(i)%area, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(bcdata(i)%surfIndex)) &
                deallocate (bcdata(i)%surfIndex, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(bcdata(i)%F)) &
                deallocate (bcdata(i)%F, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(bcdata(i)%Fv)) &
                deallocate (bcdata(i)%Fv, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(bcdata(i)%Fp)) &
                deallocate (bcdata(i)%Fp, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(bcdata(i)%T)) &
                deallocate (bcdata(i)%T, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(bcdata(i)%Tv)) &
                deallocate (bcdata(i)%Tv, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(bcdata(i)%Tp)) &
                deallocate (bcdata(i)%Tp, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(bcdata(i)%rface)) &
                deallocate (bcdata(i)%rface, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(bcdata(i)%uSlip)) &
                deallocate (bcdata(i)%uSlip, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(bcdata(i)%TNS_Wall)) &
                deallocate (bcdata(i)%TNS_Wall, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(bcdata(i)%ptInlet)) &
                deallocate (bcdata(i)%ptInlet, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(bcdata(i)%ttInlet)) &
                deallocate (bcdata(i)%ttInlet, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(bcdata(i)%htInlet)) &
                deallocate (bcdata(i)%htInlet, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(bcdata(i)%flowXdirInlet)) &
                deallocate (bcdata(i)%flowXdirInlet, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(bcdata(i)%flowYdirInlet)) &
                deallocate (bcdata(i)%flowYdirInlet, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(bcdata(i)%flowZdirInlet)) &
                deallocate (bcdata(i)%flowZdirInlet, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(bcdata(i)%rho)) &
                deallocate (bcdata(i)%rho, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(bcdata(i)%velx)) &
                deallocate (bcdata(i)%velx, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(bcdata(i)%vely)) &
                deallocate (bcdata(i)%vely, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(bcdata(i)%velz)) &
                deallocate (bcdata(i)%velz, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(bcdata(i)%ps)) &
                deallocate (bcdata(i)%ps, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(bcdata(i)%turbInlet)) &
                deallocate (bcdata(i)%turbInlet, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(bcdata(i)%normALE)) &
                deallocate (bcdata(i)%normALE, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
            if (associated(bcdata(i)%rFaceALE)) &
                deallocate (bcdata(i)%rFaceALE, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
            if (associated(bcdata(i)%uSlipALE)) &
                deallocate (bcdata(i)%uSlipALE, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
            if (associated(bcdata(i)%cellHeatFlux)) &
                deallocate (bcdata(i)%cellHeatFlux, stat=ierr)
            if (associated(bcdata(i)%nodeHeatFlux)) &
                deallocate (bcdata(i)%nodeHeatFlux, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(bcdata(i)%iBlank)) &
                deallocate (bcdata(i)%iBlank, stat=ierr)
 
            if (ierr /= 0) deallocationfailure = .true.
 
            nullify (bcdata(i)%norm)
            nullify (bcdata(i)%rface)
            nullify (bcdata(i)%F)
            nullify (bcdata(i)%Fv)
            nullify (bcdata(i)%Fp)
            nullify (bcdata(i)%T)
            nullify (bcdata(i)%Tv)
            nullify (bcdata(i)%Tp)
 
            nullify (bcdata(i)%uSlip)
            nullify (bcdata(i)%TNS_Wall)
 
            nullify (bcdata(i)%normALE)
            nullify (bcdata(i)%rfaceALE)
            nullify (bcdata(i)%uSlipALE)
            nullify (bcdata(i)%cellHeatFlux)
            nullify (bcdata(i)%nodeHeatFlux)
 
            nullify (bcdata(i)%ptInlet)
            nullify (bcdata(i)%ttInlet)
            nullify (bcdata(i)%htInlet)
            nullify (bcdata(i)%flowXdirInlet)
            nullify (bcdata(i)%flowYdirInlet)
            nullify (bcdata(i)%flowZdirInlet)
 
            nullify (bcdata(i)%turbInlet)
 
            nullify (bcdata(i)%rho)
            nullify (bcdata(i)%velx)
            nullify (bcdata(i)%vely)
            nullify (bcdata(i)%velz)
            nullify (bcdata(i)%ps)
            nullify (bcdata(i)%iblank)
 
        end do
 
        if (associated(flowdoms(nn, level, sps)%BCType)) &
            deallocate (flowdoms(nn, level, sps)%BCType, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%BCFaceID)) &
            deallocate (flowdoms(nn, level, sps)%BCFaceID, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%cgnsSubface)) &
            deallocate (flowdoms(nn, level, sps)%cgnsSubface, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%inBeg)) &
            deallocate (flowdoms(nn, level, sps)%inBeg, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%inEnd)) &
            deallocate (flowdoms(nn, level, sps)%inEnd, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%jnBeg)) &
            deallocate (flowdoms(nn, level, sps)%jnBeg, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%jnEnd)) &
            deallocate (flowdoms(nn, level, sps)%jnEnd, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%knBeg)) &
            deallocate (flowdoms(nn, level, sps)%knBeg, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%knEnd)) &
            deallocate (flowdoms(nn, level, sps)%knEnd, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%dinBeg)) &
            deallocate (flowdoms(nn, level, sps)%dinBeg, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%dinEnd)) &
            deallocate (flowdoms(nn, level, sps)%dinEnd, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%djnBeg)) &
            deallocate (flowdoms(nn, level, sps)%djnBeg, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%djnEnd)) &
            deallocate (flowdoms(nn, level, sps)%djnEnd, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%dknBeg)) &
            deallocate (flowdoms(nn, level, sps)%dknBeg, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%dknEnd)) &
            deallocate (flowdoms(nn, level, sps)%dknEnd, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%icBeg)) &
            deallocate (flowdoms(nn, level, sps)%icBeg, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%icEnd)) &
            deallocate (flowdoms(nn, level, sps)%icEnd, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%jcBeg)) &
            deallocate (flowdoms(nn, level, sps)%jcBeg, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%jcEnd)) &
            deallocate (flowdoms(nn, level, sps)%jcEnd, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%kcBeg)) &
            deallocate (flowdoms(nn, level, sps)%kcBeg, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%kcEnd)) &
            deallocate (flowdoms(nn, level, sps)%kcEnd, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%neighBlock)) &
            deallocate (flowdoms(nn, level, sps)%neighBlock, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%neighProc)) &
            deallocate (flowdoms(nn, level, sps)%neighProc, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%l1)) &
            deallocate (flowdoms(nn, level, sps)%l1, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%l2)) &
            deallocate (flowdoms(nn, level, sps)%l2, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%l3)) &
            deallocate (flowdoms(nn, level, sps)%l3, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%groupNum)) &
            deallocate (flowdoms(nn, level, sps)%groupNum, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%iblank)) &
            deallocate (flowdoms(nn, level, sps)%iblank, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%forcedRecv)) &
            deallocate (flowdoms(nn, level, sps)%forcedRecv, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%status)) &
            deallocate (flowdoms(nn, level, sps)%status, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%BCData)) &
            deallocate (flowdoms(nn, level, sps)%BCData, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%viscSubface)) &
            deallocate (flowdoms(nn, level, sps)%viscSubface, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%viscIminPointer)) &
            deallocate (flowdoms(nn, level, sps)%viscIminPointer, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%viscImaxPointer)) &
            deallocate (flowdoms(nn, level, sps)%viscImaxPointer, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%viscJminPointer)) &
            deallocate (flowdoms(nn, level, sps)%viscJminPointer, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%viscJmaxPointer)) &
            deallocate (flowdoms(nn, level, sps)%viscJmaxPointer, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%viscKminPointer)) &
            deallocate (flowdoms(nn, level, sps)%viscKminPointer, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%viscKmaxPointer)) &
            deallocate (flowdoms(nn, level, sps)%viscKmaxPointer, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%x)) &
            deallocate (flowdoms(nn, level, sps)%x, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%xOld)) &
            deallocate (flowdoms(nn, level, sps)%xOld, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%si)) &
            deallocate (flowdoms(nn, level, sps)%si, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%sj)) &
            deallocate (flowdoms(nn, level, sps)%sj, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%sk)) &
            deallocate (flowdoms(nn, level, sps)%sk, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%vol)) &
            deallocate (flowdoms(nn, level, sps)%vol, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (useskewnesscheck) then
            if (associated(flowdoms(nn, level, sps)%skew)) &
                deallocate (flowdoms(nn, level, sps)%skew, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
        end if
 
        if (associated(flowdoms(nn, level, sps)%volRef)) &
            deallocate (flowdoms(nn, level, sps)%volRef, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%volOld)) &
            deallocate (flowdoms(nn, level, sps)%volOld, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%pori)) &
            deallocate (flowdoms(nn, level, sps)%pori, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%porj)) &
            deallocate (flowdoms(nn, level, sps)%porj, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%pork)) &
            deallocate (flowdoms(nn, level, sps)%pork, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%indFamilyI)) &
            deallocate (flowdoms(nn, level, sps)%indFamilyI, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%indFamilyJ)) &
            deallocate (flowdoms(nn, level, sps)%indFamilyJ, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%indFamilyK)) &
            deallocate (flowdoms(nn, level, sps)%indFamilyK, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%factFamilyI)) &
            deallocate (flowdoms(nn, level, sps)%factFamilyI, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%factFamilyJ)) &
            deallocate (flowdoms(nn, level, sps)%factFamilyJ, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%factFamilyK)) &
            deallocate (flowdoms(nn, level, sps)%factFamilyK, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%rotMatrixI)) &
            deallocate (flowdoms(nn, level, sps)%rotMatrixI, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%rotMatrixJ)) &
            deallocate (flowdoms(nn, level, sps)%rotMatrixJ, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%rotMatrixK)) &
            deallocate (flowdoms(nn, level, sps)%rotMatrixK, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%sFaceI)) &
            deallocate (flowdoms(nn, level, sps)%sFaceI, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%sFaceJ)) &
            deallocate (flowdoms(nn, level, sps)%sFaceJ, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%sFaceK)) &
            deallocate (flowdoms(nn, level, sps)%sFaceK, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%w)) &
            deallocate (flowdoms(nn, level, sps)%w, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%wOld)) &
            deallocate (flowdoms(nn, level, sps)%wOld, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%p)) &
            deallocate (flowdoms(nn, level, sps)%p, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%aa)) &
            deallocate (flowdoms(nn, level, sps)%aa, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%gamma)) &
            deallocate (flowdoms(nn, level, sps)%gamma, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%ux)) &
            deallocate (flowdoms(nn, level, sps)%ux, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%uy)) &
            deallocate (flowdoms(nn, level, sps)%uy, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%uz)) &
            deallocate (flowdoms(nn, level, sps)%uz, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%vx)) &
            deallocate (flowdoms(nn, level, sps)%vx, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%vy)) &
            deallocate (flowdoms(nn, level, sps)%vy, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%vz)) &
            deallocate (flowdoms(nn, level, sps)%vz, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%wx)) &
            deallocate (flowdoms(nn, level, sps)%wx, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%wy)) &
            deallocate (flowdoms(nn, level, sps)%wy, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%wz)) &
            deallocate (flowdoms(nn, level, sps)%wz, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%qx)) &
            deallocate (flowdoms(nn, level, sps)%qx, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%qy)) &
            deallocate (flowdoms(nn, level, sps)%qy, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%qz)) &
            deallocate (flowdoms(nn, level, sps)%qz, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%rlv)) &
            deallocate (flowdoms(nn, level, sps)%rlv, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%rev)) &
            deallocate (flowdoms(nn, level, sps)%rev, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%s)) &
            deallocate (flowdoms(nn, level, sps)%s, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%p1)) &
            deallocate (flowdoms(nn, level, sps)%p1, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%dw)) &
            deallocate (flowdoms(nn, level, sps)%dw, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%fw)) &
            deallocate (flowdoms(nn, level, sps)%fw, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%dwOldRK)) &
            deallocate (flowdoms(nn, level, sps)%dwOldRK, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%w1)) &
            deallocate (flowdoms(nn, level, sps)%w1, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%wr)) &
            deallocate (flowdoms(nn, level, sps)%wr, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%mgIFine)) &
            deallocate (flowdoms(nn, level, sps)%mgIFine, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%mgJFine)) &
            deallocate (flowdoms(nn, level, sps)%mgJFine, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%mgKFine)) &
            deallocate (flowdoms(nn, level, sps)%mgKFine, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%mgIWeight)) &
            deallocate (flowdoms(nn, level, sps)%mgIWeight, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%mgJWeight)) &
            deallocate (flowdoms(nn, level, sps)%mgJWeight, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%mgKWeight)) &
            deallocate (flowdoms(nn, level, sps)%mgKWeight, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%mgICoarse)) &
            deallocate (flowdoms(nn, level, sps)%mgICoarse, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%mgJCoarse)) &
            deallocate (flowdoms(nn, level, sps)%mgJCoarse, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%mgKCoarse)) &
            deallocate (flowdoms(nn, level, sps)%mgKCoarse, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%iCo)) &
            deallocate (flowdoms(nn, level, sps)%iCo, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%jCo)) &
            deallocate (flowdoms(nn, level, sps)%jCo, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%kCo)) &
            deallocate (flowdoms(nn, level, sps)%kCo, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%wn)) &
            deallocate (flowdoms(nn, level, sps)%wn, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%pn)) &
            deallocate (flowdoms(nn, level, sps)%pn, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%dtl)) &
            deallocate (flowdoms(nn, level, sps)%dtl, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%radI)) &
            deallocate (flowdoms(nn, level, sps)%radI, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%radJ)) &
            deallocate (flowdoms(nn, level, sps)%radJ, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%radK)) &
            deallocate (flowdoms(nn, level, sps)%radK, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%d2Wall)) &
            deallocate (flowdoms(nn, level, sps)%d2Wall, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%bmti1)) &
            deallocate (flowdoms(nn, level, sps)%bmti1, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%bmti2)) &
            deallocate (flowdoms(nn, level, sps)%bmti2, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%bmtj1)) &
            deallocate (flowdoms(nn, level, sps)%bmtj1, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%bmtj2)) &
            deallocate (flowdoms(nn, level, sps)%bmtj2, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%bmtk1)) &
            deallocate (flowdoms(nn, level, sps)%bmtk1, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%bmtk2)) &
            deallocate (flowdoms(nn, level, sps)%bmtk2, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%bvti1)) &
            deallocate (flowdoms(nn, level, sps)%bvti1, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%bvti2)) &
            deallocate (flowdoms(nn, level, sps)%bvti2, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%bvtj1)) &
            deallocate (flowdoms(nn, level, sps)%bvtj1, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%bvtj2)) &
            deallocate (flowdoms(nn, level, sps)%bvtj2, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%bvtk1)) &
            deallocate (flowdoms(nn, level, sps)%bvtk1, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%bvtk2)) &
            deallocate (flowdoms(nn, level, sps)%bvtk2, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%globalCell)) &
            deallocate (flowdoms(nn, level, sps)%globalCell, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (associated(flowdoms(nn, level, sps)%globalNode)) &
            deallocate (flowdoms(nn, level, sps)%globalNode, stat=ierr)
        if (ierr /= 0) deallocationfailure = .true.
 
        if (equationmode == unsteady .and. useale) then
 
            ! Added by HDN
            if (associated(flowdoms(nn, level, sps)%xALE)) &
                deallocate (flowdoms(nn, level, sps)%xALE, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(flowdoms(nn, level, sps)%sIALE)) &
                deallocate (flowdoms(nn, level, sps)%sIALE, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(flowdoms(nn, level, sps)%sJALE)) &
                deallocate (flowdoms(nn, level, sps)%sJALE, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(flowdoms(nn, level, sps)%sKALE)) &
                deallocate (flowdoms(nn, level, sps)%sKALE, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(flowdoms(nn, level, sps)%sVeloIALE)) &
                deallocate (flowdoms(nn, level, sps)%sVeloIALE, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(flowdoms(nn, level, sps)%sVeloJALE)) &
                deallocate (flowdoms(nn, level, sps)%sVeloJALE, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(flowdoms(nn, level, sps)%sVeloKALE)) &
                deallocate (flowdoms(nn, level, sps)%sVeloKALE, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(flowdoms(nn, level, sps)%sFaceIALE)) &
                deallocate (flowdoms(nn, level, sps)%sFaceIALE, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(flowdoms(nn, level, sps)%sFaceJALE)) &
                deallocate (flowdoms(nn, level, sps)%sFaceJALE, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            if (associated(flowdoms(nn, level, sps)%sFaceKALE)) &
                deallocate (flowdoms(nn, level, sps)%sFaceKALE, stat=ierr)
            if (ierr /= 0) deallocationfailure = .true.
 
            ! if( associated(flowDoms(nn,level,sps)%dwALE) ) &
            !      deallocate(flowDoms(nn,level,sps)%dwALE, stat=ierr)
            ! if(ierr /= 0) deallocationFailure = .true.
            !
            ! if( associated(flowDoms(nn,level,sps)%fwALE) ) &
            !      deallocate(flowDoms(nn,level,sps)%fwALE, stat=ierr)
            ! if(ierr /= 0) deallocationFailure = .true.
        end if
 
        ! Check for errors in the deallocation.
 
        if (deallocationfailure) &
            call terminate("deallocateBlock", &
                           "Something went wrong when deallocating memory")
 
        ! Nullify the pointers of this block.
        call nullifyflowdompointers(nn, level, sps)
 
    end subroutine deallocateblock
 
    integer function setcgnsrealtype()
        !
        !       setCGNSRealType sets the cgns real type, depending on the
        !       compiler options. Note that quadrupole precision is not
        !       supported by CGNS; double precision is used instead for the
        !       CGNS IO.
        !
        use su_cgns, only: realsingle, realdouble
        implicit none
 
#ifdef USE_NO_CGNS
 
        call terminate("setCGNSRealType", &
             "Function should not be called if no cgns support &
             &is selected.")
 
#else
 
# ifdef USE_SINGLE_PRECISION
        setcgnsrealtype = realsingle
# else
        setcgnsrealtype = realdouble
# endif
 
#endif
 
    end function setcgnsrealtype
 
    subroutine returnfail(routineName, errorMessage)
        !
        !       returnFail writes an error message to standard output and
        !       sets fail flags to be returned to python.
        !
        use constants
        use communication, only: adflow_comm_world, myid
#ifndef USE_TAPENADE
        use killsignals, only: fatalfail, frompython, routinefailed
#endif
        implicit none
        !
        !      Subroutine arguments
        !
        character(len=*), intent(in) :: routineName
        character(len=*), intent(in) :: errorMessage
#ifndef USE_TAPENADE
 
        !
        !      Local parameter
        !
        integer, parameter :: maxCharLine = 55
        !
        !      Local variables
        !
        integer :: ierr, len, i2
        logical :: firstTime
 
        character(len=len_trim(errorMessage)) :: message
        character(len=8) :: integerString
 
        ! Copy the errorMessage into message. It is not possible to work
        ! with errorMessage directly, because it is modified in this
        ! routine. Sometimes a constant string is passed to this routine
        ! and some compilers simply fail then.
 
        message = errormessage
 
        ! Print a nice error message. In case of a parallel executable
        ! also the processor id is printed.
 
        print "(a)", "#"
        print "(a)", "#--------------------------- !!! Error !!! &
             &----------------------------"
 
        write (integerstring, "(i8)") myid
        integerstring = adjustl(integerstring)
 
        print "(2a)", "#* returnFail called by processor ", &
            trim(integerstring)
 
        ! Write the header of the error message.
 
        print "(2a)", "#* Run-time error in procedure ", &
            trim(routinename)
 
        ! Loop to write the error message. If the message is too long it
        ! is split over several lines.
 
        firsttime = .true.
        do
            ! Determine the remaining error message to be written.
            ! If longer than the maximum number of characters allowed
            ! on a line, it is attempted to split the message.
 
            message = adjustl(message)
            len = len_trim(message)
            i2 = min(maxcharline, len)
 
            if (i2 < len) i2 = index(message(:i2), " ", .true.) - 1
            if (i2 < 0) i2 = index(message, " ") - 1
            if (i2 < 0) i2 = len
 
            ! Write this part of the error message. If it is the first
            ! line of the message some additional stuff is printed.
 
            if (firsttime) then
                print "(2a)", "#* Error message: ", &
                    trim(message(:i2))
                firsttime = .false.
            else
                print "(2a)", "#*                ", &
                    trim(message(:i2))
            end if
 
            ! Exit the loop if the entire message has been written.
 
            if (i2 == len) exit
 
            ! Adapt the string for the next part to be written.
 
            message = message(i2 + 1:)
 
        end do
 
        ! Write the trailing message.
 
        print "(a)", "#*"
        print "(a)", "#------------------------------------------&
             &----------------------------"
        print "(a)", "#"
 
        ! Call abort and stop the program. This stop should be done in
        ! abort, but just to be sure.
 
        if (frompython) then
            routinefailed = .true.
            fatalfail = .true.
        else
            call mpi_abort(adflow_comm_world, 1, ierr)
            stop
        end if
#endif
 
    end subroutine returnfail
 
    subroutine echk(errorcode, file, line)
 
        ! Check if ierr that resulted from a petsc or MPI call is in fact an
        ! error.
        use constants
        use communication, only: adflow_comm_world, myid
        implicit none
 
        integer(kind=intType), intent(in) :: errorcode
        character(len=*), intent(in) :: file
        integer(kind=intType), intent(in) :: line
        integer :: ierr
        character(len=maxStringLen) :: errorCodeFormat, errorLineFormat
 
        errorcodeformat = "(2(A, I2))"
        errorlineformat = "(A, I5, A, A)"
 
        if (errorcode == 0) then
            return ! No error, return immediately
        else
#ifndef USE_TAPENADE
            print *, '---------------------------------------------------------------------------'
            print errorcodeformat, "PETSc or MPI Error. Error Code ", errorcode, ". Detected on Proc ", myid
            print errorlineformat, "Error at line: ", line, " in file: ", file
            print *, '---------------------------------------------------------------------------'
            call mpi_abort(adflow_comm_world, errorcode, ierr)
            stop ! Just in case
#else
            stop
#endif
        end if
 
    end subroutine echk
 
    subroutine converttolowercase(string)
        !
        !       convertToLowerCase converts the given string to lower case.
        !
        use constants
        implicit none
        !
        !      Subroutine arguments
        !
        character(len=*), intent(inout) :: string
        !
        !      Local variables
        !
        integer(kind=intType), parameter :: upperToLower = iachar("a") - iachar("A")
 
        integer(kind=intType) :: i, lenString
 
        ! Determine the length of the given string and convert the upper
        ! case characters to lower case.
 
        lenstring = len_trim(string)
        do i = 1, lenstring
            if ("A" <= string(i:i) .and. string(i:i) <= "Z") &
                string(i:i) = achar(iachar(string(i:i)) + uppertolower)
        end do
 
    end subroutine converttolowercase
 
    logical function eulerwallspresent()
 
        ! eulerWallsPresent determines whether or not inviscid walls are
        ! present in the whole grid. It first determines if these walls are
        ! present locally and performs an allReduce afterwards.
 
        use constants
        use block, only: ndom, flowdoms
        use communication, only: adflow_comm_world
        implicit none
        !
        !      Local variables.
        !
        integer(kind=intType) :: nn, i
        integer :: ierr
        logical :: localeulerwalls
 
        ! Initialize localEulerWalls to .false. and loop over the
        ! boundary subfaces of the blocks to see if Euler walls are
        ! present on this processor. As the info is the same for all
        ! spectral solutions, only the 1st needs to be considered.
 
        localeulerwalls = .false.
        do nn = 1, ndom
            do i = 1, flowdoms(nn, 1, 1)%nBocos
                if (flowdoms(nn, 1, 1)%BCType(i) == eulerwall) &
                    localeulerwalls = .true.
            end do
        end do
 
        ! Set i to 1 if Euler walls are present locally and to 0
        ! otherwise. Determine the maximum over all processors
        ! and set EulerWallsPresent accordingly.
 
        i = 0
        if (localeulerwalls) i = 1
        call mpi_allreduce(i, nn, 1, adflow_integer, mpi_max, &
                           adflow_comm_world, ierr)
 
        if (nn == 0) then
            eulerwallspresent = .false.
        else
            eulerwallspresent = .true.
        end if
 
    end function eulerwallspresent
    subroutine allocconvarrays(nIterTot)
        !
        !       allocConvArrays allocates the memory for the convergence
        !       arrays. The number of iterations allocated, nIterTot, is
        !       enough to store the maximum number of iterations specified
        !       plus possible earlier iterations read from the restart file.
        !       This routine MAY be called with data already inside of
        !       convArray and this will be saved.
        !
        use constants
        use inputtimespectral, only: ntimeintervalsspectral
        use inputio, only: storeconvinneriter
        use monitor, only: convarray, nmon, solverdataarray, solvertypearray, showcpu
        implicit none
        !
        !      Subroutine argument.
        !
        integer(kind=intType) :: nIterTot
        !
        !      Local variables.
        !
        integer :: ierr
 
        integer(kind=intType) :: nSolverMon ! number of solver monitor variables
 
        ! Return immediately if the convergence history (of the inner
        ! iterations) does not need to be stored. This logical can
        ! only be .false. for an unsteady computation.
        if (.not. storeconvinneriter) return
 
        if (allocated(convarray)) then
            deallocate (convarray)
        end if
        if (allocated(solverdataarray)) then
            deallocate (solverdataarray)
        end if
        if (allocated(solvertypearray)) then
            deallocate (solvertypearray)
        end if
 
        if (showcpu) then
            nsolvermon = 5
        else
            nsolvermon = 4
        end if
 
        allocate (convarray(0:nitertot, ntimeintervalsspectral, nmon))
        allocate (solverdataarray(0:nitertot, ntimeintervalsspectral, nsolvermon))
        allocate (solvertypearray(0:nitertot, ntimeintervalsspectral))
 
        ! Zero Array:
        convarray = zero
        solverdataarray = zero
 
    end subroutine allocconvarrays
 
    subroutine alloctimearrays(nTimeTot)
        !
        !       allocTimeArrays allocates the memory for the arrays to store
        !       the time history of the unsteady computation. The number of
        !       time steps specified is enought to store the total number of
        !       time steps of the current computation plus possible earlier
        !       computations.
        !
        use constants
        use monitor, only: timearray, timedataarray, nmon
        implicit none
        !
        !      Subroutine argument.
        !
        integer(kind=intType) :: nTimeTot
        !
        !      Local variables.
        !
        integer :: ierr
 
        ! Allocate the memory for both the time array as well as the
        ! data array.
 
        if (allocated(timearray)) then
            deallocate (timearray)
        end if
        if (allocated(timedataarray)) then
            deallocate (timedataarray)
        end if
 
        allocate (timearray(ntimetot), &
                  timedataarray(ntimetot, nmon), stat=ierr)
        if (ierr /= 0) &
             call terminate("allocTimeArrays", &
             "Memory allocation failure for timeArray &
             &and timeDataArray")
 
    end subroutine alloctimearrays
 
    subroutine getmonitorvariablenames(nvar, monitor_variables)
        !
        !  copy the names in monnames to another array so that is can be
        !  passed back up the python level
        !
        use constants
        use monitor, only: nmon, monnames
        implicit none
 
        ! save the monitor variable names into a new array
        integer(kind=intType), intent(in) :: nvar
        character, dimension(nvar, maxCGNSNameLen), intent(out) :: monitor_variables
 
        ! working variables
        character(len=maxCGNSNameLen) :: var_name
        integer(kind=intType) :: c, idx_mon
 
        do idx_mon = 1, nvar
            var_name = monnames(idx_mon)
 
            do c = 1, len(monnames(idx_mon))
                monitor_variables(idx_mon, c) = var_name(c:c)
            end do
        end do
 
    end subroutine getmonitorvariablenames
 
    subroutine getsolvertypearray(niter, nsps, type_array)
        !
        !  copy the names in sovlerTypeArray to another array so that is can be
        !  passed back up the python level
        !
        use constants
        use monitor, only: solvertypearray
        use inputtimespectral, only: ntimeintervalsspectral
        use iteration, only: itertot
        implicit none
 
        ! save the monitor variable names into a new array
        integer(kind=intType), intent(in) :: niter, nsps
        character, dimension(0:niter, ntimeintervalsspectral, maxIterTypelen), intent(out) :: type_array
 
        ! working variables
        character(len=maxIterTypelen) :: type_name
        integer(kind=intType) :: c, idx_sps, idx_iter
 
        do idx_sps = 1, ntimeintervalsspectral
            do idx_iter = 0, itertot
                type_name = solvertypearray(idx_iter, idx_sps)
 
                do c = 1, len(solvertypearray(idx_iter, idx_sps))
                    type_array(idx_iter, idx_sps, c) = type_name(c:c)
                end do
 
            end do
        end do
 
    end subroutine getsolvertypearray
 
    subroutine convergenceheader
        !
        !       convergenceHeader writes the convergence header to stdout.
        !
        use cgnsnames
        use inputphysics
        use inputunsteady
        use flowvarrefstate
        use monitor
        use iteration
        use inputiteration
        implicit none
        !
        !      Local variables.
        !
        integer(kind=intType) :: i, nCharWrite
        logical :: writeIterations
 
        ! Determine whether or not the iterations must be written.
 
        if (printiterations) then
            writeiterations = .true.
            if (equationmode == unsteady .and. &
                timeintegrationscheme == explicitrk) writeiterations = .false.
 
            ! Determine the number of characters to write.
            ! First initialize this number with the variables which are
            ! always written. This depends on the equation mode. For unsteady
            ! and spectral computations a bit more info is written.
 
            ncharwrite = 10
            if (writeiterations) ncharwrite = ncharwrite + 7 + 7 + 9 + 7 + 7 + 10
            if (equationmode == unsteady) then
                ncharwrite = ncharwrite + 7 + fieldwidth + 1
            else if (equationmode == timespectral) then
                ncharwrite = ncharwrite + 11
            end if
 
            ! Add the number of characters needed for the actual variables.
 
#ifndef USE_COMPLEX
            ! for the real version this is easy
            ncharwrite = ncharwrite + nmon * (fieldwidthlarge + 1)
#else
            ! for complex we need to differentiate between residuals and functionals
            do i = 1, nmon
                select case (monnames(i))
 
                case (cgnsl2resrho, cgnsl2resmomx, &
                      cgnsl2resmomy, cgnsl2resmomz, &
                      cgnsl2resrhoe, cgnsl2resnu, &
                      cgnsl2resk, cgnsl2resomega, &
                      cgnsl2restau, cgnsl2resepsilon, &
                      cgnsl2resv2, cgnsl2resf, 'totalR')
 
                    ! complex residuals need 9 more characters
                    ncharwrite = ncharwrite + fieldwidthlarge + 1 + 9
 
                case default
                    ! complex functionals need 25 more characters
                    ncharwrite = ncharwrite + fieldwidthlarge + 1 + 25
 
                end select
            end do
#endif
 
            if (showcpu) ncharwrite = ncharwrite + fieldwidth + 1
 
            ! Write the line of - signs. This line starts with a #, such
            ! that it is ignored by some plotting software.
 
            write (*, "(a)", advance="no") "#"
            do i = 2, ncharwrite
                write (*, "(a)", advance="no") "-"
            end do
            print "(1x)"
 
            ! Write the first line of the header. First the variables that
            ! will always be written. Some extra variables must be written
            ! for unsteady and time spectral problems.
            write (*, '("# ")', advance="no")
            write (*, "(a)", advance="no") " Grid  |"
 
            if (equationmode == unsteady) then
                write (*, "(a)", advance="no") " Time |    Time    |"
            else if (equationmode == timespectral) then
                write (*, "(a)", advance="no") " Spectral |"
            end if
 
            if (writeiterations) write (*, "(a)", advance="no") " Iter | Iter |  Iter  |   CFL   | Step | Lin  |"
            if (showcpu) write (*, "(a)", advance="no") "    Wall    |"
 
            ! Write the header for the variables to be monitored.
            do i = 1, nmon
                ! Determine the variable name and write the
                ! corresponding text.
 
                ! we do the real and complex versions separately
#ifndef USE_COMPLEX
                ! real versions print the full 16 digits so these spacings are "regular" and the same for all
                select case (monnames(i))
 
                case ("totalR")
                    write (*, "(a)", advance="no") "        totalRes        |"
 
                case (cgnsl2resrho)
                    write (*, "(a)", advance="no") "        Res rho         |"
 
                case (cgnsl2resmomx)
                    write (*, "(a)", advance="no") "        Res rhou        |"
 
                case (cgnsl2resmomy)
                    write (*, "(a)", advance="no") "        Res rhov        |"
 
                case (cgnsl2resmomz)
                    write (*, "(a)", advance="no") "        Res rhow        |"
 
                case (cgnsl2resrhoe)
                    write (*, "(a)", advance="no") "        Res rhoE        |"
 
                case (cgnsl2resnu)
                    write (*, "(a)", advance="no") "       Res nuturb       |"
 
                case (cgnsl2resk)
                    write (*, "(a)", advance="no") "       Res kturb        |"
 
                case (cgnsl2resomega)
                    write (*, "(a)", advance="no") "       Res wturb        |"
 
                case (cgnsl2restau)
                    write (*, "(a)", advance="no") "       Res tauturb      |"
 
                case (cgnsl2resepsilon)
                    write (*, "(a)", advance="no") "       Res epsturb      |"
 
                case (cgnsl2resv2)
                    write (*, "(a)", advance="no") "       Res v2turb       |"
 
                case (cgnsl2resf)
                    write (*, "(a)", advance="no") "       Res fturb        |"
 
                case (cgnscl)
                    write (*, "(a)", advance="no") "         C_lift         |"
 
                case (cgnsclp)
                    write (*, "(a)", advance="no") "        C_lift_p        |"
 
                case (cgnsclv)
                    write (*, "(a)", advance="no") "        C_lift_v        |"
 
                case (cgnscd)
                    write (*, "(a)", advance="no") "        C_drag          |"
 
                case (cgnscdp)
                    write (*, "(a)", advance="no") "        C_drag_p        |"
 
                case (cgnscdv)
                    write (*, "(a)", advance="no") "        C_drag_v        |"
 
                case (cgnscfx)
                    write (*, "(a)", advance="no") "          C_Fx          |"
 
                case (cgnscfy)
                    write (*, "(a)", advance="no") "          C_Fy          |"
 
                case (cgnscfz)
                    write (*, "(a)", advance="no") "          C_Fz          |"
 
                case (cgnscmx)
                    write (*, "(a)", advance="no") "          C_Mx          |"
 
                case (cgnscmy)
                    write (*, "(a)", advance="no") "          C_My          |"
 
                case (cgnscmz)
                    write (*, "(a)", advance="no") "          C_Mz          |"
 
                case (cgnshdiffmax)
                    write (*, "(a)", advance="no") "       |H-H_inf|        |"
 
                case (cgnsmachmax)
                    write (*, "(a)", advance="no") "        Mach_max        |"
 
                case (cgnsyplusmax)
                    write (*, "(a)", advance="no") "         Y+_max         |"
 
                case (cgnseddymax)
                    write (*, "(a)", advance="no") "        Eddyv_max       |"
 
                case (cgnssepsensor)
                    write (*, "(a)", advance="no") "        SepSensor       |"
 
                case (cgnssepsensorksarea)
                    write (*, "(a)", advance="no") "       sepSensorKsArea    |"
 
                case (cgnscavitation)
                    write (*, "(a)", advance="no") "       Cavitation       |"
 
                case (cgnsaxismoment)
                    write (*, "(a)", advance="no") "       AxisMoment       |"
 
                end select
#else
                ! complex versions print the full 16 digits for real and complex for "functionals"
                ! but shorter versions only for residuals
                select case (monnames(i))
 
                case ("totalR")
                    write (*, "(a)", advance="no") "            totalRes             |"
 
                case (cgnsl2resrho)
                    write (*, "(a)", advance="no") "            Res rho              |"
 
                case (cgnsl2resmomx)
                    write (*, "(a)", advance="no") "            Res rhou             |"
 
                case (cgnsl2resmomy)
                    write (*, "(a)", advance="no") "            Res rhov             |"
 
                case (cgnsl2resmomz)
                    write (*, "(a)", advance="no") "            Res rhow             |"
 
                case (cgnsl2resrhoe)
                    write (*, "(a)", advance="no") "            Res rhoE             |"
 
                case (cgnsl2resnu)
                    write (*, "(a)", advance="no") "           Res nuturb            |"
 
                case (cgnsl2resk)
                    write (*, "(a)", advance="no") "           Res kturb             |"
 
                case (cgnsl2resomega)
                    write (*, "(a)", advance="no") "           Res wturb             |"
 
                case (cgnsl2restau)
                    write (*, "(a)", advance="no") "           Res tauturb           |"
 
                case (cgnsl2resepsilon)
                    write (*, "(a)", advance="no") "           Res epsturb           |"
 
                case (cgnsl2resv2)
                    write (*, "(a)", advance="no") "           Res v2turb            |"
 
                case (cgnsl2resf)
                    write (*, "(a)", advance="no") "           Res fturb             |"
 
                case (cgnscl)
                    write (*, "(a)", advance="no") "                     C_lift                      |"
 
                case (cgnsclp)
                    write (*, "(a)", advance="no") "                    C_lift_p                     |"
 
                case (cgnsclv)
                    write (*, "(a)", advance="no") "                    C_lift_v                     |"
 
                case (cgnscd)
                    write (*, "(a)", advance="no") "                    C_drag                       |"
 
                case (cgnscdp)
                    write (*, "(a)", advance="no") "                    C_drag_p                     |"
 
                case (cgnscdv)
                    write (*, "(a)", advance="no") "                    C_drag_v                     |"
 
                case (cgnscfx)
                    write (*, "(a)", advance="no") "                      C_Fx                       |"
 
                case (cgnscfy)
                    write (*, "(a)", advance="no") "                      C_Fy                       |"
 
                case (cgnscfz)
                    write (*, "(a)", advance="no") "                      C_Fz                       |"
 
                case (cgnscmx)
                    write (*, "(a)", advance="no") "                      C_Mx                       |"
 
                case (cgnscmy)
                    write (*, "(a)", advance="no") "                      C_My                       |"
 
                case (cgnscmz)
                    write (*, "(a)", advance="no") "                      C_Mz                       |"
 
                case (cgnshdiffmax)
                    write (*, "(a)", advance="no") "                   |H-H_inf|                     |"
 
                case (cgnsmachmax)
                    write (*, "(a)", advance="no") "                    Mach_max                     |"
 
                case (cgnsyplusmax)
                    write (*, "(a)", advance="no") "                     Y+_max                      |"
 
                case (cgnseddymax)
                    write (*, "(a)", advance="no") "                    Eddyv_max                    |"
 
                case (cgnssepsensor)
                    write (*, "(a)", advance="no") "                    SepSensor                    |"
 
                case (cgnssepsensorksarea)
                    write (*, "(a)", advance="no") "                    sepSensorKsArea                |"
 
                case (cgnscavitation)
                    write (*, "(a)", advance="no") "                   Cavitation                    |"
 
                case (cgnsaxismoment)
                    write (*, "(a)", advance="no") "                   AxisMoment                    |"
 
                end select
#endif
 
            end do
 
            print "(1x)"
 
            ! Write the second line of the header. Most of them are empty,
            ! but some variables require a second line.
            write (*, '("# ")', advance="no")
 
            write (*, "(a)", advance="no") " level |"
 
            if (equationmode == unsteady) then
                write (*, "(a)", advance="no") " Step |            |"
            else if (equationmode == timespectral) then
                write (*, "(a)", advance="no") " Solution |"
            end if
 
            if (writeiterations) write (*, "(a)", advance="no") "      | Tot  |  Type  |         |      | Res  |"
            if (showcpu) write (*, "(a)", advance="no") " Clock (s)  |"
 
            ! Loop over the variables to be monitored and write the
            ! second line.
 
            do i = 1, nmon
 
                ! Determine the variable name and write the
                ! corresponding text.
#ifndef USE_COMPLEX
                ! real mode gets the same width for all variables
                select case (monnames(i))
 
                case (cgnshdiffmax)
                    write (*, "(a)", advance="no") "           max          |"
 
                case default
                    write (*, "(a)", advance="no") "                        |"
 
                end select
#else
                ! complex mode gets shorter spacing for residuals
                select case (monnames(i))
 
                case (cgnshdiffmax)
                    write (*, "(a)", advance="no") "                       max                       |"
 
                case (cgnsl2resrho, cgnsl2resmomx, &
                      cgnsl2resmomy, cgnsl2resmomz, &
                      cgnsl2resrhoe, cgnsl2resnu, &
                      cgnsl2resk, cgnsl2resomega, &
                      cgnsl2restau, cgnsl2resepsilon, &
                      cgnsl2resv2, cgnsl2resf, 'totalR')
                    ! residuals get a shorter line
                    write (*, "(a)", advance="no") "                                 |"
                case default
                    ! regular functionals get the long empty line
                    write (*, "(a)", advance="no") "                                                 |"
 
                end select
#endif
 
            end do
            print "(1x)"
        end if
 
        ! Write again a line of - signs (starting with a #).
 
        write (*, "(a)", advance="no") "#"
        do i = 2, ncharwrite
            write (*, "(a)", advance="no") "-"
        end do
        print "(1x)"
 
    end subroutine convergenceheader
    subroutine sumresiduals(nn, mm)
        !
        !       sumResiduals adds the sum of the residuals squared at
        !       position nn to the array monLoc at position mm. It is assumed
        !       that the arrays of blockPointers already point to the correct
        !       block.
        !
        use blockpointers
        use monitor
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: nn, mm
        !
        !      Local variables.
        !
        integer(kind=intType) :: i, j, k
 
        ! Loop over the number of owned cells of this block and
        ! accumulate the residual.
 
        do k = 2, kl
            do j = 2, jl
                do i = 2, il
#ifndef USE_COMPLEX
                    monloc(mm) = monloc(mm) + (dw(i, j, k, nn) / vol(i, j, k))**2
#else
                    ! TODO squaring the complex residual when its order 1e-200 underflows and we need a better approach here
                    ! we need to square and sum the real and complex parts separately
                    monloc(mm) = monloc(mm) + &
                                 cmplx((real(dw(i, j, k, nn) / vol(i, j, k)))**2, &
                                       (aimag(dw(i, j, k, nn) / vol(i, j, k)))**2)
#endif
                end do
            end do
        end do
 
    end subroutine sumresiduals
 
    subroutine sumallresiduals(mm)
        !
        !       sumAllResiduals adds the sum of the ALL residuals squared at
        !       to monLoc at position mm.
        !
        use blockpointers
        use monitor
        use flowvarrefstate
        use inputiteration
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: mm
        !
        !      Local variables.
        !
        integer(kind=intType) :: i, j, k, l
        real(kind=realtype) :: state_sum, ovv
 
        ! Loop over the number of owned cells of this block and
        ! accumulate the residual.
 
        do k = 2, kl
            do j = 2, jl
                do i = 2, il
                    state_sum = 0.0
                    ovv = one / vol(i, j, k)
                    do l = 1, nwf
#ifndef USE_COMPLEX
                        state_sum = state_sum + (dw(i, j, k, l) * ovv)**2
#else
                        ! TODO squaring the complex residual when its order 1e-200 underflows and we need a better approach here
                        ! we need to square and sum the real and complex parts separately
                        state_sum = state_sum + &
                                    cmplx((real(dw(i, j, k, l) * ovv))**2, &
                                          (aimag(dw(i, j, k, l) * ovv))**2)
#endif
                    end do
                    do l = nt1, nt2
                        ! l-nt1+1 will index the turbResScale properly
#ifndef USE_COMPLEX
                        state_sum = state_sum + (dw(i, j, k, l) * ovv * turbresscale(l - nt1 + 1))**2
#else
                        ! we need to square and sum the real and complex parts separately
                        state_sum = state_sum + &
                                    cmplx((real(dw(i, j, k, l) * ovv * turbresscale(l - nt1 + 1)))**2, &
                                          (aimag(dw(i, j, k, l) * ovv * turbresscale(l - nt1 + 1)))**2)
#endif
                    end do
                    monloc(mm) = monloc(mm) + state_sum
                end do
            end do
        end do
 
    end subroutine sumallresiduals
 
    subroutine unsteadyheader
        !
        !       unsteadyHeader writes a header to stdout when a new time step
        !       is started.
        !
        use constants
        use monitor, only: ntimestepsrestart, timeunsteadyrestart, timeunsteady, timestepunsteady
        use commonformats, only: strings
        implicit none
        !
        !      Local variables
        !
        character(len=7) :: integerString
        character(len=12) :: realString
 
        ! Write the time step number to the integer string and the
        ! physical time to the real string.
 
        write (integerstring, "(i7)") timestepunsteady + ntimestepsrestart
        write (realstring, "(es12.5)") timeunsteady + timeunsteadyrestart
 
        integerstring = adjustl(integerstring)
        realstring = adjustl(realstring)
 
        ! Write the header to stdout.
 
        print "(a)", "#"
        print "(a)", "#**************************************************************************"
        print "(a)", "#"
        print strings, "# Unsteady time step ", trim(integerstring), ", physical time ", trim(realstring), " seconds"
        print "(a)", "#"
        print "(a)", "#**************************************************************************"
        print "(a)", "#"
 
    end subroutine unsteadyheader
 
    subroutine getcellcenters(level, n, xCen)
 
        use constants
        use inputtimespectral, only: ntimeintervalsspectral
        use blockpointers, only: ndom, il, jl, kl, x
 
        implicit none
 
        ! Input/Output
        integer(kind=intType), intent(in) :: level, n
        real(kind=realtype), dimension(3, n), intent(out) :: xcen
 
        ! Working
        integer(kind=intType) :: i, j, k, ii, nn, sps
 
        ii = 0
        do nn = 1, ndom
            do sps = 1, ntimeintervalsspectral
                call setpointers(nn, level, sps)
 
                do k = 2, kl
                    do j = 2, jl
                        do i = 2, il
                            ii = ii + 1
 
                            xcen(:, ii) = eighth * ( &
                                          x(i - 1, j - 1, k - 1, :) + &
                                          x(i, j - 1, k - 1, :) + &
                                          x(i - 1, j, k - 1, :) + &
                                          x(i, j, k - 1, :) + &
                                          x(i - 1, j - 1, k, :) + &
                                          x(i, j - 1, k, :) + &
                                          x(i - 1, j, k, :) + &
                                          x(i, j, k, :))
                        end do
                    end do
                end do
            end do
        end do
    end subroutine getcellcenters
 
    subroutine getcellcgnsblockids(level, n, cellID)
 
        use constants
        use inputtimespectral, only: ntimeintervalsspectral
        use blockpointers, only: ndom, il, jl, kl, nbkglobal
 
        implicit none
 
        ! Input/Output
        integer(kind=intType), intent(in) :: level, n
        real(kind=realtype), dimension(n), intent(out) :: cellid
 
        ! Working
        integer(kind=intType) :: i, j, k, ii, nn, sps
 
        ii = 0
        do nn = 1, ndom
            do sps = 1, ntimeintervalsspectral
                call setpointers(nn, level, sps)
 
                do k = 2, kl
                    do j = 2, jl
                        do i = 2, il
                            ii = ii + 1
 
                            cellid(ii) = nbkglobal
 
                        end do
                    end do
                end do
            end do
        end do
    end subroutine getcellcgnsblockids
 
    subroutine getncgnszones(nZones)
        use cgnsgrid
        implicit none
        integer(kind=inttype), intent(out) :: nZones
 
        nzones = cgnsndom
 
    end subroutine getncgnszones
 
    subroutine getcgnszonename(i, zone)
        use cgnsgrid
        implicit none
        character(len=maxCGNSNameLen), intent(out) :: zone
        integer(kind=intType), intent(in) :: i
 
        zone = cgnsdoms(i)%zoneName
 
    end subroutine getcgnszonename
 
    subroutine getrotmatrix(vec1, vec2, rotMat)
        !
        !       getRotmatrix computes and returns a rotation matrix that rotates
        !       the direction of vec1 to the direction of vec2. vec1 and vec2 are
        !       arrays of size 3, and the routine returns the dense 3 by 3 rotMat
        !       matrix. The vectors don't need to be normalized.
        !
        use constants
 
        implicit none
 
        real(kind=realtype), dimension(3), intent(in) :: vec1, vec2
        real(kind=realtype), dimension(3, 3), intent(out) :: rotmat
 
        integer(kind=intType) :: i
        real(kind=realtype), dimension(3) :: vec3
        real(kind=realtype), dimension(3, 3) :: wmat, wmat2
        real(kind=realtype) :: theta, dotprod, mag1, mag2, cosine
 
        ! first compute the angle between the two velocity vectors
        ! Compute the dot product of vec1 and vec2
        dotprod = zero
        do i = 1, 3
            dotprod = dotprod + vec1(i) * vec2(i)
        end do
 
        ! Compute the magnitude of vec1 and vec2
        mag1 = mynorm2(vec1)
        mag2 = mynorm2(vec2)
 
        ! Compute the cosine of the angle between vec1 and vec2
        cosine = dotprod / (mag1 * mag2)
 
        ! Compute the angle in radians
        theta = acos(cosine)
 
        ! then compute the normal vector of the rotation plane. This catches all changes in
        ! alpha and beta, and we don't need to keep track of the individual rotations and apply
        ! them in the same order.
 
        ! Compute the cross product of vec1 and vec2 to get vec3
        vec3(1) = vec1(2) * vec2(3) - vec1(3) * vec2(2)
        vec3(2) = vec1(3) * vec2(1) - vec1(1) * vec2(3)
        vec3(3) = vec1(1) * vec2(2) - vec1(2) * vec2(1)
        ! normalize
        vec3 = vec3 / sqrt(sum(vec3**2))
 
        ! compute the rotation matrix.
        ! for this, we use the "Rodrigues' Rotation Formula" as described here:
        ! https://mathworld.wolfram.com/RodriguesRotationFormula.html
        wmat = zero
        wmat2 = zero
        rotmat = zero
 
        ! compute Wmat and Wmat^2 first
        ! W is the multiplication of a 3D levi-civita tensor with vec3.
        ! we could do this in a loop but it takes more lines of code
        ! and is more difficult to understand. so just set the values
        wmat(1, 2) = -vec3(3)
        wmat(1, 3) = vec3(2)
        wmat(2, 1) = vec3(3)
        wmat(2, 3) = -vec3(1)
        wmat(3, 1) = -vec3(2)
        wmat(3, 2) = vec3(1)
        wmat2 = matmul(wmat, wmat)
 
        ! start with identity
        do i = 1, 3
            rotmat(i, i) = one
        end do
 
        ! add the W terms
        rotmat = rotmat + sin(theta) * wmat + (1.0 - cos(theta)) * wmat2
 
    end subroutine getrotmatrix
 
    real(kind=realtype) function norm2cplx(v)
        ! if input is real:
        !     returns the norm of the input array
        !
        ! if input is complex:
        !     returns a complex number where the real part represents the norm of
        !     all the real input parts and the complex part represents the norm of
        !     all the complex input parts.
 
        use constants
 
        implicit none
 
        real(kind=realtype), dimension(:), intent(in) :: v
 
#ifdef USE_COMPLEX
        norm2cplx = cmplx(norm2(real(v)), norm2(aimag(v)))
#else
        norm2cplx = norm2(v)
#endif
 
    end function norm2cplx
 
#endif
end module utils