variableReading.F90

Go to the documentation of this file.

module variablereading
 
    use constants, only: inttype, maxcgnsnamelen, cgnsrealtype, realtype, maxstringlen
    use su_cgns, only: cgsize_t
    ! halosRead:Determines if the halos where read or not.
    logical :: halosread
 
    ! cgnsInd:  File index of the CGNS file.
    ! cgnsBase: Base of the CGNS file, always set to 1.
    ! cgnsZone: Zone ID in the CGNS file.
    ! cgnsSol:  Solution ID in the zone ID of the CGNS file.
    ! location: Location where the variables are stored in CGNS.
    !           Supported possibilities are Vertex and CellCentered.
 
    integer :: cgnsind, cgnsbase, cgnszone, cgnssol, location
 
    ! zoneNumbers: Corresponding zoneNumbers of the sorted
    !              zoneNames.
    ! zoneNames:   Zone names, sorted in increasing order, of the
    !              zones in the CGNS restart file.
    ! varNames: Variable names, sorted in increasing order,
    !           of the variables.
 
    integer(kind=intType), allocatable, dimension(:) :: zonenumbers
    character(len=maxCGNSNameLen), allocatable, dimension(:) :: zonenames
    character(len=maxCGNSNameLen), allocatable, dimension(:) :: varnames
 
    ! rangeMin(3):    Lower index in i, j and k direction of the
    !                 range to be read.
    ! rangeMax(3):    Upper index in i, j and k direction of the
    !                 range to be read.
    integer(kind=cgsize_t), dimension(3) :: rangemin, rangemax
 
    ! nVar:     Number of variables stored in the solution file.
    ! solID:    Loop variables for the number of solutions to be read.
 
    integer :: nvar
    integer(kind=intType) :: solid
 
    ! interpolSpectral:     Whether or not to interpolate the
    !                       coordinates/solutions for the time
    !                       spectral mode.
    ! copySpectral:         Whether or not to copy the solutions
    !                       for the time spectral mode.
    logical :: interpolspectral, copyspectral
 
    ! rhoScale: Scale factor for the density.
    ! velScale: Scale factor for the velocity.
    ! pScale:   Scale factor for the pressure.
    ! muScale:  Scale factor for the molecular viscosity.
 
    real(kind=realtype) :: rhoscale, velscale, pscale, muscale
 
    ! nSolsRead:            Number of solution files to read.
    ! solFiles(nSolsRead):  Names of the solution files to be read.
 
    integer(kind=intType) :: nsolsread
    character(len=maxStringLen), dimension(:), allocatable :: solfiles
 
    ! varTypes(nVar): Variable types of the variables stored.
    integer, allocatable, dimension(:) :: vartypes
 
    ! buffer(2:il,2:jl,2:kl): Buffer to read and store the cell
    !                         centered values.
    ! bufferVertex(:):        Additional buffer needed to read
    !                         vertex data and transform them into
    !                         cell centered values.
 
    real(kind=cgnsrealtype), dimension(:, :, :), allocatable :: buffer
    real(kind=cgnsrealtype), dimension(:, :, :), allocatable :: buffervertex
 
