inputParamRoutines.F90

Go to the documentation of this file.

module inputparamroutines
 
    use constants, only: maxstringlen, inttype
 
    ! Set the parameter none, which is used as a check to see
    ! whether or not some key parameters were specified.
 
    integer(kind=intType), parameter :: none = 0
 
    ! monDturb:             Whether or not the turbulent residuals
    !                        must be monitored. This must be done via
    !                        this construction, because during the
    !                        reading of the monitoring variables the
    !                        turbulence model might not be known.
    ! monitorSpecified:     Whether or not the monitoring variables
    !                        were specified.
    ! surfaceOutSpecified: Whether or not the surface output
    !                        variables were specified.
    ! volumeOutSpecified:  Whether or not the volume output
    !                        variables were specified.
    ! isoOutSpecified:      Wheter or not the isosurface output
    !                        variables were specified
    logical :: mondturb
    logical :: monitorspecified
    logical :: surfaceoutspecified
    logical :: volumeoutspecified
    logical :: isooutspecified
 
    ! liftDirSpecified: Whether or not the lift direction was
    !                     specified.
 
    logical :: liftdirspecified
 
contains
 
    subroutine checkmonitor
        !
        !       checkMonitor checks and possibly corrects the variables
        !       to be monitored during the convergence.  This depends on the
        !       governing equations to be solved. After the correction the
        !       sequence of the monitoring variable names is changed, such
        !       that the output is independent of the specified sequence.
        !       Furthermore memory is allocated for the arrays used to compute
        !       the monitoring variables and it is checked whether or not the
        !       maximum Mach number of total enthalpy difference is to be
        !       monitored.
        !
        use constants
        use cgnsnames
        use monitor, only: monnames, monglob, monloc, nmonmax, nmonsum, &
                           nmon, monmachorhmax, monref
        use inputphysics, only: equations, flowtype, turbmodel, equationmode
        use inputunsteady, only: timeintegrationscheme
        use sorting, only: qsortintegers, bsearchintegers
        use utils, only: terminate
        implicit none
        !
        !      Local variables.
        !
        integer :: ierr
 
        integer(kind=intType) :: i, ii, nn
        integer(kind=intType), dimension(:), allocatable :: sortNumber
        integer(kind=intType), dimension(:), allocatable :: tmpNumber
 
        character(len=maxCGNSNameLen), dimension(:), allocatable :: &
            tmpNames
        logical :: RKExplicit
 
        ! Find out if an explicit RK scheme is used in unsteady mode and
        ! set the logical RKExplicit accordingly. For explicit RK schemes
        ! no residuals are monitored.
 
        rkexplicit = .false.
        if (equationmode == unsteady .and. &
            timeintegrationscheme == explicitrk) rkexplicit = .true.
 
        ! If the turbulent residuals must be monitored add them to the
        ! list of monitoring names. Enough memory should have been
        ! allocated for this.
 
        if (mondturb .and. equations == ransequations .and. &
            (.not. rkexplicit)) then
 
            select case (turbmodel)
 
                ! One equation models of the spalart-allmaras family.
 
            case (spalartallmaras, spalartallmarasedwards)
                nmon = nmon + 1; nmonsum = nmonsum + 1
                monnames(nmon) = cgnsl2resnu
 
                ! Two equation models of the k-w family.
 
            case (komegawilcox, komegamodified, mentersst)
                nmon = nmon + 2; nmonsum = nmonsum + 2
                monnames(nmon - 1) = cgnsl2resk
                monnames(nmon) = cgnsl2resomega
 
                ! Two equation k-tau model.
 
            case (ktau)
                nmon = nmon + 2; nmonsum = nmonsum + 2
                monnames(nmon - 1) = cgnsl2resk
                monnames(nmon) = cgnsl2restau
 
                ! V2f model.
 
            case (v2f)
                nmon = nmon + 4; nmonsum = nmonsum + 4
                monnames(nmon - 3) = cgnsl2resk
                monnames(nmon - 2) = cgnsl2resepsilon
                monnames(nmon - 1) = cgnsl2resv2
                monnames(nmon) = cgnsl2resf
 
            end select
        end if
 
        ! Allocate the memory for sortNumber tmpNumber and tmpNames.
 
        allocate (sortnumber(nmon), tmpnumber(nmon), tmpnames(nmon), &
                  stat=ierr)
        if (ierr /= 0) &
            call terminate("checkMonitor", &
                           "Memory allocation failure for sortNumber, etc.")
 
        ! Loop over the monitoring variables, copy the name into tmpName
        ! and set a number to determine its place in the sequence. If the
        ! variable cannot be monitored for the governing equations, the
        ! priority is set to a high number, such that it will be at the
        ! end of the sorted numbers. At the same time the number of
        ! variables to be monitored, nMonSum, nMonMax and nMon, is
        ! corrected.
 
        nn = nmon
        do i = 1, nn
 
            tmpnames(i) = monnames(i)
 
            ! Determine the place in the sequence for this string.
 
            select case (monnames(i))
            case (cgnsl2resrho)
                sortnumber(i) = 1
                if (rkexplicit) then
                    sortnumber(i) = 10001
                    nmonsum = nmonsum - 1
                end if
 
            case (cgnsl2resmomx)
                sortnumber(i) = 2
                if (rkexplicit) then
                    sortnumber(i) = 10002
                    nmonsum = nmonsum - 1
                end if
 
            case (cgnsl2resmomy)
                sortnumber(i) = 3
                if (rkexplicit) then
                    sortnumber(i) = 10003
                    nmonsum = nmonsum - 1
                end if
 
            case (cgnsl2resmomz)
                sortnumber(i) = 4
                if (rkexplicit) then
                    sortnumber(i) = 10004
                    nmonsum = nmonsum - 1
                end if
 
            case (cgnsl2resrhoe)
                sortnumber(i) = 5
                if (rkexplicit) then
                    sortnumber(i) = 10005
                    nmonsum = nmonsum - 1
                end if
 
            case (cgnsl2resnu)
                sortnumber(i) = 6
                if (equations /= ransequations) then
                    sortnumber(i) = 10001
                    nmonsum = nmonsum - 1
                end if
 
            case (cgnsl2resk)
                sortnumber(i) = 7
                if (equations /= ransequations) then
                    sortnumber(i) = 10002
                    nmonsum = nmonsum - 1
                end if
 
            case (cgnsl2resomega)
                sortnumber(i) = 8
                if (equations /= ransequations) then
                    sortnumber(i) = 10003
                    nmonsum = nmonsum - 1
                end if
 
            case (cgnsl2restau)
                sortnumber(i) = 9
                if (equations /= ransequations) then
                    sortnumber(i) = 10004
                    nmonsum = nmonsum - 1
                end if
 
            case (cgnsl2resepsilon)
                sortnumber(i) = 10
                if (equations /= ransequations) then
                    sortnumber(i) = 10005
                    nmonsum = nmonsum - 1
                end if
 
            case (cgnsl2resv2)
                sortnumber(i) = 11
                if (equations /= ransequations) then
                    sortnumber(i) = 10006
                    nmonsum = nmonsum - 1
                end if
 
            case (cgnsl2resf)
                sortnumber(i) = 12
                if (equations /= ransequations) then
                    sortnumber(i) = 10007
                    nmonsum = nmonsum - 1
                end if
 
            case (cgnscl)
                sortnumber(i) = 101
                if (flowtype == internalflow) then
                    sortnumber(i) = 11001
                    nmonsum = nmonsum - 1
                end if
 
            case (cgnsclp)
                sortnumber(i) = 102
                if (flowtype == internalflow) then
                    sortnumber(i) = 11002
                    nmonsum = nmonsum - 1
                end if
 
            case (cgnsclv)
                sortnumber(i) = 103
                if (equations == eulerequations .or. &
                    flowtype == internalflow) then
                    sortnumber(i) = 11003
                    nmonsum = nmonsum - 1
                end if
 
            case (cgnscd)
                sortnumber(i) = 104
                if (flowtype == internalflow) then
                    sortnumber(i) = 11004
                    nmonsum = nmonsum - 1
                end if
 
            case (cgnscdp)
                sortnumber(i) = 105
                if (flowtype == internalflow) then
                    sortnumber(i) = 11005
                    nmonsum = nmonsum - 1
                end if
 
            case (cgnscdv)
                sortnumber(i) = 106
                if (equations == eulerequations .or. &
                    flowtype == internalflow) then
                    sortnumber(i) = 11006
                    nmonsum = nmonsum - 1
                end if
 
            case (cgnscfx)
                sortnumber(i) = 107
 
            case (cgnscfy)
                sortnumber(i) = 108
 
            case (cgnscfz)
                sortnumber(i) = 109
 
            case (cgnscmx)
                sortnumber(i) = 110
 
            case (cgnscmy)
                sortnumber(i) = 111
 
            case (cgnscmz)
                sortnumber(i) = 112
 
            case ('totalR')
                sortnumber(i) = 113
 
            case (cgnssepsensor)
                sortnumber(i) = 114
 
            case (cgnscavitation)
                sortnumber(i) = 115
 
            case (cgnsaxismoment)
                sortnumber(i) = 116
 
            case (cgnssepsensorksarea)
                sortnumber(i) = 117
 
            case (cgnshdiffmax)
                sortnumber(i) = 201
 
            case (cgnsmachmax)
                sortnumber(i) = 202
 
            case (cgnsyplusmax)
                sortnumber(i) = 203
                if (equations /= ransequations) then
                    sortnumber(i) = 12003
                    nmonmax = nmonmax - 1
                end if
 
            case (cgnseddymax)
                sortnumber(i) = 204
                if (equations /= ransequations) then
                    sortnumber(i) = 12004
                    nmonmax = nmonmax - 1
                end if
 
            case default
                call terminate("checkMonitor", "This should not happen")
            end select
 
        end do
 
        ! Set the new value of nMon, because this might have changed
        ! due to the corrections.
 
        nmon = nmonsum + nmonmax
 
        ! Copy sortNumber in tmpNumber and sort it in increasing order.
        ! Note that here nn must be used and not nMon.
 
        do i = 1, nn
            tmpnumber(i) = sortnumber(i)
        end do
 
        call qsortintegers(sortnumber, nn)
 
        ! Loop over the the number of monitoring variables and store the
        ! new sequence in monNames.
 
        do i = 1, nn
            ii = bsearchintegers(tmpnumber(i), sortnumber)
            monnames(ii) = tmpnames(i)
        end do
 
        ! Release the memory of sortNumber, tmpNumber and tmpNames.
 
        deallocate (sortnumber, tmpnumber, tmpnames, stat=ierr)
        if (ierr /= 0) &
            call terminate("checkMonitor", &
                           "Deallocation error for sortNumber, etc.")
 
        ! Allocate the memory for the monitoring variables.
 
        allocate (monloc(nmon), monglob(nmon), monref(nmon), stat=ierr)
        if (ierr /= 0) &
            call terminate("checkMonitor", &
                           "Memory allocation for monitoring variables")
 
        ! Check if the maximum Mach number or the maximum total enthalpy
        ! difference must be monitored.
 
        monmachorhmax = .false.
        do i = (nmonsum + 1), nmon
            if (monnames(i) == cgnshdiffmax .or. &
                monnames(i) == cgnsmachmax) monmachorhmax = .true.
        end do
 
    end subroutine checkmonitor
    subroutine checkoutput
        !
        !       checkOutput checks and possibly corrects the and output
        !       variables. This depends on the set of governing equations to
        !       be solved.
        !
        use constants
        use extraoutput
        use inputphysics, only: equations, equationmode, walldistanceneeded
        use inputunsteady, only: timeintegrationscheme
        use flowvarrefstate, only: kpresent, eddymodel
        implicit none
 
        ! Determine the governing equations to be solved and set the
        ! variables accordingly.
 
        select case (equations)
        case (eulerequations)
            surfwritecf = .false.
            surfwritech = .false.
            surfwriteyplus = .false.
            surfwritecfx = .false.
            surfwritecfy = .false.
            surfwritecfz = .false.
            surfwriteforceindragdir = .false.
            surfwriteforceinliftdir = .false.
 
            volwritemachturb = .false.
            volwriteeddyvis = .false.
            volwriteratioeddyvis = .false.
            volwritedist = .false.
            volwriteresturb = .false.
 
        case (nsequations)
            surfwriteyplus = .false.
 
            volwritemachturb = .false.
            volwriteeddyvis = .false.
            volwriteratioeddyvis = .false.
            volwritedist = .false.
            volwriteresturb = .false.
 
        case (ransequations)
 
            ! Check if it is possible to write a turbulent Mach
            ! number and eddy viscosity; this depends on the turbulence
            ! model used.
 
            if (.not. kpresent) volwritemachturb = .false.
            if (.not. eddymodel) then
                volwriteeddyvis = .false.
                volwriteratioeddyvis = .false.
            end if
 
            ! If a wall distance free turbulence model is used,
            ! set volWriteDist to .false.
 
            if (.not. walldistanceneeded) volwritedist = .false.
 
        end select
 
        if (equationmode == unsteady .and. &
            timeintegrationscheme == explicitrk) then
            volwriteresrho = .false.
            volwriteresmom = .false.
            volwriteresrhoe = .false.
            volwriteresturb = .false.
        end if
 
    end subroutine checkoutput
    subroutine defaultisoout
        !
        !       defaultIsoOut sets the default set of additional
        !       variables to be written to the solution file; the primitive
        !       variables are always written. This additional set depends on
        !       the governing equations to be solved.
        !
        use constants
        use extraoutput
        use inputphysics, only: equations
        implicit none
 
        ! First set the variables, which are independent from the
        ! governing equations to be solved.
        isowriterho = .false.
        isowritevx = .false.
        isowritevy = .false.
        isowritevz = .false.
        isowritep = .false.
        isowritemx = .false.
        isowritemy = .false.
        isowritemz = .false.
        isowriterhoe = .false.
        isowritetemp = .false.
        isowritecp = .false.
        isowritemach = .false.
        isowritemachturb = .false.
        isowritedist = .false.
        isowritevort = .false.
        isowritevortx = .false.
        isowritevorty = .false.
        isowritevortz = .false.
        isowriteptotloss = .false.
        isowriteresrho = .false.
        isowriteresmom = .false.
        isowriteresrhoe = .false.
        isowriteshock = .false.
        isowritefilteredshock = .false.
        ! Set the values which depend on the equations to be solved.
 
        select case (equations)
        case (eulerequations)
            isowriteeddyvis = .false.
            isowriteratioeddyvis = .false.
            isowriteresturb = .false.
 
        case (nsequations)
            isowriteeddyvis = .false.
            isowriteratioeddyvis = .false.
            isowriteresturb = .false.
 
        case (ransequations)
            isowriteeddyvis = .false.
            isowriteratioeddyvis = .false.
            isowriteresturb = .false.
        end select
 
    end subroutine defaultisoout
    subroutine defaultmonitor
        !
        !       defaultMonitor sets the default set of variables to be
        !       monitored during the convergence. This set depends on the
        !       governing equations to be solved.
        !
        use constants
        use cgnsnames
        use inputphysics, only: equations, flowtype
        use monitor, only: nmonsum, nmonmax, nmon, monnames, showcpu
        use utils, only: terminate
        implicit none
        !
        !      Local variables.
        !
        integer :: ierr
 
        ! CPU time is written to stdout.
 
        showcpu = .true.
 
        ! Determine the governing equations to be solved.
 
        select case (equations)
        case (eulerequations)
 
            ! Set the number of summation and maximum monitor variables
            ! and allocate the memory for the monitoring names.
            ! A distinction is made between internal and external flows,
            ! because cl and cd do not make a lot of sense for the former.
 
            if (flowtype == internalflow) then
 
                ! Internal flow; only the density residual is monitored.
 
                nmonsum = 1; nmonmax = 0; nmon = 1
                allocate (monnames(nmon), stat=ierr)
                if (ierr /= 0) &
                    call terminate("defaultMonitor", &
                                   "Memory allocation failure for monNames")
 
                ! Set the names for the variables to be monitored.
 
                monnames(1) = cgnsl2resrho
 
            else
 
                ! External; also lift and drag is monitored.
 
                nmonsum = 3; nmonmax = 0; nmon = 3
                allocate (monnames(nmon), stat=ierr)
                if (ierr /= 0) &
                    call terminate("defaultMonitor", &
                                   "Memory allocation failure for monNames")
 
                ! Set the names for the variables to be monitored.
 
                monnames(1) = cgnsl2resrho
                monnames(2) = cgnscl
                monnames(3) = cgnscd
 
            end if
 
        case (nsequations)
 
            ! Set the number of summation and maximum monitor variables
            ! and allocate the memory for the monitoring names.
            ! A distinction is made between internal and external flows,
            ! because cl and cd do not make a lot of sense for the former.
 
            if (flowtype == internalflow) then
 
                ! Internal flow; only the density residual is monitored.
 
                nmonsum = 1; nmonmax = 0; nmon = 1
                allocate (monnames(nmon), stat=ierr)
                if (ierr /= 0) &
                    call terminate("defaultMonitor", &
                                   "Memory allocation failure for monNames")
 
                ! Set the names for the variables to be monitored.
 
                monnames(1) = cgnsl2resrho
 
            else
 
                ! External; also lift and drag (total and viscous)
                ! is monitored.
 
                nmonsum = 4; nmonmax = 0; nmon = 4
                allocate (monnames(nmon), stat=ierr)
                if (ierr /= 0) &
                    call terminate("defaultMonitor", &
                                   "Memory allocation failure for monNames")
 
                ! Set the names for the variables to be monitored.
 
                monnames(1) = cgnsl2resrho
                monnames(2) = cgnscl
                monnames(3) = cgnscd
                monnames(4) = cgnscdv
 
            end if
 
        case (ransequations)
 
            ! Set the number of summation and maximum monitor variables
            ! and allocate the memory for the monitoring names.
            ! A distinction is made between internal and external flows,
            ! because cl and cd do not make a lot of sense for the former.
 
            if (flowtype == internalflow) then
 
                ! Internal flow; the density residual as well as the
                ! maximum values of yplus and the eddy viscosity ration
                ! are monitored.
 
                nmonsum = 1; nmonmax = 2; nmon = 3
                allocate (monnames(nmon), stat=ierr)
                if (ierr /= 0) &
                    call terminate("defaultMonitor", &
                                   "Memory allocation failure for monNames")
 
                ! Set the names for the variables to be monitored.
 
                monnames(1) = cgnsl2resrho
                monnames(2) = cgnsyplusmax
                monnames(3) = cgnseddymax
 
            else
 
                ! External; also lift and drag (total and viscous)
                ! is monitored.
 
                nmonsum = 4; nmonmax = 2; nmon = 6
                allocate (monnames(nmon), stat=ierr)
                if (ierr /= 0) &
                    call terminate("defaultMonitor", &
                                   "Memory allocation failure for monNames")
 
                ! Set the names for the variables to be monitored.
 
                monnames(1) = cgnsl2resrho
                monnames(2) = cgnscl
                monnames(3) = cgnscd
                monnames(4) = cgnscdv
                monnames(5) = cgnsyplusmax
                monnames(6) = cgnseddymax
 
            end if
 
        end select
 
    end subroutine defaultmonitor
    subroutine defaultsurfaceout
        !
        !       defaultSurfaceOut sets the default set of surface variables
        !       to be written to the solution file. This set depends on the
        !       governing equations to be solved.
        !
        use constants
        use extraoutput
        use inputphysics, only: equations
        implicit none
 
        ! First set the variables, which are independent from the
        ! governing equations to be solved.
 
        surfwriterho = .true.
        surfwritep = .false.
        surfwritetemp = .false.
        surfwritevx = .true.
        surfwritevy = .true.
        surfwritevz = .true.
        surfwritecp = .true.
        surfwritemach = .true.
 
        ! Set the values which depend on the equations to be solved.
 
        select case (equations)
        case (eulerequations)
            surfwriteptotloss = .true.
            surfwritecf = .false.
            surfwritech = .false.
            surfwriteyplus = .false.
            surfwritecfx = .false.
            surfwritecfy = .false.
            surfwritecfz = .false.
 
        case (nsequations)
            surfwriteptotloss = .false.
            surfwritecf = .true.
            surfwritech = .false.
            surfwriteyplus = .false.
            surfwritecfx = .true.
            surfwritecfy = .true.
            surfwritecfz = .true.
 
        case (ransequations)
            surfwriteptotloss = .false.
            surfwritecf = .true.
            surfwritech = .false.
            surfwriteyplus = .true.
            surfwritecfx = .true.
            surfwritecfy = .true.
            surfwritecfz = .true.
        end select
 
    end subroutine defaultsurfaceout
    subroutine defaultvolumeout
        !
        !       defaultVolumeOut sets the default set of additional
        !       variables to be written to the solution file; the primitive
        !       variables are always written. This additional set depends on
        !       the governing equations to be solved.
        !
        use constants
        use extraoutput
        use inputphysics, only: equations
        implicit none
 
        ! First set the variables, which are independent from the
        ! governing equations to be solved.
 
        volwritemx = .false.
        volwritemy = .false.
        volwritemz = .false.
        volwriterhoe = .false.
        volwritetemp = .false.
        volwritecp = .false.
        volwritemach = .false.
        volwritemachturb = .false.
        volwritedist = .false.
        volwritevort = .false.
        volwritevortx = .false.
        volwritevorty = .false.
        volwritevortz = .false.
        volwriteptotloss = .true.
        volwriteresrho = .true.
        volwriteresmom = .false.
        volwriteresrhoe = .false.
 
        ! Set the values which depend on the equations to be solved.
 
        select case (equations)
        case (eulerequations)
            volwriteeddyvis = .false.
            volwriteratioeddyvis = .false.
            volwriteresturb = .false.
 
        case (nsequations)
            volwriteeddyvis = .false.
            volwriteratioeddyvis = .false.
            volwriteresturb = .false.
 
        case (ransequations)
            volwriteeddyvis = .true.
            volwriteratioeddyvis = .true.
            volwriteresturb = .true.
        end select
 
    end subroutine defaultvolumeout
 
    subroutine dummyreadparamfile
        !
        !       This subroutine is the same as readParamFile EXCEPT it does not
        !       read the actual file. Values are set diectly from python for
        !       all the options and then this file is run.
        !
        use constants
        use inputphysics, only: cpmodel
 
        implicit none
        !
        !      Local variables
        !
        integer, parameter :: readUnit = 32
 
        ! Check if all the desired input parameters were specified and
        ! print warnings if some irrelevant ones are specified.
 
        call checkinputparam
 
        ! Read the cp curve fits from file if a variable cp model
        ! must be used.
 
        if (cpmodel == cptempcurvefits) call readcptempcurvefits
 
        ! Determine the number of governing equations and set the
        ! corresponding parameters accordingly.
 
        call setequationparameters
 
        ! Extract the multigrid info from the string.
 
        call extractmginfo
 
        ! If no monitoring variables were specified, set the default set.
        ! Idem for the surface and volume output variables.
 
        if (.not. monitorspecified) call defaultmonitor
        if (.not. surfaceoutspecified) call defaultsurfaceout
        if (.not. volumeoutspecified) call defaultvolumeout
        if (.not. isooutspecified) call defaultisoout
 
        ! Check the monitoring and output variables.
 
        call checkmonitor
        call checkoutput
 
    end subroutine dummyreadparamfile
    subroutine extractmginfo
        !
        !       extractMgInfo creates the integer array cycleStrategy from
        !       the string describing the multigrid strategy. This string
        !       either contains a predefined strategy, like sg, 2v, 4w, etc.,
        !       or a combination of -1's, 0's and 1's, which defines a user
        !       defined strategy. The integers -1, 0 and 1 have the following
        !       meaning:  0 -> perform an iteration step on the current grid.
        !                 1 -> go to next coarser grid.
        !                -1 -> go to next finer grid.
        !       For a valid cycling strategy the sum of the elements of the
        !       array should be 0.
        !
        use constants
        use inputiteration, only: cyclestrategy, nmglevels, nmgsteps, mgstartlevel, mgdescription
        use inputphysics, only: equationmode
        use inputunsteady, only: timeintegrationscheme
        use communication, only: myid
        use utils, only: converttolowercase, terminate
        implicit none
        !
        !      Local variables
        !
        integer :: stringLen, error
        integer(kind=intType) :: i, ii, nMinus, nn
        character(len=maxStringLen) :: errorMessage
 
        ! For an unsteady computation using explicit Runge-Kutta schemes
        ! overrule mgDescription to sg.
 
        if (equationmode == unsteady .and. &
            timeintegrationscheme == explicitrk) mgdescription = "sg"
 
        ! Create a lower case version of mgDescription and determine the
        ! length of the string. Note that if the string contains a user
        ! defined cycling strategy the contents (0's, 1's, -1's) is not
        ! changed by the call to convertToLowerCase
 
        mgdescription = trim(mgdescription)
        call converttolowercase(mgdescription)
        stringlen = len_trim(mgdescription)
 
        ! Check for predefined cycling strategies.
 
        if (mgdescription == "sg") then
 
            ! Single grid computation. Set the values for the parameters
            ! nMGSteps (number of steps in the cycle strategy) and
            ! cycleStrategy (the cycle strategy itself).
 
            nmgsteps = 1
            allocate (cyclestrategy(1), stat=error)
            if (error /= 0) call terminate("extractMgInfo", &
                                           "allocation error for 1 integer")
            cyclestrategy(1) = 0
 
            ! Set the number of grid levels needed by the multigrid to 1.
 
            nmglevels = 1
 
        else if (mgdescription(stringlen:stringlen) == "v") then
 
            ! Must be a v-cycle. The rest of the string should only contain
            ! digits. Check this.
 
            if (.not. digitsonlyinstring(mgdescription(:stringlen - 1))) then
                write (errormessage, *) "Invalid cycle strategy, ", &
                    mgdescription(:stringlen), ", specified"
                if (myid == 0) call terminate("extractMgInfo", errormessage)
            end if
 
            ! Read the number of levels in the cycle.
 
            read (mgdescription(:stringlen - 1), *) nmglevels
 
            ! Determine the number of steps in cycleStrategy and allocate
            ! the memory for it.
 
            nmgsteps = 4 * nmglevels - 4
            allocate (cyclestrategy(nmgsteps), stat=error)
            if (error /= 0) then
                write (errormessage, *) "Allocation error for", nmgsteps, &
                    "integers for the v-cycle ", &
                    mgdescription(:stringlen)
                call terminate("extractMgInfo", errormessage)
            end if
 
            ! Set the values of cycleStrategy.
 
            ii = 1
            do i = 1, (nmglevels - 1)
                cyclestrategy(ii) = 0
                cyclestrategy(ii + 1) = 1
                ii = ii + 2
            end do
 
            do i = 1, (nmglevels - 1)
                cyclestrategy(ii) = 0
                cyclestrategy(ii + 1) = -1
                ii = ii + 2
            end do
 
        else if (mgdescription(stringlen:stringlen) == "w") then
 
            ! Must be a w-cycle. The rest of the string should only contain
            ! digits. Check this.
 
            if (.not. digitsonlyinstring(mgdescription(:stringlen - 1))) then
                write (errormessage, *) "Invalid cycle strategy, ", &
                    mgdescription(:stringlen), ", specified"
                if (myid == 0) call terminate("extractMgInfo", errormessage)
            end if
 
            ! Read the number of levels in the cycle.
 
            read (mgdescription(:stringlen - 1), *) nmglevels
 
            ! Determine the number of steps in cycleStrategy and allocate
            ! the memory for it.
 
            nmgsteps = computenstepswcycle(nmglevels)
            allocate (cyclestrategy(nmgsteps), stat=error)
            if (error /= 0) then
                write (errormessage, *) "Allocation error for", nmgsteps, &
                    "integers for the w-cycle ", &
                    mgdescription(:stringlen)
                call terminate("extractMgInfo", errormessage)
            end if
 
            ! Set the values of cycleStrategy.
 
            ii = 1
            call setentrieswcycle(ii, nmglevels)
 
        else
 
            ! The string must be a collection of 0's, -1's and 1's to
            ! describe the cycle strategy. Get rid of the internal spaces
            ! first and determine the amount of -'s.
 
            ii = 0
            nminus = 0
            do i = 1, stringlen
                if (mgdescription(i:i) /= " ") then
                    ii = ii + 1
                    mgdescription(ii:ii) = mgdescription(i:i)
                    if (mgdescription(ii:ii) == "-") nminus = nminus + 1
                end if
            end do
            stringlen = ii
 
            ! Determine the number of steps in the cycle strategy and
            ! allocate the memory for it.
 
            nmgsteps = ii - nminus
            allocate (cyclestrategy(nmgsteps), stat=error)
            if (error /= 0) then
                write (errormessage, *) "Allocation error for", nmgsteps, &
                    "integers for the cycle strategy"
                call terminate("extractMgInfo", errormessage)
            end if
 
            ! Determine the entries for cycleStrategy.
 
            i = 1
            nn = 1
            do
                ii = i
                if (mgdescription(i:i) == "-") i = i + 1
 
                ! Determine the case we are having here.
 
                select case (mgdescription(ii:i))
                case ("0")
                    cyclestrategy(nn) = 0
                case ("1")
                    cyclestrategy(nn) = 1
                case ("-1")
                    cyclestrategy(nn) = -1
                case default
                    write (errormessage, *) "Invalid character, ", &
                         mgdescription(ii:i), &
                         ", in the string describing &
                         &cycling strategy"
                    if (myid == 0) call terminate("extractMgInfo", errormessage)
                end select
 
                ! Update i and nn
 
                i = i + 1
                nn = nn + 1
 
                ! Exit the do loop in case i is larger than stringLen.
 
                if (i > stringlen) exit
            end do
 
            ! Check if the string specified is valid and determine the
            ! maximum grid level needed in the cycle.
 
            nn = 0
            nmglevels = 0
            do i = 1, nmgsteps
                nn = nn + cyclestrategy(i)
                nmglevels = max(nn, nmglevels)
            end do
            nmglevels = nmglevels + 1
 
            if (nn /= 0 .and. myid == 0) &
                call terminate("extractMgInfo", &
                               "sum of coefficients in cycle strategy is not 0")
        end if
 
        ! Correct the value of mgStartlevel in case a nonpositive number
        ! has been specified. In that case it is set to -1, the default
        ! value.
 
        if (mgstartlevel <= 0) mgstartlevel = -1
 
        ! Determine the value of mgStartlevel. This parameter might be
        ! specified in the python script file and is checked here for
        ! consistency. If mgStartlevel has not been specified to any
        ! specific level it is set to the coarsest level in the mg cycle
        ! (starting from free stream) or to the finest level (restart).
        ! The restart is handled in python wrapper.
 
        if (mgstartlevel == -1) then
 
            ! Value has not been specified. Default value is set, see
            ! the comments above.
            mgstartlevel = nmglevels
 
        end if
 
    end subroutine extractmginfo
 
    !      ==================================================================
 
    logical function digitsonlyinstring(string)
        !
        !       digitsOnlyInString checks whether the given string contains
        !       digits only or if other character types are present. In the
        !       former case the function returns .True., otherwise .False.
        !
        implicit none
        !
        !      Subroutine argument                                              *
        !
        character(len=*), intent(in) :: string
        !
        !      Local variables
        !
        integer :: i, stringlen
 
        ! Initialize digitsOnlyInString to .True.
 
        digitsonlyinstring = .true.
 
        ! Determine the length of the string.
 
        stringlen = len_trim(string)
 
        ! Loop over the elements of the string and check if they are digits.
 
        do i = 1, stringlen
            if (string(i:i) < "0" .or. string(i:i) > "9") &
                digitsonlyinstring = .false.
        end do
 
    end function digitsonlyinstring
 
    !      ==================================================================
 
    recursive function computenstepswcycle(nLevels) result(nSteps)
        !
        !       computeNstepsWcycle is recursive function, which determines
        !       the number of entries of a w-cycle of a given level.
        !
        use constants
        use communication
        use utils, only: terminate
        implicit none
        !
        !      Result variable
        !
        integer(kind=intType) :: nsteps
        !
        !      Function argument
        !
        integer(kind=intType), intent(in) :: nlevels
        !
        !      Local variables
        !
        character(len=maxStringLen) :: errormessage
 
        ! Determine the case we are having here. For nLevels is less
        ! than 2 an error message is printed, in case nLevels is 2
        ! the recursion is broken and otherwise a recursive call is made.
 
        if (nlevels < 2) then
            write (errormessage, *) "Wrong value of nLevels", nlevels
            if (myid == 0) call terminate("computeNstepsWcycle", errormessage)
        else if (nlevels == 2) then
            nsteps = 4
        else
            nsteps = 4 + 2 * computenstepswcycle(nlevels - 1)
        end if
 
    end function computenstepswcycle
 
    !      ==================================================================
 
    recursive subroutine setentrieswcycle(counter, nLevels)
        !
        !       setEntriesWcycle is a recursive subroutine, which actually
        !       fills the entries of cycleStrategy for a w-cycle.
        !
        use constants
        use inputiteration, only: cyclestrategy
        use communication, only: myid
        use utils, only: terminate
        implicit none
        !
        !      Subroutine argument.
        !
        integer(kind=intType), intent(inout) :: counter
        integer(kind=intType), intent(in) :: nlevels
        !
        !      Local variables
        !
        character(len=maxStringLen) :: errormessage
 
        ! Determine the case we are having here. For nLevels is less
        ! than 2 an error message is printed, in case nLevels is 2
        ! the recursion is broken and otherwise a recursive call is made.
 
        if (nlevels < 2) then
 
            write (errormessage, *) "Wrong value of nLevels", nlevels
            if (myid == 0) call terminate("setEntriesWcycle", errormessage)
 
        else if (nlevels == 2) then
 
            cyclestrategy(counter) = 0
            cyclestrategy(counter + 1) = 1
            cyclestrategy(counter + 2) = 0
            cyclestrategy(counter + 3) = -1
 
            counter = counter + 4
 
        else
 
            cyclestrategy(counter) = 0
            cyclestrategy(counter + 1) = 1
            counter = counter + 2
 
            call setentrieswcycle(counter, nlevels - 1)
            call setentrieswcycle(counter, nlevels - 1)
 
            cyclestrategy(counter) = 0
            cyclestrategy(counter + 1) = -1
            counter = counter + 2
 
        end if
 
    end subroutine setentrieswcycle
    subroutine isovariables(variables)
        !
        !       isoVariables extracts from the given string the extra
        !       iso surface variables to be written to the solution file.
        !
        use constants
        use extraoutput
        use utils, only: converttolowercase, terminate
        implicit none
        !
        !      Subroutine arguments.
        !
        character(len=*), intent(inout) :: variables
        !
        !      Local variables.
        !
        integer :: nVarSpecified, pos
 
        character(len=15) :: keyword
        character(len=maxStringLen) :: errorMessage
 
        ! Convert the string variables to lower case.
 
        call converttolowercase(variables)
 
        ! Initialize all the iso output variables to .False.
        isowriterho = .false.
        isowritevx = .false.
        isowritevy = .false.
        isowritevz = .false.
        isowritep = .false.
        isowriteturb = .false.
 
        isowritemx = .false.
        isowritemy = .false.
        isowritemz = .false.
        isowriterhoe = .false.
        isowritetemp = .false.
        isowritevort = .false.
        isowritevortx = .false.
        isowritevorty = .false.
        isowritevortz = .false.
 
        isowritecp = .false.
        isowritemach = .false.
        isowritemachturb = .false.
        isowriteptotloss = .false.
 
        isowriteeddyvis = .false.
        isowriteratioeddyvis = .false.
        isowritedist = .false.
 
        isowriteresrho = .false.
        isowriteresmom = .false.
        isowriteresrhoe = .false.
        isowriteresturb = .false.
 
        isowriteshock = .false.
        isowritefilteredshock = .false.
 
        isowriteblank = .false.
 
        ! Initialize nVarSpecified to 0. This serves as a test
        ! later on.
 
        nvarspecified = 0
 
        ! Loop to extract the info from the string variables.
 
        do
            ! Condition to exit the loop.
 
            if (len_trim(variables) == 0) exit
 
            ! Locate the first occurance of the _ in the string and
            ! determine the string keyword.
 
            pos = index(variables, "_")
            if (pos == 0) then
                keyword = variables
                variables = ""
            else
                keyword = variables(:pos - 1)
                variables = variables(pos + 1:)
            end if
 
            ! Check the keyword.
 
            select case (keyword)
            case ("")
                ! Multiple occurence of "_". Just ignore it.
 
            case ("rho")
                isowriterho = .true.
                nvarspecified = nvarspecified + 1
 
            case ("vx")
                isowritevx = .true.
                nvarspecified = nvarspecified + 1
 
            case ("vy")
                isowritevy = .true.
                nvarspecified = nvarspecified + 1
 
            case ("vz")
                isowritevz = .true.
                nvarspecified = nvarspecified + 1
 
            case ("P")
                isowritep = .true.
                nvarspecified = nvarspecified + 1
 
            case ("turb")
                isowriteturb = .true.
                nvarspecified = nvarspecified + 1
 
            case ("mx")
                isowritemx = .true.
                nvarspecified = nvarspecified + 1
 
            case ("my")
                isowritemy = .true.
                nvarspecified = nvarspecified + 1
 
            case ("mz")
                isowritemz = .true.
                nvarspecified = nvarspecified + 1
 
            case ("rvx")
                isowritervx = .true.
                nvarspecified = nvarspecified + 1
 
            case ("rvy")
                isowritervy = .true.
                nvarspecified = nvarspecified + 1
 
            case ("rvz")
                isowritervz = .true.
                nvarspecified = nvarspecified + 1
 
            case ("rhoe")
                isowriterhoe = .true.
                nvarspecified = nvarspecified + 1
 
            case ("temp")
                isowritetemp = .true.
                nvarspecified = nvarspecified + 1
 
            case ("vort")
                isowritevort = .true.
                nvarspecified = nvarspecified + 1
 
            case ("vortx")
                isowritevortx = .true.
                nvarspecified = nvarspecified + 1
 
            case ("vorty")
                isowritevorty = .true.
                nvarspecified = nvarspecified + 1
 
            case ("vortz")
                isowritevortz = .true.
                nvarspecified = nvarspecified + 1
 
            case ("cp")
                isowritecp = .true.
                nvarspecified = nvarspecified + 1
 
            case ("mach")
                isowritemach = .true.
                nvarspecified = nvarspecified + 1
 
            case ("rmach")
                isowritermach = .true.
                nvarspecified = nvarspecified + 1
 
            case ("macht")
                isowritemachturb = .true.
                nvarspecified = nvarspecified + 1
 
            case ("ptloss")
                isowriteptotloss = .true.
                nvarspecified = nvarspecified + 1
 
            case ("eddy")
                isowriteeddyvis = .true.
                nvarspecified = nvarspecified + 1
 
            case ("eddyratio")
                isowriteratioeddyvis = .true.
                nvarspecified = nvarspecified + 1
 
            case ("dist")
                isowritedist = .true.
                nvarspecified = nvarspecified + 1
 
            case ("resrho")
                isowriteresrho = .true.
                nvarspecified = nvarspecified + 1
 
            case ("resmom")
                isowriteresmom = .true.
                nvarspecified = nvarspecified + 1
 
            case ("resrhoe")
                isowriteresrhoe = .true.
                nvarspecified = nvarspecified + 1
 
            case ("resturb")
                isowriteresturb = .true.
                nvarspecified = nvarspecified + 1
 
            case ("blank")
                isowriteblank = .true.
                nvarspecified = nvarspecified + 1
 
            case ("shock")
                isowriteshock = .true.
                nvarspecified = nvarspecified + 1
 
            case ("filteredshock")
                isowritefilteredshock = .true.
                nvarspecified = nvarspecified + 1
 
            case default
                pos = len_trim(keyword)
                write (errormessage, "(3a)") "Unknown extra iso output &
                     &variable, ", trim(keyword), &
                     ", specified"
                call terminate("isoVariables", errormessage)
 
            end select
 
        end do
 
        ! Set this to true regardless...it is possible no varibles were
        ! specified
        isooutspecified = .true.
 
    end subroutine isovariables
 
    subroutine monitorvariables(variables)
        !
        !       monitorVariables extracts from the given string the variables
        !       to be monitored during the convergence.
        !
        use constants
        use cgnsnames
        use communication, only: myid, adflow_comm_world
        use monitor, only: monnames, nmon, nmonmax, nmonsum, showcpu
        use utils, only: converttolowercase, terminate
        implicit none
        !
        !      Subroutine arguments.
        !
        character(len=*), intent(inout) :: variables
        !
        !      Local parameter.
        !
        integer(kind=intType), parameter :: nVarMax = 21
        !
        !      Local variables.
        !
        integer :: pos, ierr
 
        character(len=15) :: keyword
        character(len=maxStringLen) :: errorMessage
 
        character(len=maxCGNSNameLen), dimension(nVarMax) :: tmpNames
 
        logical :: monDrho, monTotalR
 
        ! Check if the monitoring names have already been allocated.
        ! This happens when multiple lines for the monitoring variables
        ! are specified in the parameter file. If this happens the last
        ! value is taken and thus release the memory of previously
        ! specified names.
 
        if (allocated(monnames)) then
            deallocate (monnames, stat=ierr)
            if (ierr /= 0) call terminate("monitorVariables", &
                                          "Deallocation error for monNames")
        end if
 
        ! Initialize monDrho, monDturb and showCPU to .false.
 
        mondrho = .false.
        montotalr = .false.
        mondturb = .false.
        showcpu = .false.
 
        ! Initialize nMonSum, nMonMax and nMon to 0.
 
        nmonsum = 0
        nmonmax = 0
        nmon = 0
 
        ! Convert the string variables to lower case.
 
        call converttolowercase(variables)
 
        ! Loop to extract the info from the string variables.
 
        do
            ! Condition to exit the loop.
 
            if (len_trim(variables) == 0) exit
 
            ! Locate the first occurance of the _ in the string and
            ! determine the string keyword.
 
            pos = index(variables, "_")
            if (pos == 0) then
                keyword = variables
                variables = ""
            else
                keyword = variables(:pos - 1)
                variables = variables(pos + 1:)
            end if
 
            ! Check the keyword.
 
            select case (keyword)
            case ("")
                ! Multiple occurence of "_". Just ignore it.
 
            case ("cpu")    ! only written to stdout.
                showcpu = .true.
 
            case ("resrho")
                mondrho = .true.
                nmon = nmon + 1; nmonsum = nmonsum + 1
                tmpnames(nmon) = cgnsl2resrho
 
            case ("resmom")
                nmon = nmon + 3; nmonsum = nmonsum + 3
                tmpnames(nmon - 2) = cgnsl2resmomx
                tmpnames(nmon - 1) = cgnsl2resmomy
                tmpnames(nmon) = cgnsl2resmomz
 
            case ("resrhoe")
                nmon = nmon + 1; nmonsum = nmonsum + 1
                tmpnames(nmon) = cgnsl2resrhoe
 
            case ("resturb")  ! special case, because the turbulence model
                ! Is not yet known. See checkMonitor.
                mondturb = .true.
 
            case ("cl")
                nmon = nmon + 1; nmonsum = nmonsum + 1
                tmpnames(nmon) = cgnscl
 
            case ("clp")
                nmon = nmon + 1; nmonsum = nmonsum + 1
                tmpnames(nmon) = cgnsclp
 
            case ("clv")
                nmon = nmon + 1; nmonsum = nmonsum + 1
                tmpnames(nmon) = cgnsclv
 
            case ("cd")
                nmon = nmon + 1; nmonsum = nmonsum + 1
                tmpnames(nmon) = cgnscd
 
            case ("cdp")
                nmon = nmon + 1; nmonsum = nmonsum + 1
                tmpnames(nmon) = cgnscdp
 
            case ("cdv")
                nmon = nmon + 1; nmonsum = nmonsum + 1
                tmpnames(nmon) = cgnscdv
 
            case ("cfx")
                nmon = nmon + 1; nmonsum = nmonsum + 1
                tmpnames(nmon) = cgnscfx
 
            case ("cfy")
                nmon = nmon + 1; nmonsum = nmonsum + 1
                tmpnames(nmon) = cgnscfy
 
            case ("cfz")
                nmon = nmon + 1; nmonsum = nmonsum + 1
                tmpnames(nmon) = cgnscfz
 
            case ("cmx")
                nmon = nmon + 1; nmonsum = nmonsum + 1
                tmpnames(nmon) = cgnscmx
 
            case ("cmy")
                nmon = nmon + 1; nmonsum = nmonsum + 1
                tmpnames(nmon) = cgnscmy
 
            case ("cmz")
                nmon = nmon + 1; nmonsum = nmonsum + 1
                tmpnames(nmon) = cgnscmz
 
            case ("hdiff")
                nmon = nmon + 1; nmonmax = nmonmax + 1
                tmpnames(nmon) = cgnshdiffmax
 
            case ("mach")
                nmon = nmon + 1; nmonmax = nmonmax + 1
                tmpnames(nmon) = cgnsmachmax
 
            case ("yplus")
                nmon = nmon + 1; nmonmax = nmonmax + 1
                tmpnames(nmon) = cgnsyplusmax
 
            case ("eddyv")
                nmon = nmon + 1; nmonmax = nmonmax + 1
                tmpnames(nmon) = cgnseddymax
 
            case ("totalr")
                montotalr = .true.
                nmon = nmon + 1; nmonsum = nmonsum + 1
                tmpnames(nmon) = 'totalR'
 
            case ("sepsensor")
                nmon = nmon + 1; nmonsum = nmonsum + 1
                tmpnames(nmon) = cgnssepsensor
 
            case ("SepSensorKsArea")
                nmon = nmon + 1; nmonsum = nmonsum + 1
                tmpnames(nmon) = cgnssepsensorksarea
 
            case ("cavitation")
                nmon = nmon + 1; nmonsum = nmonsum + 1
                tmpnames(nmon) = cgnscavitation
 
            case ("axismoment")
                nmon = nmon + 1; nmonsum = nmonsum + 1
                tmpnames(nmon) = cgnsaxismoment
 
            case default
                write (errormessage, "(3a)") "Unknown monitoring variable, ", &
                    trim(keyword), ", specified"
                if (myid == 0) &
                    call terminate("monitorVariables", errormessage)
                call mpi_barrier(adflow_comm_world, ierr)
 
            end select
 
        end do
 
        ! If the density residual was not specified to be monitored,
        ! add it to tmpNames.
 
        if (.not. mondrho) then
            nmon = nmon + 1; nmonsum = nmonsum + 1
            tmpnames(nmon) = cgnsl2resrho
        end if
 
        ! If the total residual was not specified to be monitored, add
        ! it to tmpNames.
 
        if (.not. montotalr) then
            nmon = nmon + 1; nmonsum = nmonsum + 1
            tmpnames(nmon) = "totalR"
        end if
 
        ! Allocate the memory for monNames. If the turbulent residuals
        ! must be monitored allocate some extra place.
 
        pos = nmon
        if (mondturb) pos = nmon + 4
        allocate (monnames(pos), stat=ierr)
        if (ierr /= 0) &
            call terminate("monitorVariables", &
                           "Memory allocation failure for monNames")
 
        ! Copy the monitoring names into monNames.
 
        do pos = 1, nmon
            monnames(pos) = tmpnames(pos)
        end do
 
        ! Set monitorSpecified to .true. to indicate that monitoring
        ! variables have been specified.
 
        monitorspecified = .true.
 
    end subroutine monitorvariables
    subroutine readcptempcurvefits
        !
        !       readCpTempCurveFits reads the curve fits for the cp as a
        !       function of the temperature from the file cpFile.
        !
        use constants
        use communication, only: myid, adflow_comm_world
        use cpcurvefits, only: cptrange, cptempfit, cpeint, cphint, cvn, cv0, &
                               cpnparts
        use inputio, only: cpfile
        use utils, only: terminate
        implicit none
 
        ! Local variables.
 
        integer, parameter :: readUnit = 32
 
        integer :: ios, ierr
 
        integer(kind=intType) :: nn, mm, kk, ii
        real(kind=realtype) :: t1, t2, e0
 
        character(len=2*maxStringLen) :: errorMessage
        character(len=512) :: string
 
        ! Open the file for reading and check if it went okay. If the file
        ! is not found, processor 0 prints an error message.
 
        open (unit=readunit, file=cpfile, status="old", &
              action="read", iostat=ios)
 
        if (ios /= 0) then
 
            write (errormessage, *) "Cp curve fit file ", trim(cpfile), &
                " not found."
            if (myid == 0) &
                call terminate("readCpTempCurveFits", errormessage)
 
            call mpi_barrier(adflow_comm_world, ierr)
        end if
 
        ! Skip the comment lines and read the number of parts.
        ! Check if a valid number is read.
 
        call findnextinfoline(readunit, string)
        read (string, *) cpnparts
 
        if (cpnparts <= 0) then
            if (myid == 0) &
                 call terminate("readCpTempCurveFits", &
                 "Wrong number of temperature ranges in &
                 &Cp curve fit file.")
            call mpi_barrier(adflow_comm_world, ierr)
        end if
 
        ! Allocate the memory for the variables to store the curve fit
        ! data.
 
        allocate (cptrange(0:cpnparts), cpeint(0:cpnparts), &
                  cphint(0:cpnparts), cptempfit(cpnparts), &
                  stat=ierr)
        if (ierr /= 0) &
             call terminate("readCpTempCurveFits", &
             "Memory allocation failure for cpTrange, &
             &cpEint, cpHint and cpTempFit")
 
        ! Loop over the number of temperature ranges.
 
        nranges: do nn = 1, cpnparts
 
            ! Find the next line with information and read the temperature
            ! range.
 
            call findnextinfoline(readunit, string)
            read (string, *) t1, t2
 
            ! If this is the first range, set the temperature range;
            ! otherwise check if the lower boundary equals the upper
            ! boundary of the previous range.
 
            if (nn == 1) then
                cptrange(0) = t1
                cptrange(1) = t2
            else
                cptrange(nn) = t2
 
                if (t1 /= cptrange(nn - 1)) then
                    if (myid == 0) &
                        call terminate("readCpTempCurveFits", &
                                       "Curve fit boundary not continuous")
                    call mpi_barrier(adflow_comm_world, ierr)
                end if
            end if
 
            ! Read the number of points in the fit.
 
            call findnextinfoline(readunit, string)
            read (string, *) cptempfit(nn)%nterm
 
            ! Allocate the memory for the exponents and the constants.
 
            ii = cptempfit(nn)%nterm
            allocate (cptempfit(nn)%exponents(ii), &
                      cptempfit(nn)%constants(ii), stat=ierr)
            if (ierr /= 0) &
                 call terminate("readCpTempCurveFits", &
                 "Memory allocation failure for exponents and &
                 &constants")
 
            ! Read the exponents from the file.
 
            call findnextinfoline(readunit, string)
            do ii = 1, cptempfit(nn)%nterm
 
                ! Read the exponent from the string.
 
                read (string, *) cptempfit(nn)%exponents(ii)
 
                ! Remove this value from the string if this is not the
                ! last exponent to be read.
 
                if (ii < cptempfit(nn)%nterm) then
                    ios = index(string, " ")
                    if (ios > 0) then
                        string = string(ios:)
                        string = adjustl(string)
                        string = trim(string)
                    else
                        if (myid == 0) &
                             call terminate("readCpTempCurveFits", &
                             "Not enough exponents on line; &
                             &Cp curve fit file not valid.")
                        call mpi_barrier(adflow_comm_world, ierr)
                    end if
                end if
            end do
 
            ! Read the constants from the file.
 
            call findnextinfoline(readunit, string)
            do ii = 1, cptempfit(nn)%nterm
 
                ! Read the constant from the string.
 
                read (string, *) cptempfit(nn)%constants(ii)
 
                ! Remove this value from the string if this is not the
                ! last constant to be read.
 
                if (ii < cptempfit(nn)%nterm) then
                    ios = index(string, " ")
                    if (ios > 0) then
                        string = string(ios:)
                        string = adjustl(string)
                        string = trim(string)
                    else
                        if (myid == 0) &
                             call terminate("readCpTempCurveFits", &
                             "Not enough constants on line; &
                             &Cp curve fit file not valid.")
                        call mpi_barrier(adflow_comm_world, ierr)
                    end if
                end if
            end do
 
        end do nranges
 
        ! Close the file
 
        close (unit=readunit)
        !
        !       Compute the constants eint0, such that the internal energy is
        !       a continous function of the temperature.
        !
        ! First for the first interval, such that at T = 0 Kelvin the
        ! energy is also zero.
 
        t1 = cptrange(0)
        cv0 = -one          ! cv/R = cp/R - 1.0
        e0 = -t1           ! e = integral of cv, not of cp.
 
        do ii = 1, cptempfit(1)%nterm
 
            ! Update cv0.
 
            t2 = t1**(cptempfit(1)%exponents(ii))
            cv0 = cv0 + cptempfit(1)%constants(ii) * t2
 
            ! Update e0, for which this contribution must be integrated.
            ! Take the exceptional case exponent is -1 into account.
 
            if (cptempfit(1)%exponents(ii) == -1_inttype) then
                e0 = e0 + cptempfit(1)%constants(ii) * log(t1)
            else
                t2 = t1 * t2
                e0 = e0 + cptempfit(1)%constants(ii) * t2 &
                     / (cptempfit(1)%exponents(ii) + 1)
            end if
 
        end do
 
        ! Set the value of the internal energy at the temperature T1.
        ! Cv is assumed to be constant in the temperature range 0 - T1.
        ! Idem for the internal enthalpy.
 
        cpeint(0) = cv0 * t1
        cphint(0) = cpeint(0) + t1
 
        ! Compute the integration constant for the energy.
 
        cptempfit(1)%eint0 = cpeint(0) - e0
 
        ! Loop over the other temperature ranges to compute their
        ! integration constant and the energy at the curve fit boundary.
 
        nranges2: do nn = 2, cpnparts
 
            ! Store nn-1, the previous temperature range, in mm.
 
            mm = nn - 1
 
            ! Store the temperature at the interface a bit easier.
 
            t1 = cptrange(mm)
 
            ! First compute the internal energy (scaled by r) from the
            ! previous range. Actually not the energy but the enthalpy is
            ! computed. This leads to the same integraton constant.
            ! Again check for exponent -1 when integrating.
 
            e0 = cptempfit(mm)%eint0
 
            do ii = 1, cptempfit(mm)%nterm
                if (cptempfit(mm)%exponents(ii) == -1_inttype) then
                    e0 = e0 + cptempfit(mm)%constants(ii) * log(t1)
                else
                    kk = cptempfit(mm)%exponents(ii) + 1
                    t2 = t1**kk
                    e0 = e0 + cptempfit(mm)%constants(ii) * t2 / kk
                end if
            end do
 
            ! Store the enthalpy and energy at the curve fit boundary.
            ! Remember that cp was integrated.
 
            cphint(mm) = e0
            cpeint(mm) = e0 - t1
 
            ! Substract the part coming from the integration of cp/r of
            ! the range nn.
 
            do ii = 1, cptempfit(nn)%nterm
                if (cptempfit(nn)%exponents(ii) == -1_inttype) then
                    e0 = e0 - cptempfit(nn)%constants(ii) * log(t1)
                else
                    kk = cptempfit(nn)%exponents(ii) + 1
                    t2 = t1**kk
                    e0 = e0 - cptempfit(nn)%constants(ii) * t2 / kk
                end if
            end do
 
            ! Store the integration constant for the range nn.
 
            cptempfit(nn)%eint0 = e0
 
        end do nranges2
 
        ! Compute the values of cv and the internal energy at the upper
        ! boundary of the curve fit. This is needed for the extrapolation
        ! of the energy if states occur with a higher temperature than
        ! the validness of the curve fits.
 
        ! First initialize these values.
 
        nn = cpnparts
        t1 = cptrange(nn)
        cvn = -one                       ! cv/R = cp/R - 1.0
        e0 = cptempfit(nn)%eint0 - t1  ! e = integral of cv, not of cp.
 
        do ii = 1, cptempfit(nn)%nterm
 
            ! Update cvn.
 
            t2 = t1**(cptempfit(nn)%exponents(ii))
            cvn = cvn + cptempfit(nn)%constants(ii) * t2
 
            ! Update e0, for which this contribution must be integrated.
            ! Take the exceptional case exponent is -1 into account.
 
            if (cptempfit(nn)%exponents(ii) == -1_inttype) then
                e0 = e0 + cptempfit(nn)%constants(ii) * log(t1)
            else
                e0 = e0 + cptempfit(nn)%constants(ii) * t2 * t1 &
                     / (cptempfit(nn)%exponents(ii) + 1)
            end if
 
        end do
 
        ! Store e0 correctly.
 
        cpeint(nn) = e0
        cphint(nn) = e0 + t1
 
        ! Compute the values of the integrands of cp/(R*T) at the lower
        ! and upper curve fit boundary. This is needed to compute the
        ! total pressure. This cannot be done with a single integration
        ! constant, because of the singularity at T = 0.
 
        nranges3: do nn = 1, cpnparts
 
            ! Store the temperatures of the lower and upper boundary a
            ! bit easier.
 
            t1 = cptrange(nn - 1)
            t2 = cptrange(nn)
 
            ! Initializes the integrands to zero.
 
            cptempfit(nn)%intCpovrT_1 = zero
            cptempfit(nn)%intCpovrT_2 = zero
 
            ! Loop over the number of terms of the curve fits and compute
            ! the integral cp/(r*t).
 
            do ii = 1, cptempfit(nn)%nterm
 
                ! Store the coefficient a bit easier in mm. As the integral
                ! of cp/(R*T) must be computed, this is also the exponent
                ! of the primitive function; except of course when the exponent
                ! is 0.
 
                mm = cptempfit(nn)%exponents(ii)
 
                ! Update the integrands if the temperature is larger than
                ! 0 kelvin. In case the boundary is 0 kelvin the value is not
                ! needed anyway.
 
                if (t1 > zero) then
                    if (mm == 0_inttype) then
                        cptempfit(nn)%intCpovrT_1 = cptempfit(nn)%intCpovrT_1 &
                                                    + cptempfit(nn)%constants(ii) * log(t1)
                    else
                        cptempfit(nn)%intCpovrT_1 = cptempfit(nn)%intCpovrT_1 &
                                                    + (cptempfit(nn)%constants(ii) * t1**mm) / mm
                    end if
                end if
 
                if (t2 > zero) then
                    if (mm == 0_inttype) then
                        cptempfit(nn)%intCpovrT_2 = cptempfit(nn)%intCpovrT_2 &
                                                    + cptempfit(nn)%constants(ii) * log(t2)
                    else
                        cptempfit(nn)%intCpovrT_2 = cptempfit(nn)%intCpovrT_2 &
                                                    + (cptempfit(nn)%constants(ii) * t2**mm) / mm
                    end if
                end if
 
            end do
 
        end do nranges3
 
    end subroutine readcptempcurvefits
 
    !      ==================================================================
 
    subroutine findnextinfoline(readUnit, string)
        !
        !       findNextInfoLine skips the comment lines in the given unit
        !       and finds the first line containing information.
        !
        use communication
        use utils, only: terminate
        implicit none
        !
        !      Subroutine arguments
        !
        integer, intent(in) :: readUnit
        character(len=512), intent(out) :: string
        !
        !      Local variables.
        !
        integer :: ios, ierr
 
        ! Loop to skip the comment lines.
 
        do
            read (unit=readunit, fmt="(a512)", iostat=ios) string
 
            ! Test if everything went okay.
 
            if (ios /= 0) then
                if (myid == 0) &
                    call terminate("findNextInfoLine", &
                                   "Unexpected end of Cp curve fit file")
                call mpi_barrier(adflow_comm_world, ierr)
            end if
 
            ! Get rid of the leading and trailing spaces in string.
 
            string = adjustl(string)
            string = trim(string)
 
            ! Check if this is the correct line. If so, exit
 
            if ((len_trim(string) > 0) .and. (string(:1) /= "#")) exit
        end do
 
    end subroutine findnextinfoline
 
    subroutine setequationparameters
        !
        !       setEquationParameters sets the number of variables in the
        !       governing equations, the number of turbulent variables, etc.
        !
        use constants
        use paramturb
        use turbcurvefits
        use flowvarrefstate, only: nw, nwf, nt1, nt2, nwt, viscous, &
                                   eddymodel, kpresent
        use inputphysics, only: equations, turbmodel, wallfunctions, rvfn
        implicit none
 
        ! Set the number of flow variables to 5, nt1 to 6. This is valid
        ! for all governing equations. Furthermore initialize viscous,
        ! kPresent and eddyModel to .False., which indicates an inviscid
        ! computation. For ns and rans this will be corrected.
 
        nwf = 5
        nt1 = 6
 
        viscous = .false.
        kpresent = .false.
        eddymodel = .false.
 
        ! Determine the set of governing equations to solve for and set
        ! the parameters accordingly.
 
        select case (equations)
        case (eulerequations)
            nw = 5
            nt2 = 5
 
            !===============================================================
 
        case (nsequations)
            nw = 5
            nt2 = 5
 
            viscous = .true.
 
            !===============================================================
 
        case (ransequations)
 
            viscous = .true.
 
            select case (turbmodel)
 
            case (spalartallmaras)
                nw = 6
                nt2 = 6
 
                eddymodel = .true.
                if (wallfunctions) call initcurvefitdatasa
 
                !===========================================================
 
            case (spalartallmarasedwards)
                nw = 6
                nt2 = 6
 
                eddymodel = .true.
                if (wallfunctions) call initcurvefitdatasae
 
                !===========================================================
 
            case (komegawilcox)
                nw = 7
                nt2 = 7
 
                kpresent = .true.
                eddymodel = .true.
                if (wallfunctions) call initcurvefitdatakw
 
                !===========================================================
 
            case (komegamodified)
                nw = 7
                nt2 = 7
 
                kpresent = .true.
                eddymodel = .true.
                if (wallfunctions) call initcurvefitdatakwmod
 
                !===========================================================
 
            case (mentersst)
                nw = 7
                nt2 = 7
 
                kpresent = .true.
                eddymodel = .true.
                if (wallfunctions) call initcurvefitdatasst
 
                !===========================================================
 
            case (ktau)
                nw = 7
                nt2 = 7
 
                kpresent = .true.
                eddymodel = .true.
                if (wallfunctions) call initcurvefitdataktau
 
                !===========================================================
 
            case (v2f)
                nw = 9
                nt2 = 9
 
                rvflimitk = 1.e-25_realtype
                rvflimite = 1.e-25_realtype
 
                if (rvfn == 6) then
                    rvfcmu = rvfn6cmu
                    rvfcl = rvfn6cl
                else
                    rvfcmu = rvfn1cmu
                    rvfcl = rvfn1cl
                end if
 
                kpresent = .true.
                eddymodel = .true.
                if (wallfunctions) call initcurvefitdatavf
 
            end select
 
        end select
 
        ! Determine the number of turbulent variables.
 
        nwt = nw - nwf
 
    end subroutine setequationparameters
    subroutine setstagecoeffexplicitrk
        !
        !       setStageCoeffExplicitRK determines the coefficients of the
        !       stages for the explicit Runge Kutta time integration schemes
        !       for unsteady problems.
        !
        use constants
        use inputunsteady, only: timeaccuracy, betarkunsteady, &
                                 gammarkunsteady, nrkstagesunsteady
        use utils, only: terminate
        implicit none
        !
        !      Local variables.
        !
        integer :: ierr
 
        ! Determine the number of Runge Kutta stages as a function of
        ! the accuracy.
 
        select case (timeaccuracy)
        case (firstorder)
            nrkstagesunsteady = 1
 
        case (secondorder)
            nrkstagesunsteady = 2
 
        case (thirdorder)
            nrkstagesunsteady = 3
 
        case default
            call terminate("setStageCoeffExplicitRK", &
                           "No higher order stuff yet")
        end select
 
        ! Allocate and determine betaRKUnsteady and gammaRKUnsteady.
 
        allocate (betarkunsteady(nrkstagesunsteady, nrkstagesunsteady), &
                  gammarkunsteady(nrkstagesunsteady), stat=ierr)
        if (ierr /= 0) &
             call terminate("setStageCoeffExplicitRK", &
             "Memory allocation failure for betaRKUnsteady &
             &and gammaRKUnsteady.")
 
        betarkunsteady = zero
 
        select case (timeaccuracy)
        case (firstorder)
 
            ! Just the forward Euler time integration scheme.
 
            betarkunsteady(1, 1) = 1.0_realtype
            gammarkunsteady(1) = 0.0_realtype
 
            !==============================================================
 
        case (secondorder)
 
            ! The TVD Runge Kutta scheme which allows for the maximum
            ! CFL number (1.0).
 
            betarkunsteady(1, 1) = 1.0_realtype
            betarkunsteady(2, 1) = -0.5_realtype
            betarkunsteady(2, 2) = 0.5_realtype
 
            gammarkunsteady(1) = 0.0_realtype
            gammarkunsteady(2) = 1.0_realtype
 
            !==============================================================
 
        case (thirdorder)
 
            ! Low storage (although not exploited in this implemetation)
            ! 3 stage scheme of Le and Moin.
 
            betarkunsteady(1, 1) = 8.0_realtype / 15.0_realtype
            betarkunsteady(2, 1) = -17.0_realtype / 60.0_realtype
            betarkunsteady(2, 2) = 5.0_realtype / 12.0_realtype
            betarkunsteady(3, 2) = -5.0_realtype / 12.0_realtype
            betarkunsteady(3, 3) = 3.0_realtype / 4.0_realtype
 
            gammarkunsteady(1) = 0.0_realtype
            gammarkunsteady(2) = 8.0_realtype / 15.0_realtype
            gammarkunsteady(3) = 2.0_realtype / 3.0_realtype
 
            ! The TVD Runge Kutta scheme which allows for the maximum
            ! CFL number (1.0).
 
            ! betaRKUnsteady(1,1) =  1.0_realType
            ! betaRKUnsteady(2,1) = -3.0_realType/ 4.0_realType
            ! betaRKUnsteady(2,2) =  1.0_realType/ 4.0_realType
            ! betaRKUnsteady(3,1) = -1.0_realType/12.0_realType
            ! betaRKUnsteady(3,2) = -1.0_realType/12.0_realType
            ! betaRKUnsteady(3,3) =  2.0_realType/ 3.0_realType
 
            ! gammaRKUnsteady(1)  = 0.0_realType
            ! gammaRKUnsteady(2)  = 1.0_realType
            ! gammaRKUnsteady(3)  = 0.5_realType
 
            !==============================================================
 
        case default
            call terminate("setStageCoeffExplicitRK", &
                           "No higher order stuff yet")
        end select
 
    end subroutine setstagecoeffexplicitrk
    subroutine surfacevariables(variables)
        !
        !       surfaceVariables extracts from the given string the surface
        !       variables to be written to the solution file.
        !
        use constants
        use extraoutput
        use utils, only: converttolowercase, terminate
        implicit none
        !
        !      Subroutine arguments.
        !
        character(len=*), intent(inout) :: variables
        !
        !      Local variables.
        !
        integer :: nVarSpecified, pos
 
        character(len=15) :: keyword
        character(len=maxStringLen) :: errorMessage
 
        ! Convert the string variables to lower case.
 
        call converttolowercase(variables)
 
        ! Initialize all the surface output variables to .false.
 
        surfwriterho = .false.
        surfwritep = .false.
        surfwritetemp = .false.
        surfwritevx = .false.
        surfwritevy = .false.
        surfwritevz = .false.
        surfwritervx = .false.
        surfwritervy = .false.
        surfwritervz = .false.
 
        surfwritecp = .false.
        surfwriteptotloss = .false.
        surfwritemach = .false.
        surfwritermach = .false.
 
        surfwritecf = .false.
        surfwritech = .false.
        surfwriteyplus = .false.
        surfwritecfx = .false.
        surfwritecfy = .false.
        surfwritecfz = .false.
        surfwriteforceindragdir = .false.
        surfwriteforceinliftdir = .false.
 
        surfwriteblank = .false.
        surfwritesepsensor = .false.
        surfwritesepsensorks = .false.
        surfwritesepsensorksarea = .false.
        surfwritecavitation = .false.
        surfwriteaxismoment = .false.
        surfwritegc = .false.
 
        ! Initialize nVarSpecified to 0. This serves as a test
        ! later on.
 
        nvarspecified = 0
 
        ! Loop to extract the info from the string variables.
 
        do
            ! Condition to exit the loop.
 
            if (len_trim(variables) == 0) exit
 
            ! Locate the first occurance of the _ in the string and
            ! determine the string keyword.
 
            pos = index(variables, "_")
            if (pos == 0) then
                keyword = variables
                variables = ""
            else
                keyword = variables(:pos - 1)
                variables = variables(pos + 1:)
            end if
 
            ! Check the keyword.
 
            select case (keyword)
            case ("")
                ! Multiple occurence of "_". Just ignore it.
 
            case ("rho")
                surfwriterho = .true.
                nvarspecified = nvarspecified + 1
 
            case ("p")
                surfwritep = .true.
                nvarspecified = nvarspecified + 1
 
            case ("temp")
                surfwritetemp = .true.
                nvarspecified = nvarspecified + 1
 
            case ("vx")
                surfwritevx = .true.
                nvarspecified = nvarspecified + 1
 
            case ("vy")
                surfwritevy = .true.
                nvarspecified = nvarspecified + 1
 
            case ("vz")
                surfwritevz = .true.
                nvarspecified = nvarspecified + 1
 
            case ("rvx")
                surfwritervx = .true.
                nvarspecified = nvarspecified + 1
 
            case ("rvy")
                surfwritervy = .true.
                nvarspecified = nvarspecified + 1
 
            case ("rvz")
                surfwritervz = .true.
                nvarspecified = nvarspecified + 1
 
            case ("cp")
                surfwritecp = .true.
                nvarspecified = nvarspecified + 1
 
            case ("ptloss")
                surfwriteptotloss = .true.
                nvarspecified = nvarspecified + 1
 
            case ("mach")
                surfwritemach = .true.
                nvarspecified = nvarspecified + 1
 
            case ("rmach")
                surfwritermach = .true.
                nvarspecified = nvarspecified + 1
 
            case ("cf")
                surfwritecf = .true.
                nvarspecified = nvarspecified + 1
 
            case ("ch")
                surfwritech = .true.
                nvarspecified = nvarspecified + 1
 
            case ("yplus")
                surfwriteyplus = .true.
                nvarspecified = nvarspecified + 1
 
            case ("cfx")
                surfwritecfx = .true.
                nvarspecified = nvarspecified + 1
 
            case ("cfy")
                surfwritecfy = .true.
                nvarspecified = nvarspecified + 1
 
            case ("cfz")
                surfwritecfz = .true.
                nvarspecified = nvarspecified + 1
 
            case ("forceindragdir")
                surfwriteforceindragdir = .true.
                nvarspecified = nvarspecified + 1
 
            case ("forceinliftdir")
                surfwriteforceinliftdir = .true.
                nvarspecified = nvarspecified + 1
 
            case ("blank")
                surfwriteblank = .true.
                nvarspecified = nvarspecified + 1
 
            case ("sepsensor")
                surfwritesepsensor = .true.
                nvarspecified = nvarspecified + 1
 
            case ("sepsensorks")
                surfwritesepsensorks = .true.
                nvarspecified = nvarspecified + 1
 
            case ("sepsensorksarea")
                surfwritesepsensorksarea = .true.
                nvarspecified = nvarspecified + 1
 
            case ("cavitation")
                surfwritecavitation = .true.
                nvarspecified = nvarspecified + 1
 
            case ("axismoment")
                surfwriteaxismoment = .true.
                nvarspecified = nvarspecified + 1
 
            case ("gc")
                surfwritegc = .true.
                nvarspecified = nvarspecified + 1
 
            case default
                pos = len_trim(keyword)
                write (errormessage, "(3a)") "Unknown surface output &
                     &variable, ", trim(keyword), &
                     ", specified"
                call terminate("surfaceVariables", errormessage)
 
            end select
 
        end do
 
        ! Set surfaceOutSpecified to .true. if variables were specified.
        ! If not, later on the defaults will be set.
 
        if (nvarspecified > 0) surfaceoutspecified = .true.
 
    end subroutine surfacevariables
 
    subroutine volumevariables(variables)
        !
        !       volumeVariables extracts from the given string the extra
        !       volume variables to be written to the solution file.
        !
        use constants
        use extraoutput
        use utils, only: converttolowercase, terminate
        implicit none
        !
        !      Subroutine arguments.
        !
        character(len=*), intent(inout) :: variables
        !
        !      Local variables.
        !
        integer :: nVarSpecified, pos
 
        character(len=15) :: keyword
        character(len=maxStringLen) :: errorMessage
 
        ! Convert the string variables to lower case.
 
        call converttolowercase(variables)
 
        ! Initialize all the volume output variables to .False.
 
        volwritemx = .false.
        volwritemy = .false.
        volwritemz = .false.
        volwriterhoe = .false.
        volwritetemp = .false.
        volwritevort = .false.
        volwritevortx = .false.
        volwritevorty = .false.
        volwritevortz = .false.
 
        volwritecp = .false.
        volwritemach = .false.
        volwritemachturb = .false.
        volwriteptotloss = .false.
 
        volwriteeddyvis = .false.
        volwriteratioeddyvis = .false.
        volwritedist = .false.
 
        volwriteresrho = .false.
        volwriteresmom = .false.
        volwriteresrhoe = .false.
        volwriteresturb = .false.
 
        volwriteshock = .false.
        volwritefilteredshock = .false.
 
        volwriteblank = .false.
        volwritegc = .false.
        volwritestatus = .false.
        volwriteintermittency = .false.
 
        ! Initialize nVarSpecified to 0. This serves as a test
        ! later on.
 
        nvarspecified = 0
 
        ! Loop to extract the info from the string variables.
 
        do
            ! Condition to exit the loop.
 
            if (len_trim(variables) == 0) exit
 
            ! Locate the first occurance of the _ in the string and
            ! determine the string keyword.
 
            pos = index(variables, "_")
            if (pos == 0) then
                keyword = variables
                variables = ""
            else
                keyword = variables(:pos - 1)
                variables = variables(pos + 1:)
            end if
 
            ! Check the keyword.
 
            select case (keyword)
            case ("")
                ! Multiple occurence of "_". Just ignore it.
 
            case ("mx")
                volwritemx = .true.
                nvarspecified = nvarspecified + 1
 
            case ("my")
                volwritemy = .true.
                nvarspecified = nvarspecified + 1
 
            case ("mz")
                volwritemz = .true.
                nvarspecified = nvarspecified + 1
 
            case ("rvx")
                volwritervx = .true.
                nvarspecified = nvarspecified + 1
 
            case ("rvy")
                volwritervy = .true.
                nvarspecified = nvarspecified + 1
 
            case ("rvz")
                volwritervz = .true.
                nvarspecified = nvarspecified + 1
 
            case ("rhoe")
                volwriterhoe = .true.
                nvarspecified = nvarspecified + 1
 
            case ("temp")
                volwritetemp = .true.
                nvarspecified = nvarspecified + 1
 
            case ("vort")
                volwritevort = .true.
                nvarspecified = nvarspecified + 1
 
            case ("vortx")
                volwritevortx = .true.
                nvarspecified = nvarspecified + 1
 
            case ("vorty")
                volwritevorty = .true.
                nvarspecified = nvarspecified + 1
 
            case ("vortz")
                volwritevortz = .true.
                nvarspecified = nvarspecified + 1
 
            case ("cp")
                volwritecp = .true.
                nvarspecified = nvarspecified + 1
 
            case ("mach")
                volwritemach = .true.
                nvarspecified = nvarspecified + 1
 
            case ("rmach")
                volwritermach = .true.
                nvarspecified = nvarspecified + 1
 
            case ("macht")
                volwritemachturb = .true.
                nvarspecified = nvarspecified + 1
 
            case ("ptloss")
                volwriteptotloss = .true.
                nvarspecified = nvarspecified + 1
 
            case ("eddy")
                volwriteeddyvis = .true.
                nvarspecified = nvarspecified + 1
 
            case ("eddyratio")
                volwriteratioeddyvis = .true.
                nvarspecified = nvarspecified + 1
 
            case ("dist")
                volwritedist = .true.
                nvarspecified = nvarspecified + 1
 
            case ("resrho")
                volwriteresrho = .true.
                nvarspecified = nvarspecified + 1
 
            case ("resmom")
                volwriteresmom = .true.
                nvarspecified = nvarspecified + 1
 
            case ("resrhoe")
                volwriteresrhoe = .true.
                nvarspecified = nvarspecified + 1
 
            case ("resturb")
                volwriteresturb = .true.
                nvarspecified = nvarspecified + 1
 
            case ("shock")
                volwriteshock = .true.
                nvarspecified = nvarspecified + 1
 
            case ("filteredshock")
                volwritefilteredshock = .true.
                nvarspecified = nvarspecified + 1
 
            case ("blank")
                volwriteblank = .true.
                nvarspecified = nvarspecified + 1
 
            case ("gc")
                volwritegc = .true.
                nvarspecified = nvarspecified + 1
 
            case ("status")
                volwritestatus = .true.
                nvarspecified = nvarspecified + 1
 
            case ("intermittency")
                volwriteintermittency = .true.
                nvarspecified = nvarspecified + 1
 
            case default
                pos = len_trim(keyword)
                write (errormessage, "(3a)") "Unknown extra volume output &
                     &variable, ", trim(keyword), &
                     ", specified"
                call terminate("volumeVariables", errormessage)
 
            end select
 
        end do
 
        ! Set volumeOutSpecified to .true. if variables were specified.
        ! If not, later on the defaults will be set.
 
        if (nvarspecified > 0) volumeoutspecified = .true.
 
    end subroutine volumevariables
 
    subroutine checkinputparam
        !
        !       checkInputParam checks if all necessary data has been
        !       specified. If some key data is missing an error message will
        !       be printed and the program will exit. Key data depends on the
        !       case to be solved. E.g. for the Navier Stokes equations it is
        !       necessary to specify the Reynolds number, but for Euler this
        !       can be omitted.
        !       Furthermore warnings are printed in case parameters have been
        !       specified that are ignored, e.g. Mach number for internal flow
        !       computations.
        !       Note that only processor 0 prints warning and error messages,
        !       such that the output does not become messy.
        !
        use constants
        !  --------- Bare imports...too many to list -------
        use inputdiscretization
        use inputio
        use inputiteration
        use inputmotion
        use inputoverset
        use inputparallel
        use inputphysics
        use inputtimespectral
        use inputunsteady
        use inputadjoint
        use inputtsstabderiv
        !  ------------------------------------------------
        use communication, only: myid, adflow_comm_world
        use iteration, only: coeftime, coeftimeale, coefmeshale, &
                             oldsolwritten, nalemeshes, nalesteps, noldlevels
        use monitor, only: ntimestepsrestart
        use utils, only: terminate
        implicit none
        !
        !      Local variables
        !
        integer :: ierr
 
        integer(kind=intType) :: nn, oldSolWrittenSize
 
        real(kind=realtype) :: veclength, dot
 
        logical :: gridPrecisionWarning, solPrecisionWarning
 
        !       Discretization parameters. Check if the key parameters have
        !       been specified and set some coarse grid parameters in case
        !       these have not been specified.
        !
        if (spacediscr == none) then
            if (myid == 0) &
                call terminate("checkInputParam", &
                               "Discretization scheme not specified")
            call mpi_barrier(adflow_comm_world, ierr)
        end if
 
        if (spacediscrcoarse == none) spacediscrcoarse = spacediscr
 
        if (riemanncoarse == none) riemanncoarse = riemann
 
        ! Set dirScaling to .false. if a scheme other than scalar
        ! dissipation is used.
 
        if (spacediscr /= dissscalar) dirscaling = .false.
 
        ! Determine whether or not the spectral radIi are needed for
        ! the flux computations.
 
        radiineededfine = .false.
        if (spacediscr == dissscalar) radiineededfine = .true.
 
        radiineededcoarse = .false.
        if (spacediscrcoarse == dissscalar) radiineededcoarse = .true.
        !
        !       IO parameters. Check if the grid file has been specified
        !       Possibly correct the
        !       value of restart. Note that restart got the default value of
        !       .true. in case no restart file has been specified it is now
        !       set to false. Set the names of the solution files if not
        !       specified and check if a cp curve fit file has been specified
        !       if curve fits must be used.
        !       If the code has been compiled without cgns check that the file
        !       format is not cgns.
        !       Overwrite storeConvInnerIter to .true. if this is not an
        !       unsteady computation.
        !
        if (gridfile == "") then
            if (myid == 0) &
                call terminate("checkInputParam", "Grid file not specified")
            call mpi_barrier(adflow_comm_world, ierr)
        end if
 
        if (newgridfile == "") then
            newgridfile = "NewGrid.cgns"
        end if
 
        if (solfile == "") then
            solfile = "SolADflow.cgns"
        end if
 
        if (surfacesolfile == "") &
            surfacesolfile = trim(solfile)//"Surface"
 
        if (cpmodel == cptempcurvefits .and. cpfile == "") then
            if (myid == 0) &
                call terminate("checkInputParam", &
                               "Cp curve fit file not specified")
            call mpi_barrier(adflow_comm_world, ierr)
        end if
 
#ifdef USE_NO_CGNS
 
        if (fileformatread == cgnsformat .or. &
            fileformatwrite == cgnsformat) then
            if (myid == 0) &
                call terminate("checkInputParam", &
                               "cgns support disabled during compile time")
            call mpi_barrier(adflow_comm_world, ierr)
        end if
 
#endif
 
        if (equationmode == unsteady) then
            if (timeintegrationscheme == explicitrk) &
                storeconvinneriter = .false.
        end if
        !
        !       Iteration parameters. Check if the key parameters have specified
        !       been and set some coarse grid parameters in case these
        !       have not been specified.
        !
        if (equationmode == unsteady .and. &
            timeintegrationscheme == explicitrk) then
            smoother = none
        else
            if (smoother == none) then
                if (myid == 0) &
                    call terminate("checkInputParam", "Smoother not specified")
                call mpi_barrier(adflow_comm_world, ierr)
            end if
 
            if (ncycles < 0) then
                if (myid == 0) &
                     call terminate("checkInputParam", &
                     "Number of multigrid cycles not or wrongly &
                     &specified")
                call mpi_barrier(adflow_comm_world, ierr)
            end if
 
            if (cfl < zero) then
                if (myid == 0) &
                    call terminate("checkInputParam", &
                                   "cfl number not or wrongly specified")
                call mpi_barrier(adflow_comm_world, ierr)
            end if
 
            if (l2conv <= zero .or. l2conv >= one) then
                if (myid == 0) &
                     call terminate("checkInputParam", &
                     "Relative L2 norm for convergence must be a &
                     & number between 0 and 1.")
                call mpi_barrier(adflow_comm_world, ierr)
            end if
 
            if (l2convcoarse <= zero .or. l2convcoarse >= one) then
                if (myid == 0) &
                     call terminate("checkInputParam", &
                     "Relative L2 norm for convergence coarse grid &
                     &must be a number between 0 and 1.")
                call mpi_barrier(adflow_comm_world, ierr)
            end if
        end if
        !
        !       Grid motion parameters. These can only be specified for an
        !       external flow problem.
        !
        if (flowtype == internalflow .and. gridmotionspecified) then
            if (myid == 0) &
                 call terminate("checkInputParam", &
                 "Grid motion specified for an internal flow; &
                 &this is not possible")
            call mpi_barrier(adflow_comm_world, ierr)
        end if
        !
        !       Physics parameters. Check if the key parameters have been
        !       specified and set the unit vector for the free-stream velocity.
        !
        if (equations == none) then
            if (myid == 0) &
                call terminate("checkInputParam", "Equations not specified")
            call mpi_barrier(adflow_comm_world, ierr)
        end if
 
        if (equationmode == none) then
            if (myid == 0) &
                call terminate("checkInputParam", "Mode not specified")
            call mpi_barrier(adflow_comm_world, ierr)
        end if
 
        if (flowtype == none) then
            if (myid == 0) &
                call terminate("checkInputParam", "Flow type not specified")
            call mpi_barrier(adflow_comm_world, ierr)
        end if
 
        if (mach < zero .and. flowtype == externalflow) then
            if (myid == 0) &
                call terminate("checkInputParam", &
                               "Mach not or wrongly specified")
            call mpi_barrier(adflow_comm_world, ierr)
        end if
 
        if (equations == ransequations .and. turbmodel == none) then
            if (myid == 0) &
                call terminate("checkInputParam", &
                               "Turbulence model not specified")
            call mpi_barrier(adflow_comm_world, ierr)
        end if
 
        ! Create a unit vector for the free stream velocity. It is checked
        ! if the vector specified is a valid one. If not processor 0 prints
        ! an error message. Only for external flows.
 
        if (flowtype == externalflow) then
            veclength = sqrt(veldirfreestream(1) * veldirfreestream(1) &
                             + veldirfreestream(2) * veldirfreestream(2) &
                             + veldirfreestream(3) * veldirfreestream(3))
            if (veclength < eps) then
                if (myid == 0) &
                     call terminate("checkInputParam", &
                     "Free stream velocity direction wrongly &
                     &specified")
                call mpi_barrier(adflow_comm_world, ierr)
            end if
 
            veclength = one / veclength
            veldirfreestream(1) = veldirfreestream(1) * veclength
            veldirfreestream(2) = veldirfreestream(2) * veclength
            veldirfreestream(3) = veldirfreestream(3) * veclength
        else
            ! Internal flow; simply reset the velocity direction. The value
            ! will be determined later from the inflow boundary conditions.
 
            veldirfreestream(1) = one
            veldirfreestream(2) = zero
            veldirfreestream(3) = zero
        end if
 
        ! Set the drag direction to the velocity direction.
 
        dragdirection = veldirfreestream
 
        ! Check the lift direction if it was specified for an external
        ! flow. Otherwise set the default direction.
 
        if (liftdirspecified .and. flowtype == externalflow) then
 
            ! Create a unit vector. Perform the same check as for
            ! for the free stream velocity direction.
 
            veclength = sqrt(liftdirection(1) * liftdirection(1) &
                             + liftdirection(2) * liftdirection(2) &
                             + liftdirection(3) * liftdirection(3))
            if (veclength < eps) then
                if (myid == 0) &
                    call terminate("checkInputParam", &
                                   "Lift direction wrongly specified")
                call mpi_barrier(adflow_comm_world, ierr)
            end if
 
            veclength = one / veclength
            liftdirection(1) = liftdirection(1) * veclength
            liftdirection(2) = liftdirection(2) * veclength
            liftdirection(3) = liftdirection(3) * veclength
 
            ! Check the orthogonality with the drag direction.
 
            dot = liftdirection(1) * dragdirection(1) &
                  + liftdirection(2) * dragdirection(2) &
                  + liftdirection(3) * dragdirection(3)
 
            if (abs(dot) > 1.e-3_realtype) then
                if (myid == 0) &
                     call terminate("checkInputParam", &
                     "Lift direction not orthogonal to &
                     &free-stream")
                call mpi_barrier(adflow_comm_world, ierr)
            end if
 
        else
 
            ! Lift direction not specified. Set the default direction.
            ! It will have a zero component in the y-direction and a positive
            ! one in the z-direction.
 
            liftdirection(1) = -dragdirection(3)
            liftdirection(2) = zero
            liftdirection(3) = dragdirection(1)
 
            if (liftdirection(3) < zero) then
                liftdirection(1) = -liftdirection(1)
                liftdirection(3) = -liftdirection(3)
            end if
        end if
 
        ! Set the Mach number for the coefficients equal to the Mach
        ! number if it was not specified. For internal flow field this
        ! will again be changed in initFlo.
 
        if (machcoef < zero) machcoef = mach
        !
        !       Time spectral parameters. They only need to be specified for a
        !       time spectral computation.
        !
        testspectral: if (equationmode == timespectral) then
 
            ! Check if the number of time intervals was specified.
 
            if (ntimeintervalsspectral < 0) then
                if (myid == 0) &
                     call terminate("checkInputParam", &
                     "Number time intervals spectral not or &
                     &wrongly specified")
                call mpi_barrier(adflow_comm_world, ierr)
            end if
 
            ! If an unsteady restart solution file must be written, check
            ! if the corresponding time step has been specified.
 
            if (writeunsteadyrestartspectral) then
                if (dtunsteadyrestartspectral <= zero) then
                    if (myid == 0) &
                         call terminate("checkInputParam", &
                         "Time step (in sec) for unsteady restart &
                         &not or wrongly specified.")
                    call mpi_barrier(adflow_comm_world, ierr)
                end if
            end if
 
            ! If solution files (for postprocessing) must be written,
            ! check if the number has been specified.
 
            if (writeunsteadyvolspectral .or. &
                writeunsteadysurfspectral) then
                if (nunsteadysolspectral <= 0) then
                    if (myid == 0) &
                         call terminate("checkInputParam", &
                         "Number of unsteady solution files &
                         &not or wrongly specified.")
                    call mpi_barrier(adflow_comm_world, ierr)
                end if
            end if
 
            else testspectral
 
            ! No spectral method. Set nTimeIntervalsSpectral to 1.
 
            ntimeintervalsspectral = 1
 
        end if testspectral
        !
        !       Unsteady parameters. They only need to be specified for an
        !       unsteady computation.
        !
        testunsteady: if (equationmode == unsteady) then
 
            ! Physical time step parameters.
 
            if (ntimestepsfine < 0) then
                if (myid == 0) &
                     call terminate("checkInputParam", &
                     "Number of unsteady time steps fine grid &
                     &not or wrongly specified")
                call mpi_barrier(adflow_comm_world, ierr)
            end if
 
            if (ntimestepscoarse < 0) ntimestepscoarse = ntimestepsfine
 
            if (deltat < 0) then
                if (myid == 0) &
                     call terminate("checkInputParam", &
                     "Unsteady time step (in sec) &
                     &not or wrongly specified")
                call mpi_barrier(adflow_comm_world, ierr)
            end if
 
            ! Check if the rigid body rotation parameters are consistent.
            ! The polynomial rotation coefficients.
 
            if (degreepolxrot >= 0 .and. &
                .not. allocated(coefpolxrot)) then
                if (myid == 0) &
                     call terminate("checkInputParam", &
                     "Polynomial coefficients x-rotation &
                     &not specified")
                call mpi_barrier(adflow_comm_world, ierr)
            end if
 
            if (degreepolyrot >= 0 .and. &
                .not. allocated(coefpolyrot)) then
                if (myid == 0) &
                     call terminate("checkInputParam", &
                     "Polynomial coefficients y-rotation &
                     &not specified")
                call mpi_barrier(adflow_comm_world, ierr)
            end if
 
            if (degreepolzrot >= 0 .and. &
                .not. allocated(coefpolzrot)) then
                if (myid == 0) &
                     call terminate("checkInputParam", &
                     "Polynomial coefficients z-rotation &
                     &not specified")
                call mpi_barrier(adflow_comm_world, ierr)
            end if
 
            ! The fourier rotation coefficients.
 
            if (degreefourxrot >= 0 .and. &
                .not. allocated(coscoeffourxrot)) then
                if (myid == 0) &
                     call terminate("checkInputParam", &
                     "Fourier cosine coefficients x-rotation &
                     &not specified")
                call mpi_barrier(adflow_comm_world, ierr)
            end if
 
            if (degreefourxrot >= 1 .and. &
                .not. allocated(sincoeffourxrot)) then
                if (myid == 0) &
                     call terminate("checkInputParam", &
                     "Fourier sine coefficients x-rotation &
                     &not specified")
                call mpi_barrier(adflow_comm_world, ierr)
            end if
 
            if (degreefouryrot >= 0 .and. &
                .not. allocated(coscoeffouryrot)) then
                if (myid == 0) &
                     call terminate("checkInputParam", &
                     "Fourier cosine coefficients y-rotation &
                     &not specified")
                call mpi_barrier(adflow_comm_world, ierr)
            end if
 
            if (degreefouryrot >= 1 .and. &
                .not. allocated(sincoeffouryrot)) then
                if (myid == 0) &
                     call terminate("checkInputParam", &
                     "Fourier sine coefficients y-rotation &
                     &not specified")
                call mpi_barrier(adflow_comm_world, ierr)
            end if
 
            if (degreefourzrot >= 0 .and. &
                .not. allocated(coscoeffourzrot)) then
                if (myid == 0) &
                     call terminate("checkInputParam", &
                     "Fourier cosine coefficients z-rotation &
                     &not specified")
                call mpi_barrier(adflow_comm_world, ierr)
            end if
 
            if (degreefourzrot >= 1 .and. &
                .not. allocated(sincoeffourzrot)) then
                if (myid == 0) &
                     call terminate("checkInputParam", &
                     "Fourier sine coefficients z-rotation &
                     &not specified")
                call mpi_barrier(adflow_comm_world, ierr)
            end if
 
        end if testunsteady
        !
        !                             Warning messages.
        !
        ! Check for an invisid problem if the Reynolds number is specified.
        ! If so, print a Warning that this info is ignored.
 
        if (myid == 0 .and. equations == eulerequations .and. &
            reynolds > zero) then
 
            print "(a)", "#"
            print "(a)", "#                      Warning"
            print "(a)", "# Reynolds number specified for the Euler &
                 &equations."
            print "(a)", "# This information is ignored."
            print "(a)", "#"
 
        end if
 
        ! Check if the Mach and Reynolds number are specified for an
        ! internal flow problem. If so, print a Warning message that this
        ! info is ignored.
 
        if (flowtype == internalflow) then
 
            ! Check whether a viscous or an inviscid problem is to be solved.
            ! For an inviscid problem you do not want to mention that the
            ! Reynolds number is ignored, because this has already been
            ! taken care of.
 
            if ((equations == nsequations .or. &
                 equations == ransequations) .and. &
                mach > zero .and. reynolds > zero) then
 
                ! Viscous problem, where both the Mach and Reynolds were
                ! specified. Processor 0 prints the Warning.
 
                if (myid == 0) then
                    print "(a)", "#"
                    print "(a)", "#                      Warning"
                    print "(a)", "# Mach and Reynolds number specified &
                         &for an internal flow problem."
                    print "(a)", "# This information is ignored."
                    print "(a)", "#"
                end if
 
            else if (mach > zero) then
 
                ! The Mach number has been specified. Processor 0 prints
                ! a Warning.
 
                if (myid == 0) then
                    print "(a)", "#"
                    print "(a)", "#                      Warning"
                    print "(a)", "# Mach number specified for an internal &
                         &flow problem."
                    print "(a)", "# This information is ignored."
                    print "(a)", "#"
                end if
 
            end if
 
        end if
 
        ! For a steady computation possible specified rigid body
        ! rotation info is ignored. Processor 0 will print the Warning.
 
        if (degreepolxrot >= 0 .or. degreepolyrot >= 0 .or. &
            degreepolzrot >= 0 .or. degreefourxrot >= 0 .or. &
            degreefouryrot >= 0 .or. degreefourzrot >= 0) then
 
            if (equationmode == steady .and. myid == 0) then
                print "(a)", "#"
                print "(a)", "#                      Warning"
                print "(a)", "# Rigid body rotation info specified for &
                     &a steady computation."
                print "(a)", "# This information is ignored."
                print "(a)", "#"
            end if
        end if
 
        ! Print warning messages if the precision to be written
        ! is larger than the precision used in the computation.
 
        gridprecisionwarning = .false.
        solprecisionwarning = .false.
 
#ifdef USE_SINGLE_PRECISION
        if (precisiongrid == precisiondouble) gridprecisionwarning = .true.
        if (precisionsol == precisiondouble) solprecisionwarning = .true.
#endif
 
        if (gridprecisionwarning .and. myid == 0) then
            print "(a)", "#"
            print "(a)", "#                      Warning"
            print "(a)", "# Precision of the grid file to write is &
                 &bigger than used in the computation."
            print "(a)", "# This does not make sense and is a waste &
                 &of disk space"
            print "(a)", "#"
        end if
 
        if (solprecisionwarning .and. myid == 0) then
            print "(a)", "#"
            print "(a)", "#                      Warning"
            print "(a)", "# Precision of the solution file to write is &
                 &bigger than used in the computation."
            print "(a)", "# This does not make sense and is a waste &
                 &of disk space"
            print "(a)", "#"
        end if
        !
        !       Wall functions can only be used if the RANS equations are to
        !       be solved. If no wall functions are used the wall offset is
        !       set to zero.
        !
        if (equations /= ransequations) wallfunctions = .false.
        if (.not. wallfunctions) walloffset = zero
        !
        !       Check whether or not the wall distance is needed for the
        !       turbulence model.
        !
        if (equations == ransequations) then
 
            ! RANS simulation. Determine if the turbulence model is
            ! wall distance free. Note that updateWallDistanceUnsteady is
            ! NOT overruled, because this is just the case for which this
            ! parameter was intended.
 
            select case (turbmodel)
            case (komegawilcox, komegamodified, ktau)
 
                ! Wall distance free turbulence models.
 
                walldistanceneeded = .false.
 
                !=============================================================
 
            case default
 
                ! The turbulence model needs the wall distance
 
                walldistanceneeded = .true.
 
            end select
 
        else
 
            ! Laminar or inviscid computation. Simply initialize the
            ! logicals for the wall distance to .false.
 
            walldistanceneeded = .false.
            updatewalldistanceunsteady = .false.
 
        end if
        !
        !       Parallelization parameters. Set the minimum load imbalance to
        !       3 percent to avoid any problems.
        !
        loadimbalance = max(loadimbalance, 0.03_realtype)
        !
        !       Some default parameters, which depend on other parameters.
        !       Only if these have not been specified of course.
        !
        if (nsgstartup < 0) nsgstartup = 0
        if (ncyclescoarse < 0) ncyclescoarse = ncycles
        if (cflcoarse < zero) cflcoarse = cfl
        if (betaturb < zero) betaturb = alfaturb
 
        if (turbrelax == turbrelaxnotdefined) then
            turbrelax = turbrelaximplicit
            if (turbmodel == v2f) turbrelax = turbrelaxexplicit
        end if
 
        ! V2f should only be solved with explicit underrelaxation.
 
        if (equations == ransequations .and. turbmodel == v2f .and. &
            turbrelax == turbrelaximplicit) then
 
            turbrelax = turbrelaxexplicit
 
            if (myid == 0) then
                print "(a)", "#"
                print "(a)", "#                      Warning"
                print "(a)", "# Implicit underrelaxation specified for &
                     &the v2f model."
                print "(a)", "# This is overwritten to explicit &
                     &underrelaxation."
                print "(a)", "#"
            end if
 
        end if
 
        if (nsavevolume <= 0) then
            select case (equationmode)
            case (steady, timespectral)
                nsavevolume = ncycles + ncyclescoarse + nsgstartup + 1
 
            case (unsteady)
                nsavevolume = ntimestepsfine + ntimestepscoarse &
                              + ntimestepsrestart + 1
            end select
        end if
 
        if (nsavesurface <= 0) nsavesurface = nsavevolume
 
        if (eddyvisinfratio < zero) then
 
            ! Default value depends on the turbulence model.
 
            select case (turbmodel)
 
            case (spalartallmaras, spalartallmarasedwards)
                eddyvisinfratio = 0.009_realtype
 
            case default
                eddyvisinfratio = 0.1_realtype
 
            end select
        end if
        !
        !       Determine the number of old grid levels needed for the BDF
        !       time integration of unsteady problems and allocate the memory
        !       for the coefficients. The actual values are not yet set,
        !       because in the first (and possibly second) time step a reduced
        !       order must be used, because the older states are not available
        !       yet. Also allocate the memory for the logicals to indicate
        !       whether or not old solutions have been written.
        !       If a Runge Kutta scheme must be used for the time integration,
        !       either explicit or implicit, a separate routine is called to
        !       set all the necessary variables.
        !
        select case (timeintegrationscheme)
        case (bdf, md)
 
            ! First check if the accuracy is okay.
 
            if (timeaccuracy > thirdorder) then
                if (myid == 0) then
                    print "(a)", "#"
                    print "(a)", "#                      Warning"
                    print "(a)", "# Maximum third order possible for BDF."
                    print "(a)", "# Order has been reduced to third."
                    print "(a)", "#"
                end if
 
                timeaccuracy = thirdorder
            end if
 
            ! Determine the accuracy and set nOldLevels accordingly.
 
            select case (timeaccuracy)
            case (firstorder)
                noldlevels = 1
 
            case (secondorder)
                noldlevels = 2
 
            case (thirdorder)
                noldlevels = 3
            end select
 
            ! Allocate the memory for coefTime.
            if (allocated(coeftime)) deallocate (coeftime)
            allocate (coeftime(0:noldlevels), stat=ierr)
            if (ierr /= 0) &
                call terminate("checkInputParam", &
                               "Memory allocation error for coefTime")
 
            ! Determine the accuracy and set ALE parameters accordingly.
            if (useale) then
                select case (timeaccuracy)
                case (firstorder)
                    nalemeshes = 1
                    nalesteps = 2
 
                case (secondorder)
                    nalemeshes = 2
                    nalesteps = 4
 
                case (thirdorder)
                    call terminate("checkInputParam", &
                                   "ALE can only use 1st and 2nd order time accuracy")
                end select
 
                if (allocated(coeftimeale)) deallocate (coeftimeale)
                allocate (coeftimeale(1:nalesteps), stat=ierr)
                if (ierr /= 0) &
                    call terminate("checkInputParam", &
                                   "Memory allocation error for coefTimeALE")
 
                if (allocated(coefmeshale)) deallocate (coefmeshale)
                allocate (coefmeshale(1:nalemeshes, 2), stat=ierr)
                if (ierr /= 0) &
                    call terminate("checkInputParam", &
                                   "Memory allocation error for coefMeshALE")
            end if
 
            !===============================================================
 
        case (explicitrk)
            noldlevels = 1
            call setstagecoeffexplicitrk
 
        end select
 
        ! Set the logicals whether or not the old solutions have been
        ! written. Note that this is only used for the second and
        ! higher order BDF schemes. However it is allocated with a
        ! minimum size of 1 to avoid problems.
 
        oldsolwrittensize = max(noldlevels - 1_inttype, 1_inttype)
 
        !check allocations for multipile succesive calls
        if (allocated(oldsolwritten)) deallocate (oldsolwritten)
        allocate (oldsolwritten(oldsolwrittensize), stat=ierr)
        if (ierr /= 0) &
            call terminate("checkInputParam", &
                           "Memory allocation error for oldSolWritten")
 
        do nn = 1, oldsolwrittensize
            oldsolwritten(nn) = .false.
        end do
        !
        !       Determine the values of the runge kutta parameters, depending
        !       on the number of stages specified.
        !
        ! Limit the number of stages between 1 and 6 and allocate the
        ! memory.
 
        nrkstages = min(6_inttype, max(1_inttype, nrkstages))
 
        !check allocations for multipile succesive calls
        if (allocated(etark)) deallocate (etark)
        if (allocated(cdisrk)) deallocate (cdisrk)
 
        allocate (etark(nrkstages), cdisrk(nrkstages), stat=ierr)
        if (ierr /= 0) &
            call terminate("checkInputParam", &
                           "Memory allocation error for etaRK and cdisRK")
 
        ! Determine the case we are having here.
 
        select case (nrkstages)
        case (1_inttype)
            etark(1) = one
 
            cdisrk(1) = one
 
        case (2_inttype)
            etark(1) = 0.2222_realtype
            etark(2) = one
 
            cdisrk(1) = one
            cdisrk(2) = one
 
        case (3_inttype)
            etark(1) = 0.2846_realtype
            etark(2) = 0.6067_realtype
            etark(3) = one
 
            cdisrk(1) = one
            cdisrk(2) = one
            cdisrk(3) = one
 
        case (4_inttype)
            etark(1) = 0.33333333_realtype
            etark(2) = 0.26666667_realtype
            etark(3) = 0.55555555_realtype
            etark(4) = one
 
            cdisrk(1) = one
            cdisrk(2) = half
            cdisrk(3) = zero
            cdisrk(4) = zero
 
        case (5_inttype)
            etark(1) = fourth
            etark(2) = 0.16666667_realtype !1/6
            etark(3) = 0.37500000_realtype !3/8
            etark(4) = half
            etark(5) = one
 
            cdisrk(1) = one
            cdisrk(2) = zero
            cdisrk(3) = 0.56_realtype
            cdisrk(4) = zero
            cdisrk(5) = 0.44_realtype
 
        case (6_inttype)
            etark(1) = 0.0722_realtype
            etark(2) = 0.1421_realtype
            etark(3) = 0.2268_realtype
            etark(4) = 0.3425_realtype
            etark(5) = 0.5349_realtype
            etark(6) = one
 
            cdisrk(1) = one
            cdisrk(2) = one
            cdisrk(3) = one
            cdisrk(4) = one
            cdisrk(5) = one
            cdisrk(6) = one
        end select
        !
        !       To avoid any problems later on, allocate the memory for the
        !       rigid body motion parameters if these values were not present
        !       in the parameter file.
        !
        if (.not. allocated(coefpolxrot)) then
            allocate (coefpolxrot(0:0), stat=ierr)
            if (ierr /= 0) &
                call terminate("checkInputParam", &
                               "Memory allocation failure for coefPolXRot")
            coefpolxrot = zero
        end if
 
        if (.not. allocated(coefpolyrot)) then
            allocate (coefpolyrot(0:0), stat=ierr)
            if (ierr /= 0) &
                call terminate("checkInputParam", &
                               "Memory allocation failure for coefPolYRot")
            coefpolyrot = zero
        end if
 
        if (.not. allocated(coefpolzrot)) then
            allocate (coefpolzrot(0:0), stat=ierr)
            if (ierr /= 0) &
                call terminate("checkInputParam", &
                               "Memory allocation failure for coefPolZRot")
            coefpolzrot = zero
        end if
 
        if (.not. allocated(coscoeffourxrot)) then
            allocate (coscoeffourxrot(0:0), stat=ierr)
            if (ierr /= 0) &
                 call terminate("checkInputParam", &
                 "Memory allocation failure for &
                 &cosCoefFourXRot")
            coscoeffourxrot = zero
        end if
 
        if (.not. allocated(sincoeffourxrot)) then
            allocate (sincoeffourxrot(1), stat=ierr)
            if (ierr /= 0) &
                 call terminate("checkInputParam", &
                 "Memory allocation failure for &
                 &sinCoefFourXRot")
            sincoeffourxrot = zero
        end if
 
        if (.not. allocated(coscoeffouryrot)) then
            allocate (coscoeffouryrot(0:0), stat=ierr)
            if (ierr /= 0) &
                 call terminate("checkInputParam", &
                 "Memory allocation failure for &
                 &cosCoefFourYRot")
            coscoeffouryrot = zero
        end if
 
        if (.not. allocated(sincoeffouryrot)) then
            allocate (sincoeffouryrot(1), stat=ierr)
            if (ierr /= 0) &
                 call terminate("checkInputParam", &
                 "Memory allocation failure for &
                 &sinCoefFourYRot")
            sincoeffouryrot = zero
        end if
 
        if (.not. allocated(coscoeffourzrot)) then
            allocate (coscoeffourzrot(0:0), stat=ierr)
            if (ierr /= 0) &
                 call terminate("checkInputParam", &
                 "Memory allocation failure for &
                 &cosCoefFourZRot")
            coscoeffourzrot = zero
        end if
 
        if (.not. allocated(sincoeffourzrot)) then
            allocate (sincoeffourzrot(1), stat=ierr)
            if (ierr /= 0) &
                 call terminate("checkInputParam", &
                 "Memory allocation failure for &
                 &sinCoefFourZRot")
            sincoeffourzrot = zero
        end if
 
        ! Allocate the memory for cpmin_family. We had to wait until
        ! nTimeIntervalsSpectral was set.
        if (.not. allocated(cpmin_family)) then
            allocate (cpmin_family(ntimeintervalsspectral), stat=ierr)
            if (ierr /= 0) &
                 call terminate("checkInputParam", &
                 "Memory allocation failure for &
                 &cpmin_family")
            cpmin_family = zero
        end if
 
        ! Allocate the memory for sepsenmaxfamily. We had to wait until
        ! nTimeIntervalsSpectral was set.
        if (.not. allocated(sepsenmaxfamily)) then
            allocate (sepsenmaxfamily(ntimeintervalsspectral), stat=ierr)
            if (ierr /= 0) &
                 call terminate("checkInputParam", &
                 "Memory allocation failure for &
                 &sepSenMaxFamily")
            sepsenmaxfamily = zero
        end if
 
    end subroutine checkinputparam
    subroutine setdefaultvalues
        !
        !       setDefaultValues sets the default values for the input
        !       parameters where-ever possible. The parameters that must be
        !       set by the user are initialized such a check can be performed
        !       later.
        !
        use constants
 
        !  --------- Bare imports...too many to list -------
        use inputdiscretization
        use inputio
        use inputiteration
        use inputmotion
        use inputoverset
        use inputparallel
        use inputphysics
        use inputtimespectral
        use inputunsteady
        use inputadjoint
        use inputtsstabderiv
        !  ------------------------------------------------
        use flowvarrefstate, only: lref, lrefspecified, pref, rhoref, &
                                   tinfdim, tref
        use iteration, only: noldsolavail, timespectralgridsnotwritten
        use monitor, only: monmasssliding, ntimestepsrestart, timeunsteadyrestart
        use killsignals, only: fatalfail, routinefailed
        use adjointpetsc, only: adjointpetscvarsallocated, adjointpetscpreprocvarsallocated
        use inputcostfunctions
        implicit none
 
        ! Initialize boundary condition warning print outs
        printbcwarnings = .true.
 
        ! Initialize monitoring the turbulent residuals as well as the
        ! monitoring of mass flow of the sliding interfaces to .false.
 
        mondturb = .false.
        monmasssliding = .false.
 
        ! Initialize the logicals to check whether or not monitoring,
        ! surface output and volume output variables were specified to
        ! .false.
 
        monitorspecified = .false.
        surfaceoutspecified = .false.
        volumeoutspecified = .false.
        isooutspecified = .false.
        !
        !       Set the default values for the discretization parameters.
        !
        spacediscr = none                  ! Serves as a check later on.
        orderturb = firstorder            ! First order discretization.
        ! Of turbulent advective terms.
        riemann = roe
        limiter = nolimiter            ! No limiter in upwind schemes.
        precond = noprecond            ! No preconditioning.
 
        eulerwallbctreatment = normalmomentum   ! Normal momentum equation is
        ! Used to determine ghost
        ! cell pressure.
 
        viscwallbctreatment = constantpressure  ! Normal momentum equation is
        ! Used to determine ghost
        ! cell pressure.
 
        outflowtreatment = constantextrapol ! Constant extrapolation at
        ! outflow boundaries.
 
        spacediscrcoarse = none             ! Serves as a check. If nothing
        riemanncoarse = none             ! is specified the fine grid
        ! parameter is taken.
 
        nonmatchtreatment = nonconservative ! Non conservative treatment
        ! of non-matching block to
        ! block boundaries.
 
        vortexcorr = .false.                ! No vortex correction is
        ! applied.
 
        vis2 = half
        vis4 = one / 64.0_realtype
        vis2coarse = half
 
        dirscaling = .true.                 ! Apply isotropic directional
        adis = two * third              ! scaling in the artificial
        ! dissipation schemes.
 
        hscalinginlet = .false.             ! No total enthalpy scaling.
 
        kappacoef = third
        !
        !       Set the default values for the IO-parameters.
 
        gridfile = ""          ! Serves as a check later on.
 
        checkrestartsol = .true.     ! Restart solution is checked for
        ! correct nonDimensionalization.
 
        newgridfile = ""             ! This will be corrected later on
        solfile = ""             ! if nothing is specified. The
        ! default names depend on the
        ! format used
 
        surfacesolfile = ""          ! This will be corrected later if no
        ! surface solution file is specified.
 
        storerindlayer = .true.     ! No halo cells in solution files.
 
        autoparameterupdate = .true. ! Update the input parameter file
        ! when a restart file is written.
        writecoormeter = .false.     ! Use original coordinate units
        ! when writing solution files.
 
        cpfile = ""                  ! Serves as a check later on.
 
        storeconvinneriter = .false. ! Do not store the convergence of
        ! iterations(inner iterations in unsteady mode).
 
#ifdef USE_SINGLE_PRECISION
        precisiongrid = precisionsingle   ! Default IO precision depends
        precisionsol = precisionsingle   ! on the default floating
        ! point type used. Note that
#else
        precisiongrid = precisiondouble   ! for quadrupole precision the
        precisionsol = precisiondouble   ! IO takes place in double
        ! precision.
#endif
 
        ! Surface solution defaults to single precision
        precisionsurfgrid = precisionsingle
        precisionsurfsol = precisionsingle
 
        !
        !       Set the default values for the iteration parameters.
        !
        ncycles = -1    ! Serves as a check later on.
        nsgstartup = 0    ! No single grid startup iterations.
        nsubiterturb = 0    ! No additional turbulent subiterations.
        nupdatebleeds = 50    ! Update the bleeds every 50 iterations.
 
        nsavevolume = 1      ! Only save at the end of the computation.
        nsavesurface = 1
 
        smoother = none
        nrkstages = 5
        nsubiterations = 1
 
        !resAveraging =  noResAveraging ! No residual averaging.
        resaveraging = noresaveraging
        smoop = 1.5_realtype
 
        turbtreatment = decoupled      ! Decoupled solver for the
        ! turbulent equations
        turbsmoother = adi            ! solved using an adi scheme.
        freezeturbsource = .true.          ! Freeze the coarse grid source
        ! terms for a coupled solver.
        turbrelax = turbrelaxnotdefined    ! Will be set later, depending
        ! on the turbulence model.
 
        cfl = -one                         ! Serves as a check later on.
 
        relaxbleeds = 0.1_realtype         ! Relaxation factor for the
        ! bleed boundary conditions.
 
        alfaturb = 0.8_realtype
        betaturb = -one                    ! Serves as a check later on.
 
        l2conv = 1.e-6_realtype     ! Six orders of magnitude for
        ! convergence.
        l2convcoarse = 1.e-2_realtype      ! Only two on coarse grids in
        ! full mg.
 
        maxl2deviationfactor = 1_realtype
        ncyclescoarse = -1             ! If these parameters are not
        cflcoarse = -one           ! specified the corresponding fine
        ! grid values are taken.
 
        fcoll = one       ! No relaxation when restricting the residuals.
 
        mgboundcorr = bcdirichlet0 ! Zero out the boundary halo's for
        ! the multigrid corrections.
 
        mgstartlevel = -1    ! Start at the coarsest grid of the mg cycle
        ! when no restart is performed.
        mgdescription = "sg" ! Single grid computation.
        !
        !       Set the default values for the motion parameters,
        !       i.e. no motion.
        !
        ! Translation data.
 
        ! Rotation data.
 
        rotpoint = zero
 
        degreepolxrot = -1      ! -1, because the start index is 0.
        degreepolyrot = -1
        degreepolzrot = -1
 
        degreefourxrot = -1     ! -1, because the start index is 0,
        ! at least of the cosine part.
        degreefouryrot = -1
        degreefourzrot = -1
 
        omegafourxrot = zero
        omegafouryrot = zero
        omegafourzrot = zero
 
        ! The logical to determine whether or not a motion is specified.
        ! Initialize it to .false.
 
        gridmotionspecified = .false.
        !
        !       Set the default values for the parallel parameters.
        !
        loadimbalance = 0.1_realtype  ! Allow 10 percent load imbalance.
        splitblocks = .true.        ! Allow the splitting of blocks to
        ! obtain a better load balancing.
        loadbalanceiter = 2           ! Do two iterations
        !
        !       Set the default values for the physics parameters.
        !
        equations = none       ! These are parameters that must be
        equationmode = none        ! specified. If not, the program
        flowtype = none        ! exits.
        turbmodel = none
 
        cpmodel = cpconstant       ! Constant cp.
 
        turbprod = strain          ! Strain is used in the production
        ! term of transport turbulence models.
 
        wallfunctions = .false.    ! No wall functions used.
 
        mach = -one            ! Both parameters must be specified
        reynolds = -one            ! for external flows. The -1. serves
        ! as a check later on.
 
        machcoef = -one            ! If not specified MachCoef will
        ! be set to Mach.
 
        veldirfreestream(1) = one  ! Free stream velocity
        veldirfreestream(2) = zero ! is specified in the
        veldirfreestream(3) = zero ! x-axis direction.
 
        liftdirspecified = .false. ! Lift direction not specified.
 
        reynoldslength = one
        tinfdim = 288.15_realtype
        gammaconstant = 1.4_realtype
        rgasdim = 287.87_realtype
 
        prandtl = 0.72_realtype
        prandtlturb = 0.90_realtype
        pklim = 20.0_realtype
        walloffset = zero
 
        ssuthdim = 110.55_realtype
        musuthdim = 1.716e-5_realtype
        tsuthdim = 273.15_realtype
 
        rvfn = 1                             ! Version 1 of the v2f
        ! model is used.
        rvfb = .true.                        ! An upper bound is used
        ! in the v2f scales.
        eddyvisinfratio = -one               ! Default value depends on
        ! the turbulence model.
        turbintensityinf = 0.001_realtype
 
        surfaceref = one
        lengthref = one
 
        pointref(1) = zero
        pointref(2) = zero
        pointref(3) = zero
 
        momentaxis(1, 1) = zero
        momentaxis(1, 2) = one
        momentaxis(2, 1) = zero
        momentaxis(2, 2) = zero
        momentaxis(3, 1) = zero
        momentaxis(3, 2) = zero
 
        !
        !       Set the default values for the time spectral parameters.
        !
        ntimeintervalsspectral = -1   ! Serves as a check later on.
 
        nunsteadysolspectral = -1       ! Serves as a check later on.
 
        writeunsteadyvolspectral = .false. ! No writing of the files
        writeunsteadysurfspectral = .false. ! for postprocessing.
 
        writeunsteadyrestartspectral = .false. ! No writing of an unsteady
        ! mode restart file.
 
        dtunsteadyrestartspectral = -one    ! Is checked later on.
        !
        !       Set the default values for the unsteady parameters.
        !
        timeaccuracy = secondorder  ! Second order time accuracy.
 
        ntimestepscoarse = -1       ! Serves as a check later on.
        ntimestepsfine = -1       ! Serves as a check later on.
 
        deltat = -one               ! Serves as a check later on.
 
        useale = .true.             ! Use the ALE scheme by default.
 
        updatewalldistanceunsteady = .true.  ! This default value is
        ! overruled for models that
        ! are wall distance free.
        !
        !       The reference state variables. Set them to -1, such that they
        !       can be checked later on.
        !
        pref = -one
        rhoref = -one
        tref = -one
        !
        !       The conversion factor of the grid units to meters. Default 1.
        !
        lref = one
        lrefspecified = .false.
        !
        !       Initialization of some unsteady restart parameters. These will
        !       be overwritten when an actual unsteady restart is performed.
        !
        noldsolavail = 1
        ntimestepsrestart = 0
        timeunsteadyrestart = zero
        !
        !       Variables needed for the writing of grid and solution files.
        !
        timespectralgridsnotwritten = .true.
 
        ! Additional Paramters Requiring Defaults
        printiterations = .true.
        routinefailed = .false.
        fatalfail = .false.
        lumpeddiss = .false.
        approxsa = .false.
        useapproxwalldistance = .false.
        updatewallassociations = .false.
        recomputeoverlapmatrix = .true.
        cfllimit = 3.0
        adjointpetscvarsallocated = .false.
        adjointpetscpreprocvarsallocated = .false.
        usematrixfreedrdw = .false.
    end subroutine setdefaultvalues
 
    subroutine initializeisosurfacevariables(values, nValues)
        !
        !       isoVariables extracts from the given string the extra
        !       iso surface variables to be written to the solution file.
        !
        use constants
        use extraoutput, only: isovalues, isosurfacenames, nisosurface
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: nValues
        real(kind=realtype), dimension(nValues), intent(in) :: values
 
        ! Basically just copy into module
        if (allocated(isovalues)) then
            deallocate (isovalues)
        end if
 
        if (allocated(isosurfacenames)) then
            deallocate (isosurfacenames)
        end if
 
        nisosurface = nvalues
        allocate (isovalues(nisosurface))
        allocate (isosurfacenames(nisosurface))
 
        isovalues = values
 
    end subroutine initializeisosurfacevariables
 
    subroutine setisosurfacevariable(variable, iVar)
 
        ! Set variable to iVar. initializeIsoSurfaceVariables MUST be called
        ! first with the desired number of values to set.
 
        use constants
        use cgnsnames
        use extraoutput
        use communication, only: myid
        use utils, only: echk
        implicit none
        !
        !      Subroutine arguments.
        !
        character(len=*), intent(in) :: variable
        integer(kind=intType) :: iVar
 
        select case (variable)
        case ("rho")
            isosurfacenames(ivar) = cgnsdensity
        case ("vx")
            isosurfacenames(ivar) = cgnsvelx
        case ("vy")
            isosurfacenames(ivar) = cgnsvely
        case ("vz")
            isosurfacenames(ivar) = cgnsvelz
        case ("P")
            isosurfacenames(ivar) = cgnspressure
        case ("mx")
            isosurfacenames(ivar) = cgnsmomx
        case ("my")
            isosurfacenames(ivar) = cgnsmomy
        case ("mz")
            isosurfacenames(ivar) = cgnsmomz
        case ("rvx")
            isosurfacenames(ivar) = cgnsrelvelx
        case ("rvy")
            isosurfacenames(ivar) = cgnsrelvely
        case ("rvz")
            isosurfacenames(ivar) = cgnsrelvelz
        case ("rhoe")
            isosurfacenames(ivar) = cgnsenergy
        case ("temp")
            isosurfacenames(ivar) = cgnstemp
        case ("vort")
            isosurfacenames(ivar) = cgnsvortmagn
        case ("vortx")
            isosurfacenames(ivar) = cgnsvortx
        case ("vorty")
            isosurfacenames(ivar) = cgnsvorty
        case ("vortz")
            isosurfacenames(ivar) = cgnsvortz
        case ("cp")
            isosurfacenames(ivar) = cgnscp
        case ("mach")
            isosurfacenames(ivar) = cgnsmach
        case ("rmach")
            isosurfacenames(ivar) = cgnsrelmach
        case ("macht")
            isosurfacenames(ivar) = cgnsmachturb
        case ("ptloss")
            isosurfacenames(ivar) = cgnsptotloss
        case ("eddy")
            isosurfacenames(ivar) = cgnseddy
        case ("eddyratio")
            isosurfacenames(ivar) = cgnseddyratio
        case ("dist")
            isosurfacenames(ivar) = cgnswalldist
        case ("resrho")
            isosurfacenames(ivar) = cgnsresrho
        case ("shock")
            isosurfacenames(ivar) = cgnsshock
        case ("filteredShock")
            isosurfacenames(ivar) = cgnsfilteredshock
        case default
 
            if (myid == 0) Then
                print *, 'Error: ', variable, 'cannot be used as an isoSurface'
            end if
            call echk(-99, __file__, __line__)
        end select
    end subroutine setisosurfacevariable
 
end module inputparamroutines