contains
    subroutine readdensity(nTypeMismatch)
        !
        !       readDensity reads the density from the given place in the
        !       cgns file. If the density itself is not stored (unlikely),
        !       then it is tried to construct the density from other
        !       variables. If this is not possible an error message is printed
        !       and the program will stop.
        !       It is assumed that the pointers in blockPointers already
        !       point to the correct block.
        !
        use constants
        use cgnsnames
        use blockpointers, only: w, nbklocal
        use iomodule, only: iovar
        use utils, only: setcgnsrealtype, terminate
        use sorting, only: bsearchstrings
        implicit none
        !
        !      Subroutine argument.
        !
        integer(kind=intType), intent(inout) :: nTypeMismatch
        !
        !      Local variables
        !
        integer :: realTypeCGNS
 
        integer(kind=intType) :: i, j, k, nn, po, ip, jp, kp
        integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, kBeg, kEnd
 
        ! Set the cell range to be copied from the buffer.
 
        ibeg = lbound(buffer, 1); iend = ubound(buffer, 1)
        jbeg = lbound(buffer, 2); jend = ubound(buffer, 2)
        kbeg = lbound(buffer, 3); kend = ubound(buffer, 3)
 
        ! Set the cgns real type and abbreviate the solution variable and
        ! the pointer offset to improve readability.
 
        realtypecgns = setcgnsrealtype()
 
        po = iovar(nbklocal, solid)%pointerOffset
        w => iovar(nbklocal, solid)%w
 
        ! Find out if the density is present in the solution file.
 
        nn = bsearchstrings(cgnsdensity, varnames)
        if (nn > 0) then
 
            ! Density is present. First determine whether or not a type
            ! mismatch occurs. If so, update nTypeMismatch.
 
            if (realtypecgns /= vartypes(nn)) &
                ntypemismatch = ntypemismatch + 1
 
            ! Read the density from the restart file and store it in buffer.
 
            call readrestartvariable(varnames(nn))
 
            ! Copy the variables from buffer into w. Multiply by the
            ! scaling factor to obtain to correct nondimensional value and
            ! take the possible pointer offset into account.
 
            do k = kbeg, kend
                kp = k + po
                do j = jbeg, jend
                    jp = j + po
                    do i = ibeg, iend
                        ip = i + po
                        w(ip, jp, kp, irho) = buffer(i, j, k) * rhoscale
                    end do
                end do
            end do
 
            ! Density is read, so a return can be made.
 
            return
 
        end if
 
        ! Not able to determine the density.
        ! Print an error message and exit.
 
        call terminate("readDensity", &
             "Not able to retrieve density from the &
             &variables in the restart file.")
 
    end subroutine readdensity
 
    subroutine readenergy(nTypeMismatch)
        !
        !       readEnergy reads the energy variable from the given place in
        !       the cgns file. If the energy is not stored then it is tried to
        !       construct it from the pressure, density and velocities. If it
        !       is not possible to create the energy an error message is
        !       printed and the program will stop. It is assumed that the
        !       pointers in blockPointers already point to the correct block.
        !
        use constants
        use cgnsnames
        use blockpointers, only: w, nbklocal
        use iomodule, only: iovar
        use utils, only: setcgnsrealtype, terminate
        use sorting, only: bsearchstrings
        use flowutils, only: etot
        use flowvarrefstate, only: kpresent
        implicit none
        !
        !      Subroutine argument.
        !
        integer(kind=intType), intent(inout) :: nTypeMismatch
        !
        !      Local variables
        !
        integer :: realTypeCGNS
 
        integer(kind=intType) :: i, j, k, nn, po, ip, jp, kp
        integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, kBeg, kEnd
 
        real(kind=realtype) :: vvx, vvy, vvz, dummyk, pres, rhoinv
 
        ! Set the cell range to be copied from the buffer.
 
        ibeg = lbound(buffer, 1); iend = ubound(buffer, 1)
        jbeg = lbound(buffer, 2); jend = ubound(buffer, 2)
        kbeg = lbound(buffer, 3); kend = ubound(buffer, 3)
 
        ! Set the cgns real type and abbreviate the solution variable and
        ! the pointer offset to improve readability.
 
        realtypecgns = setcgnsrealtype()
 
        po = iovar(nbklocal, solid)%pointerOffset
        w => iovar(nbklocal, solid)%w
 
        ! Find out if the total energy is present in the solution file.
 
        nn = bsearchstrings(cgnsenergy, varnames)
 
        testrhoepresent: if (nn > 0) then
 
            ! Total energy is present. First determine whether or not a type
            ! mismatch occurs. If so, update nTypeMismatch.
 
            if (realtypecgns /= vartypes(nn)) &
                ntypemismatch = ntypemismatch + 1
 
            ! Read the energy from the restart file and store it in buffer.
 
            call readrestartvariable(varnames(nn))
 
            ! Copy the variables from buffer into w. Multiply by the scaling
            ! factor to obtain to correct non-dimensional value and take the
            ! possible pointer offset into account.
 
            do k = kbeg, kend
                kp = k + po
                do j = jbeg, jend
                    jp = j + po
                    do i = ibeg, iend
                        ip = i + po
                        w(ip, jp, kp, irhoe) = buffer(i, j, k) * pscale
                    end do
                end do
            end do
 
            ! Energy has been read, so a return can be made.
 
            return
 
        end if testrhoepresent
 
        ! Total energy is not present. Check for the pressure.
 
        nn = bsearchstrings(cgnspressure, varnames)
 
        testpressure: if (nn > 0) then
 
            ! Pressure is present. First determine whether or not a type
            ! mismatch occurs. If so, update nTypeMismatch.
 
            if (realtypecgns /= vartypes(nn)) &
                ntypemismatch = ntypemismatch + 1
 
            ! Read the pressure from the restart file and store it in buffer.
 
            call readrestartvariable(varnames(nn))
 
            ! Compute the total energy. This depends whether or not
            ! a turbulent kinetic energy is present. Take the possible
            ! pointer offset into account.
            ! As this routine is only called to construct the states in
            ! the past for a time accurate computation, the momentum is
            ! stored and not the velocity.
 
            if (kpresent) then
 
                do k = kbeg, kend
                    kp = k + po
                    do j = jbeg, jend
                        jp = j + po
                        do i = ibeg, iend
                            ip = i + po
                            rhoinv = one / w(ip, jp, kp, irho)
                            vvx = w(ip, jp, kp, imx) * rhoinv
                            vvy = w(ip, jp, kp, imy) * rhoinv
                            vvz = w(ip, jp, kp, imz) * rhoinv
                            pres = buffer(i, j, k) * pscale
                            call etot(w(ip, jp, kp, irho), vvx, vvy, vvz, pres, &
                                      w(ip, jp, kp, itu1), w(ip, jp, kp, irhoe), &
                                      kpresent)
                        end do
                    end do
                end do
 
            else
 
                dummyk = zero
 
                do k = kbeg, kend
                    kp = k + po
                    do j = jbeg, jend
                        jp = j + po
                        do i = ibeg, iend
                            ip = i + po
                            rhoinv = one / w(ip, jp, kp, irho)
                            vvx = w(ip, jp, kp, imx) * rhoinv
                            vvy = w(ip, jp, kp, imy) * rhoinv
                            vvz = w(ip, jp, kp, imz) * rhoinv
                            pres = buffer(i, j, k) * pscale
                            call etot(w(ip, jp, kp, irho), vvx, vvy, vvz, pres, &
                                      dummyk, w(ip, jp, kp, irhoe), kpresent)
                        end do
                    end do
                end do
 
            end if
 
            ! Energy has been created. So a return can be made.
 
            return
 
        end if testpressure
 
        ! Energy could not be created. Terminate.
 
        call terminate("readEnergy", &
                       "Energy could not be created")
 
    end subroutine readenergy
 
    subroutine readpressure(nTypeMismatch)
        !
        !       readPressure reads the pressure variable from the given place
        !       in the cgns file. If the pressure itself is not present it is
        !       tried to construct if from other variables. In that case it is
        !       assumed that the density, velocity and turbulent variables are
        !       already stored in the pointer variable w.
        !       If it is not possible to create the pressure an error message
        !       is printed and the program will stop.
        !       It is assumed that the pointers in blockPointers already
        !       point to the correct block.
        !
        use constants
        use cgnsnames
        use blockpointers, only: w, nbklocal
        use iomodule, only: iovar
        use utils, only: setcgnsrealtype, terminate
        use flowutils, only: computepressure
        use sorting, only: bsearchstrings
        implicit none
        !
        !      Subroutine argument.
        !
        integer(kind=intType), intent(inout) :: nTypeMismatch
        !
        !      Local variables
        !
        integer :: realTypeCGNS
 
        integer(kind=intType) :: i, j, k, nn, po, ip, jp, kp
        integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, kBeg, kEnd
 
        ! Set the cell range to be copied from the buffer.
 
        ibeg = lbound(buffer, 1); iend = ubound(buffer, 1)
        jbeg = lbound(buffer, 2); jend = ubound(buffer, 2)
        kbeg = lbound(buffer, 3); kend = ubound(buffer, 3)
 
        ! Set the cgns real type and abbreviate the solution variable and
        ! the pointer offset to improve readability.
 
        realtypecgns = setcgnsrealtype()
 
        po = iovar(nbklocal, solid)%pointerOffset
        w => iovar(nbklocal, solid)%w
 
        ! Find out if the pressure is present in the solution file.
 
        nn = bsearchstrings(cgnspressure, varnames)
        if (nn > 0) then
 
            ! Pressure is present. First determine whether or not a type
            ! mismatch occurs. If so, update nTypeMismatch.
 
            if (realtypecgns /= vartypes(nn)) &
                ntypemismatch = ntypemismatch + 1
 
            ! Read the pressure from the restart file and store
            ! it in buffer.
 
            call readrestartvariable(varnames(nn))
 
            ! Copy the variables from buffer into the position of rhoE
            ! in w. Multiply by the pressure scale factor to obtain the
            ! correct nondimensional value and take the possible pointer
            ! offset into account.
 
            do k = kbeg, kend
                kp = k + po
                do j = jbeg, jend
                    jp = j + po
                    do i = ibeg, iend
                        ip = i + po
                        w(ip, jp, kp, irhoe) = buffer(i, j, k) * pscale
                    end do
                end do
            end do
 
            ! Pressure is read, so a return can be made.
 
            return
 
        end if
 
        ! Pressure is not present. Check for the total energy.
 
        nn = bsearchstrings(cgnsenergy, varnames)
        if (nn > 0) then
 
            ! Total energy is present. First determine whether or not a type
            ! mismatch occurs. If so, update nTypeMismatch.
 
            if (realtypecgns /= vartypes(nn)) &
                ntypemismatch = ntypemismatch + 1
 
            ! Read the total energy from the restart file and store
            ! it in buffer.
 
            call readrestartvariable(varnames(nn))
 
            ! Copy the variables from buffer into w. Multiply by the
            ! pressure scale factor to obtain the correct nondimensional
            ! value and take the possible pointer offset into account.
 
            do k = kbeg, kend
                kp = k + po
                do j = jbeg, jend
                    jp = j + po
                    do i = ibeg, iend
                        ip = i + po
                        w(ip, jp, kp, irhoe) = buffer(i, j, k) * pscale
                    end do
                end do
            end do
 
            ! Compute the pressure from energy, density and velocities.
            ! This will still be stored in the irhoE position of w.
 
            call computepressure(ibeg, iend, jbeg, jend, kbeg, kend, po)
 
            ! Pressure is constructed, so a return can be made.
 
            return
 
        end if
 
        ! Not able to determine the pressure.
        ! Print an error message and exit.
 
        call terminate("readPressure", &
             "Not able to retrieve the pressure from &
             &the variables in the restart file.")
 
    end subroutine readpressure
 
    subroutine readturbeddyvis(nTypeMismatch, eddyVisPresent)
        !
        !       readTurbEddyVis tries to read the eddy viscosity from the
        !       restart file.
        !
        use constants
        use cgnsnames
        use blockpointers, only: w, nbklocal, il, jl, kl, rev
        use utils, only: setcgnsrealtype, terminate
        use sorting, only: bsearchstrings
        use flowvarrefstate, only: muinf
        use inputphysics, only: eddyvisinfratio
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(inout) :: nTypeMismatch
        logical, intent(out) :: eddyVisPresent
        !
        !      Local variables.
        !
        integer :: realTypeCGNS
 
        integer(kind=intType) :: i, j, k, nn
 
        ! Set the cgns real type
 
        realtypecgns = setcgnsrealtype()
 
        ! Check if the eddy viscosity is present. If so, read it.
 
        nn = bsearchstrings(cgnseddy, varnames)
 
        if (nn > 0) then
 
            ! Eddy viscosity is present. Determine if a type mismatch
            ! occured, read the buffer from the file and set
            ! eddyVisPresent to .true.
 
            if (realtypecgns /= vartypes(nn)) &
                ntypemismatch = ntypemismatch + 1
 
            call readrestartvariable(varnames(nn))
 
            eddyvispresent = .true.
 
            ! Scale the eddy viscosity such that it contains the
            ! correct nonDimensional value.
 
            do k = 2, kl
                do j = 2, jl
                    do i = 2, il
                        rev(i, j, k) = muscale * buffer(i, j, k)
                    end do
                end do
            end do
 
            ! Eddy viscosity has been read, so make a return.
 
            return
 
        end if
 
        ! Eddy viscosity is not present. Check if the eddy viscosity
        ! ratio is present. If so read it and construct the eddy
        ! viscosity from it.
 
        nn = bsearchstrings(cgnseddyratio, varnames)
 
        if (nn > 0) then
 
            ! Eddy viscosity ratio is present. Determine if a type
            ! mismatch occured, read the buffer from the file and set
            ! eddyVisPresent to .true.
 
            if (realtypecgns /= vartypes(nn)) &
                ntypemismatch = ntypemismatch + 1
 
            call readrestartvariable(varnames(nn))
 
            eddyvispresent = .true.
 
            ! Multiply the eddy viscosity by the laminar viscosity such
            ! that it contains the correct nonDimensional value.
            ! As the laminar viscosity is not yet know, use the free
            ! stream viscosity.
 
            do k = 2, kl
                do j = 2, jl
                    do i = 2, il
                        rev(i, j, k) = muinf * buffer(i, j, k)
                    end do
                end do
            end do
 
            ! Eddy viscosity has been read, so make a return.
 
            return
 
        end if
 
        ! Eddy viscosity cannot be constructed. Set it to the
        ! free stream eddy viscosity.
 
        do k = 2, kl
            do j = 2, jl
                do i = 2, il
                    rev(i, j, k) = muinf * eddyvisinfratio
                end do
            end do
        end do
 
        ! Eddy viscosity is not present, so set it to .false.
 
        eddyvispresent = .false.
 
    end subroutine readturbeddyvis
 
    subroutine readturbkwtype(nTypeMismatch)
        !
        !       readTurbKwType reads or constructs the k and omega values
        !       for two-equations turbulence models of the k-omega type.
        !       If no information could be retrieved some engineering guess of
        !       the turbulent variables is made.
        !
        use constants
        use cgnsnames
        use communication, only: myid
        use blockpointers, only: w, nbklocal, il, jl, kl, rlv, rev
        use iomodule, only: iovar
        use utils, only: setcgnsrealtype, terminate
        use sorting, only: bsearchstrings
        use flowvarrefstate, only: muinf
        use inputphysics, only: turbmodel
        use turbutils, only: initkomega
        implicit none
        !
        !      Subroutine argument.
        !
        integer(kind=intType), intent(inout) :: nTypeMismatch
        !
        !      Local variables.
        !
        integer :: realTypeCGNS
 
        integer(kind=intType) :: i, j, k, nn, po, ip, jp, kp
        integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, kBeg, kEnd
 
        real(kind=realtype) :: nuscale, kscale, omegascale, val
 
        logical :: turbKPresent, omegaPresent, eddyVisPresent
 
        ! Set the cell range to be copied from the buffer.
 
        ibeg = lbound(buffer, 1); iend = ubound(buffer, 1)
        jbeg = lbound(buffer, 2); jend = ubound(buffer, 2)
        kbeg = lbound(buffer, 3); kend = ubound(buffer, 3)
 
        ! Set the cgns real type and abbreviate the solution variable and
        ! the pointer offset to improve readability.
 
        realtypecgns = setcgnsrealtype()
 
        po = iovar(nbklocal, solid)%pointerOffset
        w => iovar(nbklocal, solid)%w
 
        ! Compute the scales for nu, k and omega.
 
        nuscale = muscale / rhoscale
        kscale = velscale**2
        omegascale = kscale / nuscale
 
        ! First check if k is present.
 
        turbkpresent = .false.
 
        nn = bsearchstrings(cgnsturbk, varnames)
 
        if (nn > 0) then
 
            ! K is present. First determine whether or not a type
            ! mismatch occurs. If so, update nTypeMismatch.
 
            if (realtypecgns /= vartypes(nn)) &
                ntypemismatch = ntypemismatch + 1
 
            ! Read k from the restart file and store it in buffer.
 
            call readrestartvariable(varnames(nn))
 
            ! Copy the variables from buffer into w. Multiply by the scale
            ! factor to obtain the correct non-dimensional value and take
            ! the possible pointer offset into account.
 
            do k = kbeg, kend
                kp = k + po
                do j = jbeg, jend
                    jp = j + po
                    do i = ibeg, iend
                        ip = i + po
                        w(ip, jp, kp, itu1) = kscale * buffer(i, j, k)
                    end do
                end do
            end do
 
            ! Set turbKPresent to .true.
 
            turbkpresent = .true.
 
        end if
 
        ! Check if omega is present.
 
        omegapresent = .false.
 
        nn = bsearchstrings(cgnsturbomega, varnames)
 
        if (nn > 0) then
 
            ! Omega is present. First determine whether or not a type
            ! mismatch occurs. If so, update nTypeMismatch.
 
            if (realtypecgns /= vartypes(nn)) &
                ntypemismatch = ntypemismatch + 1
 
            ! Read omega from the restart file and store it in buffer.
 
            call readrestartvariable(varnames(nn))
 
            ! Copy the variables from buffer into w. Multiply by the scale
            ! factor to obtain the correct non-dimensional value and take
            ! the possible pointer offset into account.
 
            do k = kbeg, kend
                kp = k + po
                do j = jbeg, jend
                    jp = j + po
                    do i = ibeg, iend
                        ip = i + po
                        w(ip, jp, kp, itu2) = omegascale * buffer(i, j, k)
                    end do
                end do
            end do
 
            ! Set omegaPresent to .true.
 
            omegapresent = .true.
 
        end if
 
        ! If omega is not present, check if tau is present and
        ! initialize omega accordingly.
 
        if (.not. omegapresent) then
 
            nn = bsearchstrings(cgnsturbtau, varnames)
 
            if (nn > 0) then
 
                ! Tau is present. First determine whether or not a type
                ! mismatch occurs. If so, update nTypeMismatch.
 
                if (realtypecgns /= vartypes(nn)) &
                    ntypemismatch = ntypemismatch + 1
 
                ! Read tau from the restart file and store it in buffer.
 
                call readrestartvariable(varnames(nn))
 
                ! Transform tau to omega and copy the variables from buffer
                ! into w. Multiply by the scale factor to obtain the correct
                ! non-dimensional value and take the possible pointer offset
                ! into account.
 
                do k = kbeg, kend
                    kp = k + po
                    do j = jbeg, jend
                        jp = j + po
                        do i = ibeg, iend
                            ip = i + po
 
                            val = buffer(i, j, k)
                            w(ip, jp, kp, itu2) = omegascale / max(eps, val)
                        end do
                    end do
                end do
 
                ! Set omegaPresent to .true.
 
                omegapresent = .true.
 
            end if
 
        end if
 
        ! Check if both variables were present.
        ! If so go to the check to transform omega to tau.
 
        if (turbkpresent .and. omegapresent) goto 1001
 
        ! K and omega are not both present. It is tried to construct
        ! their values with the information that is present.
 
        ! Try to read the eddy viscosity.
 
        call readturbeddyvis(ntypemismatch, eddyvispresent)
 
        ! The eddy viscosity is either known or still initialized
        ! to the free stream value. In any case determine the
        ! situation we are dealing with and try to initialize k and
        ! omega accordingly.
 
        if (turbkpresent) then
 
            ! K is present. Compute omega using the eddy viscosity.
            ! Assume that the standard k-omega formula is also valid
            ! for the SST-model.
            ! Take the possible pointer offset into account.
 
            do k = kbeg, kend
                kp = k + po
                do j = jbeg, jend
                    jp = j + po
                    do i = ibeg, iend
                        ip = i + po
                        w(ip, jp, kp, itu2) = w(ip, jp, kp, irho) * w(ip, jp, kp, itu1) &
                                              / rev(i, j, k)
                    end do
                end do
            end do
 
            ! Print a warning that omega was not present and has been
            ! constructed. Only processor 0 does this for block 1.
 
            if ((myid == 0) .and. (nbklocal == 1)) then
 
                print "(a)", "#"
                print "(a)", "#                 Warning"
                print "(a)", "# Omega is not present in the restart file."
                if (eddyvispresent) then
                    print "(a)", "# It is initialized using the turbulent &
                         &kinetic energy and eddy viscosity."
                else
                    print "(a)", "# It is initialized using the turbulent &
                         &kinetic energy and free stream eddy &
                         &viscosity."
                end if
                print "(a)", "#"
 
            end if
 
            ! K and omega are initialized.
            ! Go to the check to transform omega to tau.
 
            goto 1001
 
        end if
 
        if (omegapresent) then
 
            ! Omega is present. Compute k using the eddy viscosity.
            ! Assume that the standard k-omega formula is also valid
            ! for the SST-model.
            ! Take the possible pointer offset into account.
 
            do k = kbeg, kend
                kp = k + po
                do j = jbeg, jend
                    jp = j + po
                    do i = ibeg, iend
                        ip = i + po
                        w(ip, jp, kp, itu1) = rev(i, j, k) * w(ip, jp, kp, itu2) &
                                              / w(ip, jp, kp, irho)
                    end do
                end do
            end do
 
            ! Print a warning that k was not present and has been
            ! constructed. Only processor 0 does this for block 1.
 
            if ((myid == 0) .and. (nbklocal == 1)) then
 
                print "(a)", "#"
                print "(a)", "#                 Warning"
                print "(a)", "# Turbulent kinetic energy is not present &
                     &in the restart file."
                if (eddyvispresent) then
                    print "(a)", "# It is initialized using omega and &
                         &the eddy viscosity."
                else
                    print "(a)", "# It is initialized using omega and &
                         &the free stream eddy viscosity."
                end if
                print "(a)", "#"
 
            end if
 
            ! K and omega are initialized.
            ! Go to the check to transform omega to tau.
 
            goto 1001
 
        end if
 
        ! Both k and omega are not present. Use a guess for omega
        ! and compute k using the known value of the eddy viscosity.
        ! As the laminar viscosity is not yet known, set it to the
        ! free-stream value.
 
        rlv = muinf
        call initkomega(po)
 
        ! Print a warning that both k and omega are not present in
        ! the restart file. Only processor 0 does this for block 1.
 
        if ((myid == 0) .and. (nbklocal == 1)) then
 
            print "(a)", "#"
            print "(a)", "#                 Warning"
            print "(a)", "# The turbulent kinetic energy and omega are &
                 &not present in the restart file."
            if (eddyvispresent) then
                print "(a)", "# They have been initialized using the &
                     &eddy viscosity."
            else
                print "(a)", "# The default initialization has been used."
            end if
 
            print "(a)", "#"
 
        end if
 
        ! For the k-tau model omega must be transformed to tau.
        ! Take the possible pointer offset into account.
 
1001    select case (turbmodel)
 
        case (ktau)
 
            do k = kbeg, kend
                kp = k + po
                do j = jbeg, jend
                    jp = j + po
                    do i = ibeg, iend
                        ip = i + po
                        w(ip, jp, kp, itu2) = one / w(ip, jp, kp, itu2)
                    end do
                end do
            end do
 
        end select
 
    end subroutine readturbkwtype
 
    subroutine readturbsa(nTypeMismatch)
        !
        !       readTurbSA reads or constructs the nu tilde transport
        !       variable for the Spalart-Allmaras type turbulence models.
        !       If no information could be retrieved some engineering guess of
        !       the turbulent variables is made.
        !
        use constants
        use cgnsnames
        use communication, only: myid
        use blockpointers, only: w, nbklocal, rlv, rev
        use iomodule, only: iovar
        use utils, only: setcgnsrealtype, terminate
        use sorting, only: bsearchstrings
        use flowvarrefstate, only: muinf, winf
        use turbutils, only: sanuknowneddyratio
        implicit none
        !
        !      Subroutine argument.
        !
        integer(kind=intType), intent(inout) :: nTypeMismatch
        !
        !      Local variables.
        !
        integer :: realTypeCGNS
 
        integer(kind=intType) :: i, j, k, nn, po, ip, jp, kp
        integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, kBeg, kEnd
 
        real(kind=realtype) :: nuscale, ratio, nu
 
        logical :: eddyVisPresent
 
        ! Set the cell range to be copied from the buffer.
 
        ibeg = lbound(buffer, 1); iend = ubound(buffer, 1)
        jbeg = lbound(buffer, 2); jend = ubound(buffer, 2)
        kbeg = lbound(buffer, 3); kend = ubound(buffer, 3)
 
        ! Set the cgns real type and abbreviate the solution variable and
        ! the pointer offset to improve readability.
        ! Also compute the kinematic viscosity scale.
 
        realtypecgns = setcgnsrealtype()
 
        po = iovar(nbklocal, solid)%pointerOffset
        w => iovar(nbklocal, solid)%w
 
        nuscale = muscale / rhoscale
 
        ! Check if the nu tilde variable is present.
 
        nn = bsearchstrings(cgnsturbsanu, varnames)
 
        nutildepresent: if (nn > 0) then
 
            ! Nu tilde is present. First determine whether or not a type
            ! mismatch occurs. If so, update nTypeMismatch.
 
            if (realtypecgns /= vartypes(nn)) &
                ntypemismatch = ntypemismatch + 1
 
            ! Read nu tilde from the restart file and store
            ! it in buffer.
 
            call readrestartvariable(varnames(nn))
 
            ! Copy the variables from buffer into w and take the possible
            ! pointer offset into account.
 
            do k = kbeg, kend
                kp = k + po
                do j = jbeg, jend
                    jp = j + po
                    do i = ibeg, iend
                        ip = i + po
                        w(ip, jp, kp, itu1) = nuscale * buffer(i, j, k)
                    end do
                end do
            end do
 
            ! Variable is read, so a return can be made.
 
            return
 
        end if nutildepresent
 
        ! NuTilde is not present. Try to construct the eddy viscosity.
 
        call readturbeddyvis(ntypemismatch, eddyvispresent)
 
        ! Check if the eddy viscosity has been constructed.
 
        eddypresent: if (eddyvispresent) then
 
            ! Eddy viscosity is present. As the laminar viscosity is not
            ! yet known, set it to the free-stream value.
 
            rlv = muinf
 
            ! Compute nuTilde from the known ratio of eddy and laminar
            ! viscosity. Take the possible pointer offset into account.
 
            do k = kbeg, kend
                kp = k + po
                do j = jbeg, jend
                    jp = j + po
                    do i = ibeg, iend
                        ip = i + po
 
                        ! Compute the eddy viscosity ratio and the laminar
                        ! kinematic viscosity and call the function to
                        ! compute the nu tilde variable.
 
                        ratio = rev(i, j, k) / rlv(i, j, k)
                        nu = rlv(i, j, k) / w(ip, jp, kp, irho)
                        w(ip, jp, kp, itu1) = sanuknowneddyratio(ratio, nu)
 
                    end do
                end do
            end do
 
            ! Print a warning that nu tilde has been constructed and
            ! not read. Only processor 0 does this for block 1.
 
            if ((myid == 0) .and. (nbklocal == 1)) then
 
                print "(a)", "#"
                print "(a)", "#                 Warning"
                print "(a)", "# Nu tilde for Spalart-Allmaras model not &
                     &present in the restart file."
                print "(a)", "# Variable has been reconstructed from &
                     &the eddy viscosity ratio."
                print "(a)", "#"
 
            end if
 
            ! Variable is constructed, so a return can be made.
 
            return
 
        end if eddypresent
 
        ! No turbulence info is present in the restart file.
        ! Initialize nu tilde to the free stream value.
        ! Take the possible pointer offset into account.
 
        do k = kbeg, kend
            kp = k + po
            do j = jbeg, jend
                jp = j + po
                do i = ibeg, iend
                    ip = i + po
                    w(ip, jp, kp, itu1) = winf(itu1)
                end do
            end do
        end do
 
        ! Print a warning that nu tilde has been set to the
        ! free stream values. Only processor 0 does this for block 1.
 
        if ((myid == 0) .and. (nbklocal == 1)) then
 
            print "(a)", "#"
            print "(a)", "#                 Warning"
            print "(a)", "# No turbulent info present in the restart file."
            print "(a)", "# Nu tilde for Spalart-Allmaras model has &
                 &been set to the free stream value."
            print "(a)", "#"
 
        end if
 
    end subroutine readturbsa
    subroutine readturbv2f(nTypeMismatch)
        !
        !       readTurbV2f reads or constructs the four transport variables
        !       for the v2f model. If no information could be retrieved some
        !       engineering guess of the turbulent variables is made.
        !
        use constants
        use cgnsnames
        use communication, only: myid
        use blockpointers, only: w, nbklocal
        use iomodule, only: iovar
        use utils, only: setcgnsrealtype, terminate
        use sorting, only: bsearchstrings
        use flowvarrefstate, only: winf, nt1, nt2
        implicit none
        !
        !      Subroutine argument.
        !
        integer(kind=intType), intent(inout) :: nTypeMismatch
        !
        !      Local variables.
        !
        integer :: realTypeCGNS, itu
        integer, dimension(4) :: indW
 
        integer(kind=intType) :: i, j, k, ii, nn, po, ip, jp, kp
        integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, kBeg, kEnd
 
        real(kind=realtype) :: nuscale, kscale, epsscale, fscale
 
        real(kind=realtype), dimension(4) :: turbscale
 
        character(len=maxCGNSNameLen), dimension(4) :: namesVar
 
        ! Set the cell range to be copied from the buffer.
 
        ibeg = lbound(buffer, 1); iend = ubound(buffer, 1)
        jbeg = lbound(buffer, 2); jend = ubound(buffer, 2)
        kbeg = lbound(buffer, 3); kend = ubound(buffer, 3)
 
        ! Set the cgns real type and abbreviate the solution variable and
        ! the pointer offset to improve readability.
 
        realtypecgns = setcgnsrealtype()
 
        po = iovar(nbklocal, solid)%pointerOffset
        w => iovar(nbklocal, solid)%w
 
        ! Set the names and indices for the four variables.
 
        indw(1) = itu1; namesvar(1) = cgnsturbk
        indw(2) = itu2; namesvar(2) = cgnsturbepsilon
        indw(3) = itu3; namesvar(3) = cgnsturbv2
        indw(4) = itu4; namesvar(4) = cgnsturbf
 
        ! Compute the scales for nu, k, epsilon and f; v2 has the same
        ! scaling as k.
 
        nuscale = muscale / rhoscale
        kscale = velscale**2
        fscale = kscale / nuscale
        epsscale = kscale * fscale
 
        turbscale(1) = kscale
        turbscale(2) = epsscale
        turbscale(3) = kscale
        turbscale(4) = fscale
 
        ! Loop over the four variables of the v2f model.
 
        varloop: do ii = 1, 4
 
            ! Find the index of the variable in the solution file and check
            ! if it is present. If not exit the loop.
 
            nn = bsearchstrings(namesvar(ii), varnames)
 
            if (nn == 0) exit
 
            ! Variable is present. First determine whether or not a type
            ! mismatch occurs. If so, update nTypeMismatch.
 
            if (realtypecgns /= vartypes(nn)) &
                ntypemismatch = ntypemismatch + 1
 
            ! Read the variable from the restart file and store
            ! it in buffer.
 
            call readrestartvariable(varnames(nn))
 
            ! Copy the variables from buffer into w.
            ! Take the possible pointer offset into account.
 
            itu = indw(ii)
 
            do k = kbeg, kend
                kp = k + po
                do j = jbeg, jend
                    jp = j + po
                    do i = ibeg, iend
                        ip = i + po
                        w(ip, jp, kp, itu) = turbscale(ii) * buffer(i, j, k)
                    end do
                end do
            end do
 
        end do varloop
 
        ! Check if all variables were present. If not, set all turbulence
        ! variables to the free-stream values.
 
        testpresent: if (ii <= 4) then
 
            ! Not all variables are present. Set all 4 to the free-stream
            ! values. Take the possible pointer offset into account.
 
            do ii = nt1, nt2
                do k = kbeg, kend
                    kp = k + po
                    do j = jbeg, jend
                        jp = j + po
                        do i = ibeg, iend
                            ip = i + po
                            w(ip, jp, kp, ii) = winf(ii)
                        end do
                    end do
                end do
            end do
 
            ! Print a warning that the turbulence has been initialized to
            ! the free-stream. Only processor 0 does this for block 1.
 
            if ((myid == 0) .and. (nbklocal == 1)) then
 
                print "(a)", "#"
                print "(a)", "#                 Warning"
                print "(a)", "# Not all turbulence variables are present &
                     &for the v2f model."
                print "(a)", "# They have been initialized to the free &
                     &stream values."
                print "(a)", "#"
 
            end if
 
        end if testpresent
 
    end subroutine readturbv2f
 
    subroutine readturbvar(nTypeMismatch)
        !
        !       readTurbvar controls the reading of the turbulent variables
        !       for a restart. It calls the routine, which corresponds to the
        !       turbulence model used.
        !
        use constants
        use communication, only: myid, adflow_comm_world
        use inputphysics, only: equations, turbmodel
        use utils, only: terminate
        implicit none
        !
        !      Subroutine argument.
        !
        integer(kind=intType), intent(inout) :: nTypeMismatch
        !
        !      Local variables.
        !
        integer :: ierr
 
        ! Check if the rans equations must be solved. If not return.
 
        if (equations /= ransequations) return
 
        ! Determine the turbulence model to be used and call the
        ! appropriate subroutine.
 
        select case (turbmodel)
 
        case (spalartallmaras, spalartallmarasedwards)
            call readturbsa(ntypemismatch)
 
            ! !===============================================================
 
            ! case (komegaWilcox, komegaModified, menterSST, ktau)
            !   call readTurbKwType(nTypeMismatch)
 
            ! !===============================================================
 
            ! case (v2f)
            !   call readTurbV2f(nTypeMismatch)
 
            !===============================================================
 
        case default
            if (myid == 0) &
                 call terminate("readTurbvar", "Restart not implemented &
                 &for this turbulence model.")
            call mpi_barrier(adflow_comm_world, ierr)
 
        end select
 
    end subroutine readturbvar
 
    subroutine readxmomentum(nTypeMismatch)
        !
        !       readXmomentum reads the x-momentum variable from the given
        !       place in the cgns file. If the x-momentum itself is not stored
        !       then it is tried to construct it from the x-velocity and
        !       density; it is assumed that the latter is already stored in
        !       the pointer variable w.
        !       If it is not possible to create the x-velocity an error
        !       message is printed and the program will stop.
        !       It is assumed that the pointers in blockPointers already
        !       point to the correct block.
        !
        use constants
        use cgnsnames
        use blockpointers, only: w, nbklocal
        use iomodule, only: iovar
        use utils, only: setcgnsrealtype, terminate
        use sorting, only: bsearchstrings
 
        implicit none
        !
        !      Subroutine argument.
        !
        integer(kind=intType), intent(inout) :: nTypeMismatch
        !
        !      Local variables
        !
        integer :: realTypeCGNS
 
        integer(kind=intType) :: i, j, k, nn, mm, po, ip, jp, kp
        integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, kBeg, kEnd
 
        real(kind=realtype) :: momscale
 
        ! Set the cell range to be copied from the buffer.
 
        ibeg = lbound(buffer, 1); iend = ubound(buffer, 1)
        jbeg = lbound(buffer, 2); jend = ubound(buffer, 2)
        kbeg = lbound(buffer, 3); kend = ubound(buffer, 3)
 
        ! Compute the momentum scaling factor, set the cgns real type and
        ! abbreviate the solution variable and the pointer offset to
        ! improve readability.
 
        momscale = rhoscale * velscale
        realtypecgns = setcgnsrealtype()
 
        po = iovar(nbklocal, solid)%pointerOffset
        w => iovar(nbklocal, solid)%w
 
        ! Find out if the X-momentum is present in the solution file.
 
        nn = bsearchstrings(cgnsmomx, varnames)
 
        testmxpresent: if (nn > 0) then
 
            ! X-momentum is present. First determine whether or not a type
            ! mismatch occurs. If so, update nTypeMismatch.
 
            if (realtypecgns /= vartypes(nn)) &
                ntypemismatch = ntypemismatch + 1
 
            ! Read the x-momentum from the restart file and store it in buffer.
 
            call readrestartvariable(varnames(nn))
 
            ! Copy the variables from buffer into w. Multiply by the scale
            ! factor to obtain the correct non-dimensional value and take
            ! the possible pointer offset into account.
 
            do k = kbeg, kend
                kp = k + po
                do j = jbeg, jend
                    jp = j + po
                    do i = ibeg, iend
                        ip = i + po
                        w(ip, jp, kp, imx) = buffer(i, j, k) * momscale
                    end do
                end do
            end do
 
            ! X-momentum is read, so a return can be made.
 
            return
 
        end if testmxpresent
 
        ! X-momentum is not present. Check for x-velocity.
 
        nn = bsearchstrings(cgnsvelx, varnames)
 
        testvxpresent: if (nn > 0) then
 
            ! X-velocity is present. First determine whether or not a type
            ! mismatch occurs. If so, update nTypeMismatch.
 
            if (realtypecgns /= vartypes(nn)) &
                ntypemismatch = ntypemismatch + 1
 
            ! Read the x-velocity from the restart file and store it in buffer.
 
            call readrestartvariable(varnames(nn))
 
            ! Copy the variables from buffer into w. Multiply by the
            ! density and velocity scaling factor to obtain to correct
            ! non-dimensional value. Take the possible pointer offset
            ! into account.
 
            do k = kbeg, kend
                kp = k + po
                do j = jbeg, jend
                    jp = j + po
                    do i = ibeg, iend
                        ip = i + po
                        w(ip, jp, kp, imx) = buffer(i, j, k) * w(ip, jp, kp, irho) * velscale
                    end do
                end do
            end do
 
            ! X-momentum is constructed, so a return can be made.
 
            return
 
        end if testvxpresent
 
        ! X-momentum could not be created. Terminate.
 
        call terminate("readXmomentum", &
                       "X-Momentum could not be created")
 
    end subroutine readxmomentum
 
    subroutine readxvelocity(nTypeMismatch)
        !
        !       readXvelocity reads the x-velocity variable from the given
        !       place in the cgns file. If the x-velocity itself is not stored
        !       then it is tried to construct it from the x-momentum and
        !       density; it is assumed that the latter is already stored in
        !       the pointer variable w.
        !       If it is not possible to create the x-velocity an error
        !       message is printed and the program will stop.
        !       It is assumed that the pointers in blockPointers already
        !       point to the correct block.
        !
        use constants
        use cgnsnames
        use blockpointers, only: w, nbklocal
        use iomodule, only: iovar
        use utils, only: setcgnsrealtype, terminate
        use sorting, only: bsearchstrings
        !
        !      Subroutine argument.
        !
        integer(kind=intType), intent(inout) :: nTypeMismatch
        !
        !      Local variables
        !
        integer :: realTypeCGNS
 
        integer(kind=intType) :: i, j, k, nn, po, ip, jp, kp
        integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, kBeg, kEnd
 
        real(kind=realtype) :: scale
 
        ! Set the cell range to be copied from the buffer.
 
        ibeg = lbound(buffer, 1); iend = ubound(buffer, 1)
        jbeg = lbound(buffer, 2); jend = ubound(buffer, 2)
        kbeg = lbound(buffer, 3); kend = ubound(buffer, 3)
 
        ! Set the cgns real type and abbreviate the solution variable and
        ! the pointer offset to improve readability.
 
        realtypecgns = setcgnsrealtype()
 
        po = iovar(nbklocal, solid)%pointerOffset
        w => iovar(nbklocal, solid)%w
 
        ! Find out if the x-velocity is present in the solution file.
 
        nn = bsearchstrings(cgnsvelx, varnames)
        if (nn > 0) then
 
            ! X-velocity is present. First determine whether or not a type
            ! mismatch occurs. If so, update nTypeMismatch.
 
            if (realtypecgns /= vartypes(nn)) &
                ntypemismatch = ntypemismatch + 1
 
            ! Read the x-velocity from the restart file and store it
            ! in buffer.
 
            call readrestartvariable(varnames(nn))
 
            ! Copy the variables from buffer into w. Multiply by the scale
            ! factor to obtain the correct nondimensional value and take
            ! the possible pointer offset into account.
 
            do k = kbeg, kend
                kp = k + po
                do j = jbeg, jend
                    jp = j + po
                    do i = ibeg, iend
                        ip = i + po
                        w(ip, jp, kp, ivx) = buffer(i, j, k) * velscale
                    end do
                end do
            end do
 
            ! X-velocity is read, so a return can be made.
 
            return
 
        end if
 
        ! X-velocity not present. Check for x-momentum.
 
        nn = bsearchstrings(cgnsmomx, varnames)
        if (nn > 0) then
 
            ! X-momentum is present. First determine whether or not a type
            ! mismatch occurs. If so, update nTypeMismatch.
 
            if (realtypecgns /= vartypes(nn)) &
                ntypemismatch = ntypemismatch + 1
 
            ! Read the x-momentum from the restart file and store
            ! it in buffer.
 
            call readrestartvariable(varnames(nn))
 
            ! Construct the x-velocity; it is assumed that the density is
            ! already stored in w. Multiply by the momentum scale factor
            ! to obtain the correct nondimensional value and take the
            ! possible pointer offset into account.
 
            scale = rhoscale * velscale
 
            do k = kbeg, kend
                kp = k + po
                do j = jbeg, jend
                    jp = j + po
                    do i = ibeg, iend
                        ip = i + po
                        w(ip, jp, kp, ivx) = buffer(i, j, k) * scale / w(ip, jp, kp, irho)
                    end do
                end do
            end do
 
            ! X-velocity is constructed, so a return can be made.
 
            return
 
        end if
 
        ! Not able to determine the x-velocity.
        ! Print an error message and exit.
 
        call terminate("readXvelocity", &
             "Not able to retrieve x-velocity from the &
             &variables in the restart file.")
 
    end subroutine readxvelocity
 
    subroutine readymomentum(nTypeMismatch)
        !
        !       readYmomentum reads the y-momentum variable from the given
        !       place in the cgns file. If the y-momentum itself is not stored
        !       then it is tried to construct it from the y-velocity and
        !       density; it is assumed that the latter is already stored in
        !       the pointer variable w.
        !       If it is not possible to create the y-velocity an error
        !       message is printed and the program will stop.
        !       It is assumed that the pointers in blockPointers already
        !       point to the correct block.
        !
        use constants
        use cgnsnames
        use blockpointers, only: w, nbklocal
        use iomodule, only: iovar
        use utils, only: setcgnsrealtype, terminate
        use sorting, only: bsearchstrings
        implicit none
        !
        !      Subroutine argument.
        !
        integer(kind=intType), intent(inout) :: nTypeMismatch
        !
        !      Local variables
        !
        integer :: realTypeCGNS
 
        integer(kind=intType) :: i, j, k, nn, po, ip, jp, kp
        integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, kBeg, kEnd
 
        real(kind=realtype) :: momscale
 
        ! Set the cell range to be copied from the buffer.
 
        ibeg = lbound(buffer, 1); iend = ubound(buffer, 1)
        jbeg = lbound(buffer, 2); jend = ubound(buffer, 2)
        kbeg = lbound(buffer, 3); kend = ubound(buffer, 3)
 
        ! Compute the momentum scaling factor, set the cgns real type and
        ! abbreviate the solution variable and the pointer offset to
        ! improve readability.
 
        momscale = rhoscale * velscale
        realtypecgns = setcgnsrealtype()
 
        po = iovar(nbklocal, solid)%pointerOffset
        w => iovar(nbklocal, solid)%w
 
        ! Find out if the Y-momentum is present in the solution file.
 
        nn = bsearchstrings(cgnsmomy, varnames)
 
        testmypresent: if (nn > 0) then
 
            ! Y-momentum is present. First determine whether or not a type
            ! mismatch occurs. If so, update nTypeMismatch.
 
            if (realtypecgns /= vartypes(nn)) &
                ntypemismatch = ntypemismatch + 1
 
            ! Read the y-momentum from the restart file and store it in buffer.
 
            call readrestartvariable(varnames(nn))
 
            ! Copy the variables from buffer into w. Multiply by the scale
            ! factor to obtain the correct non-dimensional value and take
            ! the possible pointer offset into account.
 
            do k = kbeg, kend
                kp = k + po
                do j = jbeg, jend
                    jp = j + po
                    do i = ibeg, iend
                        ip = i + po
                        w(ip, jp, kp, imy) = buffer(i, j, k) * momscale
                    end do
                end do
            end do
 
            ! Y-momentum is read, so a return can be made.
 
            return
 
        end if testmypresent
 
        ! Y-momentum is not present. Check for y-velocity.
 
        nn = bsearchstrings(cgnsvely, varnames)
 
        testvypresent: if (nn > 0) then
 
            ! Y-velocity is present. First determine whether or not a type
            ! mismatch occurs. If so, update nTypeMismatch.
 
            if (realtypecgns /= vartypes(nn)) &
                ntypemismatch = ntypemismatch + 1
 
            ! Read the y-velocity from the restart file and store it in buffer.
 
            call readrestartvariable(varnames(nn))
 
            ! Copy the variables from buffer into w. Multiply by the
            ! density and velocity scaling factor to obtain to correct
            ! non-dimensional value. Take the possible pointer offset
            ! into account.
 
            do k = kbeg, kend
                kp = k + po
                do j = jbeg, jend
                    jp = j + po
                    do i = ibeg, iend
                        ip = i + po
                        w(ip, jp, kp, imy) = buffer(i, j, k) * w(ip, jp, kp, irho) * velscale
                    end do
                end do
            end do
 
            ! Y-momentum is constructed, so a return can be made.
 
            return
 
        end if testvypresent
 
        ! Y-momentum could not be created. Terminate.
 
        call terminate("readYmomentum", &
                       "Y-Momentum could not be created")
 
    end subroutine readymomentum
 
    subroutine readyvelocity(nTypeMismatch)
        !
        !       readYvelocity reads the y-velocity variable from the given
        !       place in the cgns file. If the y-velocity itself is not stored
        !       then it is tried to construct it from the y-momentum and
        !       density; it is assumed that the latter is already stored in
        !       the pointer variable w.
        !       If it is not possible to create the y-velocity an error
        !       message is printed and the program will stop.
        !       It is assumed that the pointers in blockPointers already
        !       point to the correct block.
        !
        use constants
        use cgnsnames
        use blockpointers, only: w, nbklocal
        use iomodule, only: iovar
        use utils, only: setcgnsrealtype, terminate
        use sorting, only: bsearchstrings
 
        implicit none
        !
        !      Subroutine argument.
        !
        integer(kind=intType), intent(inout) :: nTypeMismatch
        !
        !      Local variables
        !
        integer :: realTypeCGNS
 
        integer(kind=intType) :: i, j, k, nn, po, ip, jp, kp
        integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, kBeg, kEnd
 
        real(kind=realtype) :: scale
 
        ! Set the cell range to be copied from the buffer.
 
        ibeg = lbound(buffer, 1); iend = ubound(buffer, 1)
        jbeg = lbound(buffer, 2); jend = ubound(buffer, 2)
        kbeg = lbound(buffer, 3); kend = ubound(buffer, 3)
 
        ! Set the cgns real type and abbreviate the solution variable and
        ! the pointer offset to improve readability.
 
        realtypecgns = setcgnsrealtype()
 
        po = iovar(nbklocal, solid)%pointerOffset
        w => iovar(nbklocal, solid)%w
 
        ! Find out if the y-velocity is present in the solution file.
 
        nn = bsearchstrings(cgnsvely, varnames)
        if (nn > 0) then
 
            ! Y-velocity is present. First determine whether or not a type
            ! mismatch occurs. If so, update nTypeMismatch.
 
            if (realtypecgns /= vartypes(nn)) &
                ntypemismatch = ntypemismatch + 1
 
            ! Read the y-velocity from the restart file and store it
            ! in buffer.
 
            call readrestartvariable(varnames(nn))
 
            ! Copy the variables from buffer into w. Multiply by the scale
            ! factor to obtain the correct nondimensional value and take
            ! the possible pointer offset into account.
 
            do k = kbeg, kend
                kp = k + po
                do j = jbeg, jend
                    jp = j + po
                    do i = ibeg, iend
                        ip = i + po
                        w(ip, jp, kp, ivy) = buffer(i, j, k) * velscale
                    end do
                end do
            end do
 
            ! Y-velocity is read, so a return can be made.
 
            return
 
        end if
 
        ! Y-velocity not present. Check for y-momentum.
 
        nn = bsearchstrings(cgnsmomy, varnames)
        if (nn > 0) then
 
            ! Y-momentum is present. First determine whether or not a type
            ! mismatch occurs. If so, update nTypeMismatch.
 
            if (realtypecgns /= vartypes(nn)) &
                ntypemismatch = ntypemismatch + 1
 
            ! Read the y-momentum from the restart file and store
            ! it in buffer.
 
            call readrestartvariable(varnames(nn))
 
            ! Construct the y-velocity; it is assumed that the density is
            ! already stored in w. Multiply by the momentum scale factor
            ! to obtain the correct nondimensional value and take the
            ! possible pointer offset into account.
 
            scale = rhoscale * velscale
 
            do k = kbeg, kend
                kp = k + po
                do j = jbeg, jend
                    jp = j + po
                    do i = ibeg, iend
                        ip = i + po
                        w(ip, jp, kp, ivy) = buffer(i, j, k) * scale / w(ip, jp, kp, irho)
                    end do
                end do
            end do
 
            ! Y-velocity is constructed, so a return can be made.
 
            return
 
        end if
 
        ! Not able to determine the y-velocity.
        ! Print an error message and exit.
 
        call terminate("readYvelocity", &
             "Not able to retrieve y-velocity from the &
             &variables in the restart file.")
 
    end subroutine readyvelocity
    subroutine readzmomentum(nTypeMismatch)
        !
        !       readZmomentum reads the z-momentum variable from the given
        !       place in the cgns file. If the z-momentum itself is not stored
        !       then it is tried to construct it from the z-velocity and
        !       density; it is assumed that the latter is already stored in
        !       the pointer variable w.
        !       If it is not possible to create the z-velocity an error
        !       message is printed and the program will stop.
        !       It is assumed that the pointers in blockPointers already
        !       point to the correct block.
        !
        use constants
        use cgnsnames
        use blockpointers, only: w, nbklocal
        use iomodule, only: iovar
        use utils, only: setcgnsrealtype, terminate
        use sorting, only: bsearchstrings
        implicit none
        !
        !      Subroutine argument.
        !
        integer(kind=intType), intent(inout) :: nTypeMismatch
        !
        !      Local variables
        !
        integer :: realTypeCGNS
 
        integer(kind=intType) :: i, j, k, nn, po, ip, jp, kp
        integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, kBeg, kEnd
 
        real(kind=realtype) :: momscale
 
        ! Set the cell range to be copied from the buffer.
 
        ibeg = lbound(buffer, 1); iend = ubound(buffer, 1)
        jbeg = lbound(buffer, 2); jend = ubound(buffer, 2)
        kbeg = lbound(buffer, 3); kend = ubound(buffer, 3)
 
        ! Compute the momentum scaling factor, set the cgns real type and
        ! abbreviate the solution variable and the pointer offset to
        ! improve readability.
 
        momscale = rhoscale * velscale
        realtypecgns = setcgnsrealtype()
 
        po = iovar(nbklocal, solid)%pointerOffset
        w => iovar(nbklocal, solid)%w
 
        ! Find out if the Z-momentum is present in the solution file.
 
        nn = bsearchstrings(cgnsmomz, varnames)
 
        testmzpresent: if (nn > 0) then
 
            ! Z-momentum is present. First determine whether or not a type
            ! mismatch occurs. If so, update nTypeMismatch.
 
            if (realtypecgns /= vartypes(nn)) &
                ntypemismatch = ntypemismatch + 1
 
            ! Read the z-momentum from the restart file and store it in buffer.
 
            call readrestartvariable(varnames(nn))
 
            ! Copy the variables from buffer into w. Multiply by the scale
            ! factor to obtain the correct non-dimensional value and take
            ! the possible pointer offset into account.
 
            do k = kbeg, kend
                kp = k + po
                do j = jbeg, jend
                    jp = j + po
                    do i = ibeg, iend
                        ip = i + po
                        w(ip, jp, kp, imz) = buffer(i, j, k) * momscale
                    end do
                end do
            end do
 
            ! Z-momentum is read, so a return can be made.
 
            return
 
        end if testmzpresent
 
        ! Z-momentum is not present. Check for z-velocity.
 
        nn = bsearchstrings(cgnsvelz, varnames)
 
        testvzpresent: if (nn > 0) then
 
            ! Z-velocity is present. First determine whether or not a type
            ! mismatch occurs. If so, update nTypeMismatch.
 
            if (realtypecgns /= vartypes(nn)) &
                ntypemismatch = ntypemismatch + 1
 
            ! Read the z-velocity from the restart file and store it in buffer.
 
            call readrestartvariable(varnames(nn))
 
            ! Copy the variables from buffer into w. Multiply by the
            ! density and velocity scaling factor to obtain to correct
            ! non-dimensional value. Take the possible pointer offset
            ! into account.
 
            do k = kbeg, kend
                kp = k + po
                do j = jbeg, jend
                    jp = j + po
                    do i = ibeg, iend
                        ip = i + po
                        w(ip, jp, kp, imz) = buffer(i, j, k) * w(ip, jp, kp, irho) * velscale
                    end do
                end do
            end do
 
            ! Z-momentum is constructed, so a return can be made.
 
            return
 
        end if testvzpresent
 
        ! Z-momentum could not be created. Terminate.
 
        call terminate("readZmomentum", &
                       "Z-Momentum could not be created")
 
    end subroutine readzmomentum
    subroutine readzvelocity(nTypeMismatch)
        !
        !       readZvelocity reads the z-velocity variable from the given
        !       place in the cgns file. If the z-velocity itself is not stored
        !       then it is tried to construct it from the z-momentum and
        !       density; it is assumed that the latter is already stored in
        !       the pointer variable w.
        !       If it is not possible to create the z-velocity an error
        !       message is printed and the program will stop.
        !       It is assumed that the pointers in blockPointers already
        !       point to the correct block.
        !
        use constants
        use cgnsnames
        use blockpointers, only: w, nbklocal
        use iomodule, only: iovar
        use utils, only: setcgnsrealtype, terminate
        use sorting, only: bsearchstrings
        implicit none
        !
        !      Subroutine argument.
        !
        integer(kind=intType), intent(inout) :: nTypeMismatch
        !
        !      Local variables
        !
        integer :: realTypeCGNS
 
        integer(kind=intType) :: i, j, k, nn, po, ip, jp, kp
        integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, kBeg, kEnd
 
        real(kind=realtype) :: scale
 
        ! Set the cell range to be copied from the buffer.
 
        ibeg = lbound(buffer, 1); iend = ubound(buffer, 1)
        jbeg = lbound(buffer, 2); jend = ubound(buffer, 2)
        kbeg = lbound(buffer, 3); kend = ubound(buffer, 3)
 
        ! Set the cgns real type and abbreviate the solution variable and
        ! the pointer offset to improve readability.
 
        realtypecgns = setcgnsrealtype()
 
        po = iovar(nbklocal, solid)%pointerOffset
        w => iovar(nbklocal, solid)%w
 
        ! Find out if the z-velocity is present in the solution file.
 
        nn = bsearchstrings(cgnsvelz, varnames)
        if (nn > 0) then
 
            ! Z-velocity is present. First determine whether or not a type
            ! mismatch occurs. If so, update nTypeMismatch.
 
            if (realtypecgns /= vartypes(nn)) &
                ntypemismatch = ntypemismatch + 1
 
            ! Read the z-velocity from the restart file and store it
            ! in buffer.
 
            call readrestartvariable(varnames(nn))
 
            ! Copy the variables from buffer into w. Multiply by the scale
            ! factor to obtain the correct nondimensional value and take
            ! the possible pointer offset into account.
 
            do k = kbeg, kend
                kp = k + po
                do j = jbeg, jend
                    jp = j + po
                    do i = ibeg, iend
                        ip = i + po
                        w(ip, jp, kp, ivz) = buffer(i, j, k) * velscale
                    end do
                end do
            end do
 
            ! Z-velocity is read, so a return can be made.
 
            return
 
        end if
 
        ! Z-velocity not present. Check for z-momentum.
 
        nn = bsearchstrings(cgnsmomz, varnames)
        if (nn > 0) then
 
            ! Z-momentum is present. First determine whether or not a type
            ! mismatch occurs. If so, update nTypeMismatch.
 
            if (realtypecgns /= vartypes(nn)) &
                ntypemismatch = ntypemismatch + 1
 
            ! Read the z-momentum from the restart file and store
            ! it in buffer.
 
            call readrestartvariable(varnames(nn))
 
            ! Construct the z-velocity; it is assumed that the density is
            ! already stored in w. Multiply by the momentum scale factor
            ! to obtain the correct nondimensional value and take the
            ! possible pointer offset into account.
 
            scale = rhoscale * velscale
 
            do k = kbeg, kend
                kp = k + po
                do j = jbeg, jend
                    jp = j + po
                    do i = ibeg, iend
                        ip = i + po
                        w(ip, jp, kp, ivz) = buffer(i, j, k) * scale / w(ip, jp, kp, irho)
                    end do
                end do
            end do
 
            ! Z-velocity is constructed, so a return can be made.
 
            return
 
        end if
 
        ! Not able to determine the z-velocity.
        ! Print an error message and exit.
 
        call terminate("readZvelocity", &
             "Not able to retrieve z-velocity from the &
             &variables in the restart file.")
 
    end subroutine readzvelocity
 
    subroutine readtimehistory(fileIDs)
        !
        !       readTimeHistory attempts to read the time history of an
        !       unsteady computation from the given cgns restart file.
        !       If present it will be stored in the arrays timeArray and
        !       timeDataArray, for which memory is allocated.
        !
        use constants
        use cgnsnames
        use su_cgns
        use inputunsteady, only: ntimestepsfine
        use monitor, only: timedataarray, nmon, ntimestepsrestart, &
                           timeunsteadyrestart, monnames, timearray
        use sorting, only: qsortstrings, bsearchstrings
        use utils, only: setcgnsrealtype, terminate, alloctimearrays
        use commonformats, only: strings
        implicit none
        !
        !      Subroutine arguments.
        !
        integer, dimension(nSolsRead), intent(in) :: fileIDs
 
        !
        !      Local variables.
        !
        integer :: ierr, realTypeCGNS, dummyInt
        integer :: i, nConv, nDim
        integer(kind=cgsize_t) :: nSize(1)
 
        integer(kind=intType) :: j, ii, nn
 
        integer(kind=intType), dimension(:), allocatable :: ind
 
        character(len=maxCGNSNameLen) :: cgnsName
        character(len=maxCGNSNameLen), allocatable, dimension(:) :: &
            convNames, tmpNames
 
        logical :: allConvInfo
 
        ! Store the file ID and the base a bit easier. Note that the time
        ! history only needs to be present in the first solution file.
 
        cgnsind = fileids(1)
        cgnsbase = 1
 
        ! Set the cgns real type.
 
        realtypecgns = setcgnsrealtype()
 
        ! Check if the time history is present by trying to read it.
 
        call cg_biter_read_f(cgnsind, cgnsbase, cgnsname, &
                             dummyint, ierr)
        if (ierr /= all_ok) then
 
            ! No time history present. Set nTimeStepsRestart and
            ! timeUnsteadyRestart to zero, allocate the memory for the
            ! time history of the monitoring variables, print a warning
            ! and return.
 
            ntimestepsrestart = 0
            timeunsteadyrestart = zero
 
            call alloctimearrays(ntimestepsfine)
 
            print "(a)", "#"
            print "(a)", "#                 Warning"
            print "(a)", "# No time history found in restart file."
            print "(a)", "# Starting at timestep 1 on the finest level."
            print "(a)", "#"
 
            return
 
        end if
 
        ! Store the number of old time levels.
 
        ntimestepsrestart = dummyint
 
        ! Go to the place in the cgns file where the time history
        ! should be stored.
 
        call cg_goto_f(cgnsind, cgnsbase, ierr, &
                       "BaseIterativeData_t", 1, "end")
        if (ierr /= all_ok) &
            call terminate("readTimeHistory", &
                           "Something wrong when calling cg_goto_f")
 
        ! Find out how many convergence variables are stored.
 
        call cg_narrays_f(nconv, ierr)
        if (ierr /= all_ok) &
            call terminate("readTimeHistory", &
                           "Something wrong when calling cg_narrays_f")
 
        ! Allocate the memory for convNames, tmpNames and ind.
 
        allocate (convnames(nconv), tmpnames(nconv), ind(nconv), &
                  stat=ierr)
        if (ierr /= 0) &
            call terminate("readTimeHistory", &
                           "Memory allocation failure for convNames, etc.")
 
        ! Read the names of the convergence variables. Store them in
        ! tmpNames as well. Furthermore check the dimension of the
        ! data stored.
 
        do i = 1, nconv
            call cg_array_info_f(i, convnames(i), dummyint, ndim, &
                                 nsize, ierr)
            if (ierr /= all_ok) &
                call terminate("readConvHistory", &
                               "Something wrong when calling cg_array_info_f")
 
            if (ndim /= 1) then
                print "(a)", "#"
                print "(a)", "#                 Warning"
                print strings, "# Dimension of time history for ", trim(convnames(i)), " is not 1."
                print "(a)", "# Information is ignored."
                print "(a)", "#"
 
                ! Screw up the string such that it does not correspond to
                ! a legal name. It is appended, because it is important that
                ! all strings differ.
 
                convnames(i) = convnames(i)//"#$@&^!#$%!"
            end if
 
            if (nsize(1) /= ntimestepsrestart) then
                print "(a)", "#"
                print "(a)", "#                 Warning"
                print strings, "# Inconsistent time history for ", trim(convnames(i)), "."
                print "(a)", "# Displayed information might be incorrect."
                print "(a)", "#"
            end if
 
            ! Copy the name in tmpNames for the sorting.
 
            tmpnames(i) = convnames(i)
        end do
 
        ! Sort convNames in increasing order.
 
        nn = nconv
        call qsortstrings(convnames, nn)
 
        ! Find the numbers for the just sorted convergence names.
 
        do i = 1, nconv
            ii = bsearchstrings(tmpnames(i), convnames)
            ind(ii) = i
        end do
 
        ! Find out whether the old time values are present.
        ! If not the time history stored will be ignored.
 
        ii = bsearchstrings(cgnstimevalue, convnames)
        if (ii == 0) then
            print "(a)", "#"
            print "(a)", "#                 Warning"
            print "(a)", "# No time values found in the time history &
                 &in the restart file."
            print "(a)", "# The rest of the time history is ignored."
            print "(a)", "# Starting at timestep 1 on the finest level."
            print "(a)", "#"
 
            ! Set nTimeStepsRestart and timeUnsteadyRestart to 0.
 
            ntimestepsrestart = 0
            timeunsteadyrestart = zero
        end if
 
        ! Determine the total number of time levels and allocate the
        ! memory for the time history arrays.
 
        j = ntimestepsrestart + ntimestepsfine
        call alloctimearrays(j)
 
        ! If the time values are not present, jump to the place where
        ! the memory of the variables used in this routine is released.
 
        if (ii == 0) goto 99
 
        ! Read the time values.
 
        i = ind(ii)
        call cg_array_read_as_f(i, realtypecgns, timearray, ierr)
        if (ierr /= all_ok) &
             call terminate("readTimeHistory", &
             "Something wrong when calling &
             &cg_array_read_as_f")
 
        ! Set the value of timeUnsteadyRestart to the last value in
        ! timeArray.
 
        timeunsteadyrestart = timearray(ntimestepsrestart)
 
        ! Initialize allConvInfo to .true. and perform the loop over
        ! the number of monitoring variables.
 
        allconvinfo = .true.
 
        do j = 1, nmon
 
            ! Search for the monitoring name in the sorted
            ! convergence names present in the restart file.
 
            ii = bsearchstrings(monnames(j), convnames)
 
            ! Check if the name was found.
 
            if (ii == 0) then
 
                ! Name not present in the restart file. Set allConvInfo
                ! to .false. and the corresponding entries in timeDataArray
                ! to zero
 
                allconvinfo = .false.
                do i = 1, ntimestepsrestart
                    timedataarray(i, j) = zero
                end do
 
            else
 
                ! Name is present in the restart file. Read the corresponding
                ! time history.
 
                i = ind(ii)
                call cg_array_read_as_f(i, realtypecgns, &
                                        timedataarray(1, j), ierr)
                if (ierr /= all_ok) &
                     call terminate("readTimeHistory", &
                     "Something wrong when calling &
                     &cg_array_read_as_f")
            end if
 
        end do
 
        ! Print a warning in case not all the time history could
        ! be retrieved from the restart file.
 
        if (.not. allconvinfo) then
 
            print "(a)", "#"
            print "(a)", "#                 Warning"
            print "(a)", "# Not all the time history could be &
                 &retrieved from the restart file."
            print "(a)", "# Missing information is initialized to zero."
            print "(a)", "#"
 
        end if
 
99      continue
 
        ! Release the memory of the variables allocated in this routine.
 
        deallocate (convnames, tmpnames, ind, stat=ierr)
        if (ierr /= 0) &
            call terminate("readTimeHistory", &
                           "Deallocation error for convNames, etc.")
 
    end subroutine readtimehistory
 
    subroutine scalefactors(fileIDs)
        !
        !       scaleFactors determines the scale factors for the density,
        !       pressure and velocity from either the reference state in the
        !       given cgns file or they are simply set to 1.0; the latter
        !       occurs if the input parameter checkRestartSol is .false.
        !       If no reference state is present checkRestartSol is .true.
        !       An error message will be printed and the program terminates.
        !
        use constants
        use cgnsnames
        use su_cgns
        use communication, only: myid, adflow_comm_world
        use flowvarrefstate, only: pref, muref, rhoref
        use inputio, only: checkrestartsol
        use sorting, only: bsearchstrings, qsortstrings
        use utils, only: setcgnsrealtype, terminate
        implicit none
        !
        !      Subroutine arguments.
        !
        integer, dimension(nSolsRead), intent(in) :: fileIDs
 
        !      Local variables.
        !
        integer :: ierr, realTypeCGNS, typeCGNS
        integer :: i, nDim, nRef
        integer :: nsize
        integer(kind=cgsize_t) :: nsize2(1)
        integer(kind=intType) :: ii, nn
 
        integer(kind=intType), dimension(:), allocatable :: ind
 
        real(kind=cgnsrealtype) :: tmpscale
 
        character(len=maxCGNSNameLen), allocatable, dimension(:) :: &
            refNames, tmpNames
 
        logical :: muScalePresent
 
        ! Store the file ID and the base a bit easier. Note that the
        ! reference state only needs to be present in the first file.
 
        cgnsind = fileids(1)
        cgnsbase = 1
 
        ! Set the cgns real type.
 
        realtypecgns = setcgnsrealtype()
 
        ! Initialize the scale factors to 1.0, i.e. assume that the
        ! correct non-dimensional solution is stored in the restart file.
 
        rhoscale = one
        velscale = one
        pscale = one
        muscale = one
 
        ! Go to the place in the cgns file where the reference state
        ! should be stored.
 
        call cg_goto_f(cgnsind, cgnsbase, ierr, "end")
        if (ierr /= all_ok) &
            call terminate("scaleFactors", &
                           "Something wrong when calling cg_goto_f")
 
        ! Try going to the reference state node. If we get an error code,
        ! it doesn't exist.
 
        call cg_goto_f(cgnsind, cgnsbase, ierr, &
                       "ReferenceState_t", 1, "end")
        if (ierr /= all_ok) then
 
            ! Reference state does not exist. Check if the restart solution
            ! must be checked. If not, return; otherwise print an error
            ! message and terminate the execution. This error message is
            ! only printed by processor 0 to avoid a messy output.
 
            if (.not. checkrestartsol) return
 
            if (myid == 0) &
                 call terminate("scaleFactors", &
                 "Reference state not presented in restart &
                 &file. Scaling factors cannot be determined.")
 
            ! The other processors will wait until they are killed.
 
            call mpi_barrier(adflow_comm_world, ierr)
 
        end if
 
        ! Go to the reference state node.
 
        call cg_goto_f(cgnsind, cgnsbase, ierr, &
                       "ReferenceState_t", 1, "end")
        if (ierr /= all_ok) &
            call terminate("scaleFactors", &
                           "Something wrong when calling cg_goto_f")
 
        ! Found out how many reference variables are stored.
 
        call cg_narrays_f(nref, ierr)
        if (ierr /= all_ok) &
            call terminate("scaleFactors", &
                           "Something wrong when calling cg_narrays_f")
 
        ! Allocate the memory for refNames, tmpNames and ind.
 
        allocate (refnames(nref), tmpnames(nref), ind(nref), &
                  stat=ierr)
        if (ierr /= 0) &
            call terminate("scaleFactors", &
                           "Memory allocation failure for refNames, etc.")
 
        ! Read the names of the reference variables. Store them in
        ! tmpNames as well.
 
        do i = 1, nref
            call cg_array_info_f(i, refnames(i), typecgns, ndim, &
                                 nsize2, ierr)
            if (ierr /= all_ok) &
                call terminate("scaleFactors", &
                               "Something wrong when calling cg_array_info_f")
 
            ! Check the dimension and the size of the array.
            ! Both should be 1. If not, screw up the name such that it
            ! will never be found in the search later on.
 
            if (ndim /= 1 .or. nsize2(1) /= 1) &
                refnames(i) = refnames(i)//"#$@&^!#$%!"
 
            ! And copy it in tmpNames.
 
            tmpnames(i) = refnames(i)
        end do
 
        ! Sort refNames in increasing order.
 
        nn = nref
        call qsortstrings(refnames, nn)
 
        ! Find the numbers for the just sorted reference names.
 
        do i = 1, nref
            ii = bsearchstrings(tmpnames(i), refnames)
            ind(ii) = i
        end do
 
        ! Determine the scale factors if these must be determined.
 
        if (checkrestartsol) then
 
            ! Read the reference density from the restart file.
 
            ii = bsearchstrings(cgnsdensity, refnames)
            if (ii == 0) then
                if (myid == 0) &
                    call terminate("scaleFactors", &
                                   "No reference density found in restart file")
 
                ! The other processors will wait until they are killed.
 
                call mpi_barrier(adflow_comm_world, ierr)
            end if
 
            i = ind(ii)
            call cg_array_read_as_f(i, realtypecgns, tmpscale, ierr)
            if (ierr /= all_ok) &
                 call terminate("scaleFactors", &
                 "Something wrong when calling &
                 &cg_array_read_as_f")
            rhoscale = tmpscale
 
            ! Read the reference pressure from the restart file.
 
            ii = bsearchstrings(cgnspressure, refnames)
            if (ii == 0) then
                if (myid == 0) &
                     call terminate("scaleFactors", &
                     "No reference pressure found in &
                     &restart file")
 
                ! The other processors will wait until they are killed.
 
                call mpi_barrier(adflow_comm_world, ierr)
            end if
 
            i = ind(ii)
            call cg_array_read_as_f(i, realtypecgns, tmpscale, ierr)
            if (ierr /= all_ok) &
                 call terminate("scaleFactors", &
                 "Something wrong when calling &
                 &cg_array_read_as_f")
            pscale = tmpscale
 
            ! Read the reference velocity from the restart file.
 
            ii = bsearchstrings(cgnsvelocity, refnames)
            if (ii == 0) then
                if (myid == 0) &
                     call terminate("scaleFactors", &
                     "No reference velocity found in &
                     &restart file")
 
                ! The other processors will wait until they are killed.
 
                call mpi_barrier(adflow_comm_world, ierr)
            end if
 
            i = ind(ii)
            call cg_array_read_as_f(i, realtypecgns, tmpscale, ierr)
            if (ierr /= all_ok) &
                 call terminate("scaleFactors", &
                 "Something wrong when calling &
                 &cg_array_read_as_f")
            velscale = tmpscale
 
            ! Set muScalePresent to .true. to indicate that it is present
            ! and read or construct the reference molecular viscosity.
 
            muscalepresent = .true.
 
            ii = bsearchstrings(cgnsviscmol, refnames)
            if (ii > 0) then
 
                ! Scale is present; read the value.
 
                i = ind(ii)
                call cg_array_read_as_f(i, realtypecgns, tmpscale, ierr)
                if (ierr /= all_ok) &
                     call terminate("scaleFactors", &
                     "Something wrong when calling &
                     &cg_array_read_as_f")
                muscale = tmpscale
 
            else
 
                ! Try to read the kinematic viscosity.
 
                ii = bsearchstrings(cgnsvisckin, refnames)
                if (ii > 0) then
 
                    ! Scale is present; read the value and multiply it by the
                    ! density.
 
                    i = ind(ii)
                    call cg_array_read_as_f(i, realtypecgns, tmpscale, ierr)
                    if (ierr /= all_ok) &
                         call terminate("scaleFactors", &
                         "Something wrong when calling &
                         &cg_array_read_as_f")
 
                    muscale = tmpscale * rhoscale
 
                else
 
                    ! Final possibility. Try to read the length scale.
 
                    ii = bsearchstrings(cgnslength, refnames)
                    if (ii > 0) then
 
                        ! Scale is present; read the value and create the
                        ! reference viscosity.
 
                        i = ind(ii)
                        call cg_array_read_as_f(i, realtypecgns, tmpscale, ierr)
                        if (ierr /= all_ok) &
                             call terminate("scaleFactors", &
                             "Something wrong when calling &
                             &cg_array_read_as_f")
 
                        muscale = tmpscale * sqrt(pscale * rhoscale)
 
                    else
 
                        ! Set the logical muScalePresent to .false.
 
                        muscalepresent = .false.
 
                    end if
                end if
            end if
 
            ! Create the correct scaling factors for density, pressure,
            ! velocity and possibly viscosity.
 
            rhoscale = rhoscale / rhoref
            pscale = pscale / pref
            velscale = velscale / sqrt(pref / rhoref)
 
            if (muscalepresent) muscale = muscale / muref
 
        end if
 
        ! Release the memory of refNames, tmpNames and ind.
 
        deallocate (refnames, tmpnames, ind, stat=ierr)
        if (ierr /= 0) &
            call terminate("scaleFactors", &
                           "Deallocation error for convNames, etc.")
 
    end subroutine scalefactors
 
    subroutine readrestartvariable(cgnsVarName)
 
        !      readRestartVariable reads the given variable name from the
        !      cgns restart file.
        !
        use constants
        use blockpointers, only: il, jl, kl
        use utils, only: terminate, setcgnsrealtype
        use su_cgns
        implicit none
        !
        !      Subroutine arguments.
        !
        character(len=*), intent(in) :: cgnsVarName
        !
        !      Local variables.
        !
        integer :: ierr, realTypeCGNS
 
        integer(kind=intType) :: i, j, k
 
        ! Set the cgns real type.
 
        realtypecgns = setcgnsrealtype()
 
        ! Check where the solution variables are stored.
 
        locationtest: if (location == cellcenter) then
 
            ! Cell centered values. Read the values directly in the buffer.
 
            call cg_field_read_f(cgnsind, cgnsbase, cgnszone, cgnssol, &
                                 cgnsvarname, realtypecgns, rangemin, &
                                 rangemax, buffer, ierr)
            if (ierr /= all_ok) &
                call terminate("readRestartVariable", &
                               "Something wrong when calling cg_field_read_f")
            else locationtest
 
            ! Vertex centered values. First read the solution in the
            ! array bufferVertex.
 
            call cg_field_read_f(cgnsind, cgnsbase, cgnszone, cgnssol, &
                                 cgnsvarname, realtypecgns, rangemin, &
                                 rangemax, buffervertex, ierr)
            if (ierr /= all_ok) &
                call terminate("readRestartVariable", &
                               "Something wrong when calling cg_field_read_f")
 
            ! Create the cell centered values by averaging the vertex values.
 
            do k = 2, kl
                do j = 2, jl
                    do i = 2, il
                        buffer(i, j, k) = eighth * (buffervertex(i - 1, j - 1, k - 1) &
                                                    + buffervertex(i, j - 1, k - 1) &
                                                    + buffervertex(i - 1, j, k - 1) &
                                                    + buffervertex(i, j, k - 1) &
                                                    + buffervertex(i - 1, j - 1, k) &
                                                    + buffervertex(i, j - 1, k) &
                                                    + buffervertex(i - 1, j, k) &
                                                    + buffervertex(i, j, k))
                    end do
                end do
            end do
 
        end if locationtest
    end subroutine readrestartvariable
 
end module variablereading