solvers.F90

Go to the documentation of this file.

module solvers
contains
 
    subroutine solver
        !
        !       solver is the main subroutine of the solver library.
        !       It controls the full multigrid and the kill signals.
        !
        use constants
        use communication, only: myid
        use inputdiscretization, only: eulerwallbctreatment
        use inputiteration, only: mgstartlevel, printiterations
        use inputphysics, only: equationmode
        use inputunsteady, only: timeintegrationscheme
        use killsignals, only: localsignal, nosignal
        use iteration, only: changing_grid, currentlevel, exchangepressureearly, &
                             groundlevel, noldsolavail, t0solver
        use monitor, only: timeunsteady
        use section, only: nsections
        use utils, only: eulerwallspresent
        use multigrid, only: transfertofinegrid
        use partitioning, only: updatecoorfinemesh
        use commonformats, only: stringint1
        implicit none
        !
        !      Local variables.
        !
        real(kind=realtype), dimension(nSections) :: dtadvance
 
        ! If the normal momentum equation should be used to determine
        ! the pressure in the halo for inviscid walls, find out if there
        ! actually are inviscid walls. If so, set the logical
        ! exchangePressureEarly to .true.; otherwise set it to .false.
 
        if (eulerwallbctreatment == normalmomentum) then
            exchangepressureearly = eulerwallspresent()
        else
            exchangepressureearly = .false.
        end if
 
        ! Connect the kill signals with the appropriate functions.
        ! Initialize localSignal for safety.
        ! Only if signalling is supported.
 
#ifndef USE_NO_SIGNALS
        localsignal = nosignal
        call connect_signals
#endif
 
        ! Determine the reference time for the solver.
 
        t0solver = mpi_wtime()
 
        ! Set timeUnsteady to zero; this is amount of time simulated
        ! in unsteady mode.
 
        timeunsteady = zero
 
        ! Loop over the number of grid levels in the current computation.
        ! Note that the counter variable, groundLevel, is defined in
        ! the module iteration.
 
        do groundlevel = mgstartlevel, 1, -1
 
            ! Solve either the steady or the unsteady equations for this
            ! grid level. The time spectral method can be considered as
            ! a special kind of steady mode.
            select case (equationmode)
            case (steady, timespectral)
                call solvestate
            case (unsteady)
                select case (timeintegrationscheme)
                case (explicitrk)
                    call solverunsteadyexplicitrk
                end select
            end select
 
            ! If this is not the finest grid level, interpolate the
            ! solution to the next finer mesh and write a message that
            ! this is happening. Only processor 0 performs this latter
            ! task.
 
            if (groundlevel > 1) then
 
                currentlevel = groundlevel - 1
 
                if (myid == 0) then
                    if (printiterations) then
                        print "(a)", "#"
                        print stringint1, "# Going down to grid level ", currentlevel
                        print "(a)", "#"
                    end if
                end if
 
                call transfertofinegrid(.false.)
 
                ! Move the coordinates of the new fine grid level into the
                ! correct position. Only for unsteady problems with changing
                ! meshes. Note that the first argument of updateCoorFineMesh
                ! is an array with the time step for section.
 
                if (equationmode == unsteady .and. changing_grid) then
                    dtadvance = timeunsteady
                    call updatecoorfinemesh(dtadvance, 1_inttype)
                end if
 
                ! Reset nOldSolAvail to 1, such that the unsteady
                ! computation on the finer mesh starts with a lower
                ! order scheme.
 
                noldsolavail = 1
            end if
 
        end do
 
        ! Explictly set groundlevel to 1
        groundlevel = 1
 
    end subroutine solver
 
    ! ================================================
    ! Utilities for unsteady simulation
    ! ================================================
    subroutine solverunsteadyinit
        !
        !       Initialize variables related to unsteady simulation.
        !       Some are the same as those in the *solver* subroutine,
        !       while others are specific for unsteady problems.
        !
        use aleutils, only: fillcoor, setlevelale
        use constants
        use inputdiscretization, only: eulerwallbctreatment
        use iteration, only: exchangepressureearly, t0solver
        use killsignals, only: localsignal, nosignal
        use monitor, only: timeunsteady, timestepunsteady, writevolume, writesurface, writegrid
        use utils, only: eulerwallspresent
        implicit none
 
        ! BC treatment for normal momentum equation
        if (eulerwallbctreatment == normalmomentum) then
            exchangepressureearly = eulerwallspresent()
        else
            exchangepressureearly = .false.
        end if
 
        ! Connect the kill signals
#ifndef USE_NO_SIGNALS
        localsignal = nosignal
        call connect_signals
#endif
 
        ! Determine the reference time for the solver.
        t0solver = mpi_wtime()
 
        ! Set time to zero
        timeunsteady = zero
        timestepunsteady = 0
 
        ! Fill up old, xold and volold
        call fillcoor
 
        ! Set all ALE levels by initial configuration
        call setlevelale(-1_inttype)
    end subroutine solverunsteadyinit
 
    subroutine updateunsteadygeometry
        !
        !       Update quantities related to geometry due to modification of mesh
        !       That could happen when
        !       - Steady mode with mesh modification
        !       - Unsteady mode with non-moving mesh but with prescribed mesh motion
        !       - Unsteady mode with warping and/or rigidly moving mesh
        !       - Unsteady mode coupled with an external solver
        !
        use aleutils, only: storecoor, interpcoor, recovercoor, setlevelale, &
                            slipvelocitiesfinelevel_ale
        use bcdata, only: setbcdatafinegrid, setbcdatacoarsegrid
        use constants
        use inputmotion, only: gridmotionspecified
        use inputunsteady, only: deltat, updatewalldistanceunsteady, useale
        use iteration, only: changing_grid, deforming_grid, currentlevel, groundlevel, nalemeshes
        use monitor, only: timeunsteady, timeunsteadyrestart
        use partitioning, only: updatecoorfinemesh
        use preprocessingapi, only: shiftcoorandvolumes, metric, &
                                    updatecoordinatesalllevels, updatemetricsalllevels, facerotationmatrices
        use section, only: nsections
        use solverutils, only: gridvelocitiesfinelevel, gridvelocitiescoarselevels, &
                               gridvelocitiesfinelevelpart1, gridvelocitiesfinelevelpart2, &
                               normalvelocitiesalllevels, slipvelocitiesfinelevel, slipvelocitiescoarselevels
        use walldistance, only: updatewalldistancealllevels
        implicit none
        !
        !      Local variables.
        !
        integer(kind=intType) :: nn
        real(kind=realtype), dimension(nSections) :: tnewsec, deltatsec
        integer(kind=intType) :: lale
 
        testchanging: if (changing_grid .or. gridmotionspecified) then
 
            ! Set the new time for all sections
            do nn = 1, nsections
                tnewsec(nn) = timeunsteady + timeunsteadyrestart
                deltatsec(nn) = deltat
            end do
 
            ! For prescribed motion only
            if (gridmotionspecified) &
                call updatecoorfinemesh(deltatsec, 1_inttype)
 
            ! Adapt the geometric info on all grid levels needed for the
            ! current ground level and multigrid cycle.
            ! The wall distance only needs to be recomputed when the
            ! grid is changing; not when a rigid body motion is
            ! specified. Furthermore, the user can choose not to update
            ! the wall distance, because he may know a priori that the
            ! changes in geometry happen quite far away from the boundary
            ! layer. This is accomplished via updateWallDistanceUnsteady.
 
            call updatecoordinatesalllevels
            if (changing_grid .and. updatewalldistanceunsteady) &
                call updatewalldistancealllevels
            call updatemetricsalllevels
 
            ! Update the rotation matrices of the faces. Only needed
            ! on the finest grid level.
            ! For prescribed motion only
 
            if (gridmotionspecified) &
                call facerotationmatrices(currentlevel, .false.)
 
            if (useale) then
                ! Update the velocities using ALE scheme if moving mesh is present
 
                ! First update cell and surface velocity, both are vectors
                ! Only quantities in blocks are updated, and they will not
                ! be interpolated
 
                call gridvelocitiesfinelevelpart1(deforming_grid, tnewsec, 1_inttype)
 
                ! Secondly store x to a temporary variable xALE
 
                call storecoor
 
                ! Thirdly update surface normal and normal velocity
 
                aleloop: do lale = 1, nalemeshes
                    ! Interpolate mesh over latest time step for all ALE Meshes
                    call interpcoor(lale)
 
                    ! Update s[I,J,K], norm
                    call metric(groundlevel)
 
                    ! Update sFace[I,J,K]
                    call gridvelocitiesfinelevelpart2(deforming_grid, tnewsec, 1_inttype)
 
                    ! Update uSlip
                    call slipvelocitiesfinelevel_ale(deforming_grid, tnewsec, 1_inttype)
 
                    ! Update coarse level quantities to make sure multigrid is working
                    call gridvelocitiescoarselevels(1_inttype)
                    call slipvelocitiescoarselevels(1_inttype)
 
                    ! Update rFace
                    call normalvelocitiesalllevels(1_inttype)
 
                    ! Store data to *lale* ALE level
                    call setlevelale(lale)
                end do aleloop
 
                ! Lastly recover x from temporary variable
                ! Then compute data for current level
 
                call recovercoor
 
                ! Finish the rest of the update
                call metric(groundlevel)
                call gridvelocitiesfinelevelpart2(deforming_grid, tnewsec, 1_inttype)
                call slipvelocitiesfinelevel_ale(deforming_grid, tnewsec, 1_inttype)
 
            else
                ! Otherwise update the velocities naively
 
                ! Determine the velocities of the cell centers and faces
                ! for the current ground level. Note that the spectral mode
                ! is always 1 for unsteady mode.
 
                call gridvelocitiesfinelevel(deforming_grid, tnewsec, 1_inttype)
 
                ! Determine the new slip velocities on the viscous walls.
 
                call slipvelocitiesfinelevel(deforming_grid, tnewsec, 1_inttype)
            end if
 
            ! After velocity computations on finest level are done,
            ! Update those on coarser levels
 
            call gridvelocitiescoarselevels(1_inttype)
            call slipvelocitiescoarselevels(1_inttype)
 
            ! Compute the normal velocities of the boundaries, if
            ! needed for the corresponding boundary condition.
 
            call normalvelocitiesalllevels(1_inttype)
 
            ! Determine the prescribed boundary condition data for the
            ! boundary subfaces. First for the (currently) finest grid
            ! and then iteratively for the coarser grid levels.
 
            call setbcdatafinegrid(.false.)
            call setbcdatacoarsegrid
 
        end if testchanging
 
    end subroutine updateunsteadygeometry
 
    subroutine solverunsteadystep
        !
        !       Solve for next time step in unsteady simulation
        !
        use iteration, only: noldsolavail
        use solverutils, only: shiftsolution
        use utils, only: setcoeftimeintegrator
        use walldistance, only: updatewalldistancealllevels
        implicit none
 
        ! Shift the old solution for the new time step.
 
        call shiftsolution
 
        ! Set the coefficients for the time integrator.
 
        call setcoeftimeintegrator
 
        ! Solve the state for the current time step and
        ! update nOldSolAvail.
        call solvestate
        noldsolavail = noldsolavail + 1
 
    end subroutine solverunsteadystep
 
    ! ================================================
    ! The following are not interfaced with Python
    ! ================================================
 
    subroutine checkwriteunsteadyinloop
        !
        !       checkWriteUnsteadyInLoop checks if a solution must be
        !       written inside the time loop and if so write it.
        !
        use communication, only: adflow_comm_world
        use constants
        use inputio, only: liftdistributionfile, slicesolfile
        use inputiteration, only: nsavesurface, nsavevolume
        use inputmotion, only: gridmotionspecified
        use iteration, only: changing_grid, groundlevel, noldlevels, oldsolwritten
        use killsignals, only: localsignal, globalsignal, signalwrite, signalwritequit
        use monitor, only: ntimestepsrestart, timestepunsteady, timeunsteadyrestart, &
                           writevolume, writesurface, writegrid
        use tecplotio, only: writeliftdistributionfile, writeslicesfile
        use surfacefamilies, only: fullfamlist
        implicit none
        !
        !      Local variables.
        !
        integer :: ierr
 
        integer(kind=intType) :: nn
        character(len=7) :: intString
 
        ! Determine the global kill parameter if signals are supported.
 
#ifndef USE_NO_SIGNALS
        call mpi_allreduce(localsignal, globalsignal, 1, adflow_integer, &
                           mpi_max, adflow_comm_world, ierr)
#endif
 
        ! Initialize the logicals for the writing to .false.
 
        writevolume = .false.
        writesurface = .false.
        writegrid = .false.
 
        ! Check whether a solution file, either volume or surface, must
        ! be written. Only on the finest grid level in stand alone mode.
 
        !if(standAloneMode .and. groundLevel == 1) then
        if (groundlevel == 1) then
            if (mod(timestepunsteady, nsavevolume) == 0) &
                writevolume = .true.
            if (mod(timestepunsteady, nsavesurface) == 0) &
                writesurface = .true.
 
            if (globalsignal == signalwrite .or. &
                globalsignal == signalwritequit) then
                writevolume = .true.
                writesurface = .true.
            end if
 
            ! Determine whether or not a grid file must be written.
 
            if (changing_grid .or. gridmotionspecified) &
                writegrid = writevolume
 
            ! Write the solution.
 
            if (writegrid .or. writevolume .or. writesurface) &
                call writesol(fullfamlist, size(fullfamlist))
 
            ! Write the slice files for this timestep if they have been specified. TEMPORARY
            ! Write lift distribution TEMPORARY
            write (intstring, "(i4.4)") timestepunsteady + ntimestepsrestart
            intstring = adjustl(intstring)
            !call writeSlicesFile(trim(slicesolfile)//"_Timestep"//trim(intString)//".dat", .True.)
            !call writeLiftDistributionFile(trim(liftDistributionFile)//"_Timestep"//trim(intString)//".dat", .True.)
 
        end if
 
        ! Update the variable oldSolWritten.
 
        do nn = (noldlevels - 1), 2, -1
            oldsolwritten(nn) = oldsolwritten(nn - 1)
        end do
 
        oldsolwritten(1) = writevolume
 
    end subroutine checkwriteunsteadyinloop
 
    !      ==================================================================
 
    subroutine checkwriteunsteadyendloop
        !
        !       checkWriteUnsteadyEndLoop checks if a solution must be
        !       written at the end of the time loop and if so write it.
        !
        use constants
        use inputmotion, only: gridmotionspecified
        use iteration, only: changing_grid, groundlevel, noldlevels, &
                             oldsolwritten, standalonemode
        use monitor, only: writevolume, writesurface, writegrid
        use surfacefamilies, only: fullfamlist
        implicit none
        !
        !      Local variables.
        !
        integer(kind=intType) :: nn
 
        ! Write a solution. Only on the finest grid in stand alone mode
        ! and if some time steps have been performed since the last write.
 
        if (standalonemode .and. groundlevel == 1) then
 
            writevolume = .not. writevolume
            writesurface = .not. writesurface
 
            if (writevolume .or. writesurface) then
 
                ! Determine whether or not a grid file must be written.
 
                writegrid = .false.
                if (changing_grid .or. gridmotionspecified) &
                    writegrid = writevolume
 
                ! Write the solution.
 
                call writesol(fullfamlist, size(fullfamlist))
 
                ! Update the variable oldSolWritten.
 
                do nn = (noldlevels - 1), 2, -1
                    oldsolwritten(nn) = oldsolwritten(nn - 1)
                end do
 
                oldsolwritten(1) = writevolume
 
            end if
 
        end if
 
    end subroutine checkwriteunsteadyendloop
 
    ! ================================================
    ! Utilities for explicit RK solver
    ! ================================================
 
    subroutine solverunsteadyexplicitrk
        !
        !       solverUnsteadyExplicitRK solves the unsteady equations using
        !       the explicit Runge-Kutta schemes for the multigrid level
        !       groundLevel.
        !
        use constants
        use blockpointers, only: ndom, il, jl, kl, vol, dw, w, dwoldrk, p
        use communication
        use flowvarrefstate
        use inputphysics
        use inputunsteady
        use iteration
        use killsignals
        use monitor
        use utils, only: setpointers
        use flowutils, only: computepressure
        use haloexchange, only: whalo1, whalo2
        use turbutils, only: computeeddyviscosity
        use turbapi, only: turbresidual
        use turbbcroutines, only: applyallturbbc
        use utils, only: convergenceheader
        use residuals, only: initres, residual, sourceterms
        use flowutils, only: computelamviscosity
        use bcroutines, only: applyallbc
        use preprocessingapi, only: updatecoordinatesalllevels
        implicit none
        !
        !      Local parameter.
        !
        integer(kind=intType), parameter :: nWriteConvHeader = 50
        !
        !      Local variables.
        !
        integer(kind=intType) :: iter, nTimeSteps, stage
        integer(kind=intType) :: i, j, k, l, m, nn
 
        real(kind=realtype) :: tmp, timeunsteadyold
 
        character(len=7) :: numberString
 
        ! Set rkStage to 0. This variable is only relevant if Runge-Kutta
        ! smoothers are used as an iterative algorithm, but rkStage needs
        ! to be initialized.
 
        rkstage = 0
 
        ! Initializations of the write and signal parameters.
 
        writevolume = .false.
        writesurface = .false.
        writegrid = .false.
        globalsignal = nosignal
        localsignal = nosignal
 
        ! Initialize timeStepUnsteady to 0 and set the number of
        ! time steps depending on the grid level.
 
        timestepunsteady = 0
 
        ntimesteps = ntimestepscoarse
        if (groundlevel == 1) ntimesteps = ntimestepsfine
 
        ! Write a message to stdout about how many time steps will
        ! be taken on the current grid. Only processor 0 does this.
 
        if (myid == 0) then
            write (numberstring, "(i6)") ntimesteps
            numberstring = adjustl(numberstring)
 
            print "(a)", "#"
            print "(A, I1, 3(A))", "# Grid ", groundlevel, ": Performing ", trim(numberstring), &
                " explicit Runge Kutta time steps."
            print "(a)", "#"
 
            ! Also write the convergence header. Technically this is
            ! not really a convergence history for explicit RK, but
            ! the routine is called this way.
 
            call convergenceheader
        end if
 
        ! Determine and write the initial convergence info. Again not
        ! really convergence for explicit RK.
 
        call convergenceinfo
        !
        !       Loop over the number of time steps to be computed.
        !
        ! Initialize the unsteady time.
 
        timeunsteadyold = timeunsteady
        timeunsteady = timeunsteadyold + gammarkunsteady(1) * deltat
 
        timeloop: do iter = 1, ntimesteps
 
            ! Rewrite the convergence header after a certain number of
            ! iterations. Only processor 0 does this.
 
            if (myid == 0 .and. mod(iter, nwriteconvheader) == 0) &
                call convergenceheader
 
            ! Loop over the number of RK stages.
 
            stageloop: do stage = 1, nrkstagesunsteady
 
                ! Compute the residual. Note that the turbulent residual must
                ! be computed first, because it uses the array of the flow
                ! field residuals as temporary buffers.
 
                if (equations == ransequations) then
                    call initres(nt1, nt2)
                    call turbresidual
                end if
 
                call initres(1_inttype, nwf)
                call sourceterms()
                call residual
 
                ! Loop over the number of domains.
 
                domainloop: do nn = 1, ndom
 
                    ! Set the pointers to this domain. Note that there is only
                    ! one time instance for a time accurate computation.
 
                    call setpointers(nn, currentlevel, 1)
 
                    ! Step 1: Store in dw the conservative update for the
                    !         flow field variables. This means a multiplication
                    !         with deltaT/vol. Also store in w the conservative
                    !         flow field variables.
 
                    do k = 2, kl
                        do j = 2, jl
                            do i = 2, il
                                tmp = deltat / vol(i, j, k)
                                do l = 1, nw
                                    dw(i, j, k, l) = tmp * dw(i, j, k, l)
                                end do
 
                                w(i, j, k, ivx) = w(i, j, k, ivx) * w(i, j, k, irho)
                                w(i, j, k, ivy) = w(i, j, k, ivy) * w(i, j, k, irho)
                                w(i, j, k, ivz) = w(i, j, k, ivz) * w(i, j, k, irho)
 
                            end do
                        end do
                    end do
 
                    ! Step 2: Compute the new conservative flow field variables
                    !         and primitive turbulence variables.
 
                    do m = 1, (stage - 1)
                        tmp = betarkunsteady(stage, m)
                        if (tmp /= zero) then
                            do l = 1, nw
                                do k = 2, kl
                                    do j = 2, jl
                                        do i = 2, il
                                            w(i, j, k, l) = w(i, j, k, l) - tmp * dwoldrk(m, i, j, k, l)
                                        end do
                                    end do
                                end do
                            end do
                        end if
                    end do
 
                    tmp = betarkunsteady(stage, stage)
                    do l = 1, nw
                        do k = 2, kl
                            do j = 2, jl
                                do i = 2, il
                                    w(i, j, k, l) = w(i, j, k, l) - tmp * dw(i, j, k, l)
                                end do
                            end do
                        end do
                    end do
 
                    ! Step 3. Convert the conservative variables back to
                    !         primitive variables for this block. Compute the
                    !         laminar and eddy viscosities as well.
 
                    ! Convert the momentum into velocities.
                    ! Also store the total energy in the pressure, such that
                    ! the routine computePressure can be used.
 
                    do k = 2, kl
                        do j = 2, jl
                            do i = 2, il
                                tmp = one / w(i, j, k, irho)
                                w(i, j, k, ivx) = w(i, j, k, ivx) * tmp
                                w(i, j, k, ivy) = w(i, j, k, ivy) * tmp
                                w(i, j, k, ivz) = w(i, j, k, ivz) * tmp
                                p(i, j, k) = w(i, j, k, irhoe)
                            end do
                        end do
                    end do
 
                    ! Compute the pressure.
 
                    call computepressure(2_inttype, il, 2_inttype, jl, &
                                         2_inttype, kl, 0_inttype)
 
                    ! Swap the pressure and total energy, because
                    ! computePressure stores the pressure in the position of
                    ! the total energy.
 
                    do k = 2, kl
                        do j = 2, jl
                            do i = 2, il
                                tmp = p(i, j, k)
                                p(i, j, k) = w(i, j, k, irhoe)
                                w(i, j, k, irhoe) = tmp
                            end do
                        end do
                    end do
 
                    ! Compute the laminar and eddy viscosities.
 
                    call computelamviscosity(.false.)
                    call computeeddyviscosity(.false.)
 
                    ! Step 4. Store dw in dwOldRK if this is not the last
                    !         RK stage.
 
                    if (stage < nrkstagesunsteady) then
                        do l = 1, nw
                            do k = 2, kl
                                do j = 2, jl
                                    do i = 2, il
                                        dwoldrk(stage, i, j, k, l) = dw(i, j, k, l)
                                    end do
                                end do
                            end do
                        end do
                    end if
 
                end do domainloop
 
                ! Exchange the pressure if the pressure must be exchanged
                ! early. Only the first halo's are needed, thus whalo1 is
                ! called.
 
                call whalo1(currentlevel, 1_inttype, 0_inttype, .true., &
                            .false., .false.)
 
                ! Apply the boundary conditions to all blocks.
 
                call applyallbc(.true.)
                if (equations == ransequations) call applyallturbbc(.true.)
 
                ! Exchange the halo data. As we are on the fine grid
                ! the second halo is needed.
 
                call whalo2(currentlevel, 1_inttype, nw, .true., &
                            .true., .true.)
 
                ! Determine the time and initialize the geometrical and
                ! boundary info for next stage, if needed.
 
                if (stage < nrkstagesunsteady) then
                    timeunsteady = timeunsteadyold &
                                   + gammarkunsteady(stage + 1) * deltat
 
                    call initstagerk(stage + 1)
                end if
 
            end do stageloop
 
            ! Increment timeStepUnsteady and update
            ! timeUnsteady with the current time step.
 
            timestepunsteady = timestepunsteady + 1
            timeunsteady = timeunsteadyold + deltat
 
            ! Write the convergence info.
 
            call convergenceinfo
 
            ! Determine the time and initialize the geometrical and
            ! boundary info for the next time step, such that the writing
            ! of grid files with moving geometries is done correctly.
 
            timeunsteadyold = timeunsteady
            timeunsteady = timeunsteadyold + gammarkunsteady(1) * deltat
 
            call initstagerk(1_inttype)
 
            ! Determine whether or not solution files must be written.
 
            call checkwriteunsteadyinloop
 
            ! Exit the loop if the corresponding kill signal
            ! has been received.
 
            if (globalsignal == signalwritequit) exit
 
        end do timeloop
 
        ! Determine whether or not the final solution must be written.
 
        call checkwriteunsteadyendloop
 
    end subroutine solverunsteadyexplicitrk
 
    !=================================================================
 
    subroutine initstagerk(stage)
        !
        !       initStageRK performs the initialization tasks for the
        !       Runge-Kutta schemes in unsteady mode.
        !
        use inputmotion
        use inputunsteady
        use iteration
        use monitor
        use section
        use bcdata, only: setbcdatafinegrid
        use walldistance, only: updatewalldistancealllevels
        use solverutils
        use aleutils
        use preprocessingapi, only: updatemetricsalllevels, &
                                    updatecoordinatesalllevels, facerotationmatrices
        use partitioning, only: updatecoorfinemesh
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: stage
        !
        !      Local variables.
        !
        integer(kind=intType) :: nn
        real(kind=realtype) :: fact
 
        real(kind=realtype), dimension(nSections) :: tnewsec, deltatsec
 
        ! If the grid is changing a whole lot of geometric
        ! info must be adapted.
 
        testchanging: if (changing_grid .or. gridmotionspecified) then
 
            ! Determine the time step relative to the previous stage.
 
            if (stage == 1) then
                fact = deltat * (one - gammarkunsteady(nrkstagesunsteady))
            else
                fact = deltat * (gammarkunsteady(stage) &
                                 - gammarkunsteady(stage - 1))
            end if
 
            ! Set the time for all sections; also store their
            ! time step. They are the same for all sections, but all
            ! of them should be updated because of consistency.
 
            do nn = 1, nsections
                tnewsec(nn) = timeunsteady + timeunsteadyrestart
                deltatsec(nn) = fact
            end do
 
            ! Advance the coordinates 1 time step.
 
            call updatecoorfinemesh(deltatsec, 1_inttype)
 
            ! Adapt the geometric info on all grid levels needed for the
            ! current ground level and multigrid cycle.
            ! The wall distance only needs to be recomputed when the
            ! grid is changing; not when a rigid body motion is
            ! specified. Furthermore, the user can choose not to update
            ! the wall distance, because he may know a priori that the
            ! changes in geometry happen quite far away from the boundary
            ! layer. This is accomplished via updateWallDistanceUnsteady.
 
            call updatecoordinatesalllevels
            if (changing_grid .and. updatewalldistanceunsteady) &
                call updatewalldistancealllevels
 
            call updatemetricsalllevels
 
            ! Update the rotation matrices of the faces. Only needed
            ! on the finest grid level.
 
            call facerotationmatrices(currentlevel, .false.)
 
            ! Determine the velocities of the cell centers and faces
            ! for the current ground level. Note that the spectral mode
            ! is always 1 for unsteady mode.
 
            call gridvelocitiesfinelevel(deforming_grid, tnewsec, 1_inttype)
 
            ! Determine the new slip velocities on the viscous walls.
 
            call slipvelocitiesfinelevel(deforming_grid, tnewsec, 1_inttype)
 
        end if testchanging
 
        ! Determine the prescribed boundary condition data for the
        ! boundary subfaces for the (currently) finest grid.
 
        call setbcdatafinegrid(.false.)
 
    end subroutine initstagerk
 
    ! ================================================
    ! Internal utilities
    ! ================================================
    subroutine solvestate
        !
        !       solveState computes either the steady or unsteady state
        !       vector for the multigrid level groundLevel, which can be
        !       found in the module iteration.
        !
        use constants
        use communication, only: myid, adflow_comm_world
        use nksolver, only: nklsfuncevals, freestreamresset, nk_ls, &
                            nk_switchtol, usenksolver, nk_cfl, getfreestreamresidual, &
                            getcurrentresidual, nkstep, computeresidualnk
        use anksolver, only: ank_switchtol, useanksolver, ank_cfl, ankstep, destroyanksolver
        use inputio, only: forcedliftfile, forcedslicefile, forcedvolumefile, &
                           forcedsurfacefile, solfile, newgridfile, surfacesolfile, convsolfilebasename
        use inputiteration, only: cfl, cflcoarse, miniternum, ncycles, &
                                  ncyclescoarse, nmgsteps, nupdatebleeds, printiterations, rkreset, timelimit
        use iteration, only: cycling, approxtotalits, converged, cflmonitor, &
                             groundlevel, itertot, itertype, currentlevel, rhores0, totalr, t0solver, &
                             rhoresstart, totalr0, totalrfinal, totalrstart, stepmonitor, linresmonitor, ordersconverged
        use killsignals, only: globalsignal, localsignal, nosignal, routinefailed, signalwrite, &
                               signalwritequit
        use monitor, only: writegrid, writesurface, writevolume, writesoleachiter
        use utils, only: allocconvarrays, convergenceheader
        use tecplotio, only: writetecplot
        use multigrid, only: setcyclestrategy, executemgcycle
        use surfacefamilies, only: fullfamlist
        use flowvarrefstate, only: nwf
        use residuals, only: residual, initres, sourceterms
        use solverutils, only: timestep
        implicit none
        !
        !      Local parameter
        !
        integer(kind=intType), parameter :: nWriteConvHeader = 50
        !
        !      Local variables.
        !
        integer :: ierr
        integer(kind=intType) :: nMGCycles
        character(len=7) :: numberString
        logical :: absConv, relConv, firstNK, firstANK, old_writeGrid
        real(kind=realtype) :: nk_switchtol_save, curtime, ordersconvergedold
        character(len=maxStringLen) :: iterFormat
 
        ! Allocate the memory for cycling.
        if (allocated(cycling)) then
            deallocate (cycling)
        end if
        allocate (cycling(nmgsteps), stat=ierr)
 
        ! Some initializations.
 
        converged = .false.
        globalsignal = nosignal
        localsignal = nosignal
        itertot = 0
 
        ! Determine the cycling strategy for this level.
 
        call setcyclestrategy
 
        ! Determine the number of multigrid cycles, which depends
        ! on the current multigrid level.
 
        nmgcycles = ncycles
        if (groundlevel > 1) nmgcycles = ncyclescoarse
 
        ! Allocate (or reallocate) the convergence array for this solveState.
        call allocconvarrays(nmgcycles)
 
        ! Allocate space for storing hisotry of function evaluations for NK
        ! solver with non-monotone line search if necessary
        if (currentlevel == 1 .and. usenksolver .and. nk_ls == nonmonotonelinesearch) then
            allocate (nklsfuncevals(nmgcycles))
        end if
 
        ! Only compute the free stream resisudal once for efficiency on the
        ! fine grid only.
        if (groundlevel == 1) then
            if (.not. freestreamresset) then
                call getfreestreamresidual(rhores0, totalr0)
                freestreamresset = .true.
            end if
        end if
        ! Write a message. Only done by processor 0.
 
        if (myid == 0) then
 
            ! Write a message about the number of multigrid iterations
            ! to be performed.
 
#ifndef USE_COMPLEX
            iterformat = "(A, I1, 3(A), I6, A, ES10.2)"
#else
            iterformat = "(A, I1, 3(A), I6, A, 2ES10.2)"
#endif
 
            write (numberstring, "(i6)") nmgcycles
            numberstring = adjustl(numberstring)
            numberstring = trim(numberstring)
            if (printiterations) then
                print "(a)", "#"
                print iterformat, "# Grid ", groundlevel, ": Performing ", trim(numberstring), &
                    " iterations, unless converged earlier. Minimum required iteration before NK switch: ", &
                    miniternum, ". Switch to NK at totalR of: ", (nk_switchtol * totalr0)
                print "(a)", "#"
            end if
 
            if (printiterations) then
                call convergenceheader
            end if
        end if
 
        ! Initialize the approxiate iteration count. We won't count this
        ! first residual evaluation. This way on the first convergence info
        ! call, "Iter" and "Iter Total" will both display 0.
        approxtotalits = 0
 
        ! we need to re-set the orders converged to 16 as it might have been modified in the previous iteration
        ordersconverged = 16.0_realtype
 
        ! Evaluate the initial residual
        call computeresidualnk
 
        ! Need to run the time step here since the RK/DADI is expecting
        ! the rad{i,j,k} to be computed.
        call timestep(.false.)
 
        ! Extract the rhoResStart and totalRStart
        call getcurrentresidual(rhoresstart, totalrstart)
 
        ! No iteration type for first residual evaluation
        itertype = "    None"
        ! also no CFL, step size, or linear residual for this iteration
        cflmonitor = -1
        stepmonitor = -1
        linresmonitor = -1
 
        ! Determine and write the initial convergence info.
        call convergenceinfo
 
        ! Loop over the maximum number of nonlinear iterations
        firstank = .true.
        firstnk = .true.
 
        ! Save the NKSwitch tol since it may be modified in the loop
        nk_switchtol_save = nk_switchtol
 
        ! Set the converged order factor now:
        ordersconverged = log10(totalr0 / totalrstart)
        ordersconvergedold = ordersconverged
 
        nonlineariteration: do while (approxtotalits < nmgcycles)
 
            ! Update iterTot
            itertot = itertot + 1
 
            if (mod(itertot, nwriteconvheader) == 0 .and. &
                myid == 0 .and. printiterations) then
                call convergenceheader
            end if
 
            ! Defaults for the monitor
            cflmonitor = cfl
            stepmonitor = 1.0
            linresmonitor = -1
 
            ! Determine what type of update to do:
            if (currentlevel > 1) then
 
                ! Coarse grids do RK/DADI always
                call executemgcycle
                cflmonitor = cflcoarse
            else
                if (.not. useanksolver .and. .not. usenksolver .or. (itertot <= miniternum .and. rkreset)) then
 
                    ! Always RK/DADI or a RK startup. Run the MG Cycle
 
                    call executemgcycle
 
                else if (useanksolver .and. .not. usenksolver) then
 
                    ! Approx solver, no NKSolver
 
                    if (totalr > ank_switchtol * totalr0) then
 
                        call executemgcycle
 
                    else
                        call ankstep(firstank)
                        firstank = .false.
                        cflmonitor = ank_cfl
 
                    end if
 
                else if (.not. useanksolver .and. usenksolver) then
 
                    ! NK Solver no approx solver
 
                    if (totalr > nk_switchtol * totalr0) then
 
                        call executemgcycle
 
                    else
 
                        call nkstep(firstnk)
                        firstnk = .false.
                        cflmonitor = -1
 
                    end if
 
                else if (useanksolver .and. usenksolver) then
 
                    ! Both approximate and NK solver.
 
                    if (totalr > ank_switchtol * totalr0) then
 
                        call executemgcycle
 
                    else if (totalr <= ank_switchtol * totalr0 .and. &
                             totalr > nk_switchtol * totalr0) then
 
                        call ankstep(firstank)
                        firstank = .false.
                        firstnk = .true.
                        cflmonitor = ank_cfl
 
                    else
 
                        ! We free the memory for the ANK solver here because ANK solver
                        ! module depends on the NK solver module and we cannot call
                        ! destroyANKSolver within NKSolver itself.
                        if (firstnk) then
                            call destroyanksolver()
                        end if
 
                        call nkstep(firstnk)
                        firstnk = .false.
                        firstank = .true.
                        cflmonitor = -1
 
                    end if
                end if
            end if
 
            if (timelimit > zero) then
                ! Check if we ran out of time but only if we are required to use the timeLimit
                if (myid == 0) then
                    curtime = mpi_wtime() - t0solver
                end if
 
                call mpi_bcast(curtime, 1, adflow_real, 0, adflow_comm_world, ierr)
 
                if (curtime > timelimit) then
                    ! Set the iterTot to the limit directly so that convergence
                    ! info thinks we are just out of cycles
                    approxtotalits = nmgcycles
                end if
            end if
 
            ! Determine and write the convergence info.
            call convergenceinfo
 
            totalrfinal = totalr
 
            ! Update how far we are converged:
            ordersconverged = max(log10(totalr0 / totalr), ordersconvergedold)
            ordersconvergedold = ordersconverged
 
            ! Check for divergence or nan here
            if (routinefailed) then
                exit nonlineariteration
            end if
 
            ! Exit the loop if we are converged
            if (converged) then
                exit nonlineariteration
            end if
 
            ! Check if the bleed boundary conditions must be updated and
            ! do so if needed.
 
            ! Check if we've received a signal:
#ifndef USE_NO_SIGNALS
            call mpi_allreduce(localsignal, globalsignal, 1, &
                               adflow_integer, mpi_max, adflow_comm_world, &
                               ierr)
#endif
 
            if (globalsignal == signalwrite .or. writesoleachiter) then
                ! We have been told to write the solution even though we are not done iterating
 
                ! The grid must be written along with the volume solution solution
                if (writevolume) then
                    ! temporary change the writeGrid option
                    old_writegrid = writegrid
                    writegrid = .true.
                end if
 
                if (writesoleachiter) then
                    write (numberstring, "(i7)") itertot
                    numberstring = adjustl(numberstring)
                    numberstring = trim(numberstring)
                    surfacesolfile = trim(convsolfilebasename)//"_"//trim(numberstring)//"_surf.cgns"
                    newgridfile = trim(convsolfilebasename)//"_"//trim(numberstring)//"_vol.cgns"
                    solfile = trim(convsolfilebasename)//"_"//trim(numberstring)//"_vol.cgns"
                else
                    surfacesolfile = forcedsurfacefile
                    newgridfile = forcedvolumefile
                    solfile = forcedvolumefile
                end if
 
                call writesol(fullfamlist, size(fullfamlist))
 
                ! Also write potential tecplot files. Note that we are not
                ! writing the tecplot surface file so pass in an empty file
                ! name and a nonsense family list.
                call writetecplot(forcedslicefile, .true., forcedliftfile, .true., &
                                  "", .false., [0], 1)
 
                if (writevolume) then
                    ! change the writeGrid option back
                    writegrid = old_writegrid
                end if
 
                ! Reset the signal
                localsignal = nosignal
            end if
 
            if (globalsignal == signalwritequit) then
                exit nonlineariteration
            end if
 
        end do nonlineariteration
 
        ! Restore the switch tol in case it was changed
        nk_switchtol = nk_switchtol_save
 
        ! deallocate space for storing hisotry of function evaluations for NK
        ! solver with non-monotone line search if necessary
        if (currentlevel == 1 .and. usenksolver .and. nk_ls == nonmonotonelinesearch) then
            deallocate (nklsfuncevals)
        end if
 
    end subroutine solvestate
 
    subroutine convergenceinfo
        !
        !       convergenceInfo computes and writes the convergence info to
        !       standard output. In spectral mode a convergence history for
        !       every mode is written.
        !
        use constants
        use cgnsnames
        use block, only: ncellglobal
        use blockpointers, only: ndom
        use communication, only: adflow_comm_world, myid
        use inputiteration, only: printiterations, l2convcoarse, l2conv, l2convrel, &
                                  miniternum, maxl2deviationfactor, ncycles, rkreset
        use inputphysics, only: liftdirection, dragdirection, equationmode, &
                                lengthref, machcoef, surfaceref
        use flowvarrefstate, only: pref, lref, gammainf
        use inputio, only: storeconvinneriter
        use inputtimespectral, only: ntimeintervalsspectral
        use inputunsteady, only: timeintegrationscheme
        use monitor, only: monloc, monglob, nmon, nmonmax, nmonsum, monnames, timedataarray, &
                           showcpu, monref, convarray, timeunsteadyrestart, timearray, timestepunsteady, &
                           timeunsteady, ntimestepsrestart, solverdataarray, solvertypearray
        use iteration, only: groundlevel, currentlevel, itertot, itertype, approxtotalits, &
                             cflmonitor, stepmonitor, t0solver, converged, linresmonitor
        use killsignals, only: routinefailed, frompython
        use iteration, only: rhores, rhoresstart, totalr, totalrstart, totalr0
        use oversetdata, only: oversetpresent
        use utils, only: setpointers, returnfail, maxhdiffmach, maxeddyv, sumresiduals, sumallresiduals
        use genericisnan, only: myisnan
        use surfaceintegrations, only: integratesurfaces
        use zipperintegrations, only: integratezippers
        use surfacefamilies, only: fullfamlist
        implicit none
        !
        !      Local variables.
        !
        integer :: ierr, iConvStdout
 
        integer(kind=intType) :: sps, nn, mm, iConv
 
        real(kind=realtype) :: hdiffmax, machmax
        real(kind=realtype) :: eddyvismax, yplusmax, sepsensor, cavitation, axismoment
        real(kind=realtype) :: sepsensoravg(3)
 
        real(kind=realtype) :: l2convthislevel, fact
        real(kind=realtype), dimension(3) :: cfp, cfv, cmp, cmv
        real(kind=realtype) :: cmpaxis, cmvaxis
        logical :: nanOccurred, writeIterations
        logical :: absConv, relConv
        real(kind=realtype) :: localvalues(nlocalvalues)
        real(kind=realtype) :: funcvalues(ncostfunction)
        ! Determine whether or not the iterations must be written.
 
        writeiterations = .true.
        if (equationmode == unsteady .and. &
            timeintegrationscheme == explicitrk) writeiterations = .false.
 
        ! Initializations
        converged = .false.
        nanoccurred = .false.
 
        ! Set the L2 norm for convergence for this level.
 
        l2convthislevel = l2convcoarse
        if (groundlevel == 1) l2convthislevel = l2conv
 
        ! Loop over the number of spectral solutions.
 
        spectralloop: do sps = 1, ntimeintervalsspectral
 
            ! Initialize the local monitoring variables to zero.
 
            monloc = zero
 
            ! Loop over the blocks.
 
            domains: do nn = 1, ndom
 
                ! Set the pointers for this block.
 
                call setpointers(nn, groundlevel, sps)
 
                ! Compute the forces and moments for this block.  Note that
                ! we zero localValues before each call becuase we are
                ! summing into momLocal.
                localvalues = zero
                call integratesurfaces(localvalues, fullfamlist)
 
                ! Convert to coefficients for monitoring:
                fact = two / (gammainf * machcoef * machcoef &
                              * surfaceref * lref * lref * pref)
                cfp = fact * localvalues(ifp:ifp + 2)
                cfv = fact * localvalues(ifv:ifv + 2)
                fact = fact / (lengthref * lref)
                cmp = fact * localvalues(imp:imp + 2)
                cmv = fact * localvalues(imv:imv + 2)
                yplusmax = localvalues(iyplus)
                ! Determine the maximum values of the monitoring variables
                ! of this block.
 
                call maxhdiffmach(hdiffmax, machmax)
                call maxeddyv(eddyvismax)
 
                ! Loop over the number of monitoring variables.
                nmonitoringvar: do mm = 1, nmon
 
                    ! Determine the monitoring variable and act accordingly.
 
                    select case (monnames(mm))
 
                    case ('totalR')
                        call sumallresiduals(mm)
 
                    case (cgnsl2resrho)
                        call sumresiduals(irho, mm)
 
                    case (cgnsl2resmomx)
                        call sumresiduals(imx, mm)
 
                    case (cgnsl2resmomy)
                        call sumresiduals(imy, mm)
 
                    case (cgnsl2resmomz)
                        call sumresiduals(imz, mm)
 
                    case (cgnsl2resrhoe)
                        call sumresiduals(irhoe, mm)
 
                    case (cgnsl2resnu, cgnsl2resk)
                        call sumresiduals(itu1, mm)
 
                    case (cgnsl2resomega, cgnsl2restau, cgnsl2resepsilon)
                        call sumresiduals(itu2, mm)
 
                    case (cgnsl2resv2)
                        call sumresiduals(itu3, mm)
 
                    case (cgnsl2resf)
                        call sumresiduals(itu4, mm)
 
                    case (cgnscl)
                        monloc(mm) = monloc(mm) &
                                     + (cfp(1) + cfv(1)) * liftdirection(1) &
                                     + (cfp(2) + cfv(2)) * liftdirection(2) &
                                     + (cfp(3) + cfv(3)) * liftdirection(3)
 
                    case (cgnsclp)
                        monloc(mm) = monloc(mm) + cfp(1) * liftdirection(1) &
                                     + cfp(2) * liftdirection(2) &
                                     + cfp(3) * liftdirection(3)
 
                    case (cgnsclv)
                        monloc(mm) = monloc(mm) + cfv(1) * liftdirection(1) &
                                     + cfv(2) * liftdirection(2) &
                                     + cfv(3) * liftdirection(3)
 
                    case (cgnscd)
                        monloc(mm) = monloc(mm) &
                                     + (cfp(1) + cfv(1)) * dragdirection(1) &
                                     + (cfp(2) + cfv(2)) * dragdirection(2) &
                                     + (cfp(3) + cfv(3)) * dragdirection(3)
 
                    case (cgnscdp)
                        monloc(mm) = monloc(mm) + cfp(1) * dragdirection(1) &
                                     + cfp(2) * dragdirection(2) &
                                     + cfp(3) * dragdirection(3)
 
                    case (cgnscdv)
                        monloc(mm) = monloc(mm) + cfv(1) * dragdirection(1) &
                                     + cfv(2) * dragdirection(2) &
                                     + cfv(3) * dragdirection(3)
 
                    case (cgnscfx)
                        monloc(mm) = monloc(mm) + cfp(1) + cfv(1)
 
                    case (cgnscfy)
                        monloc(mm) = monloc(mm) + cfp(2) + cfv(2)
 
                    case (cgnscfz)
                        monloc(mm) = monloc(mm) + cfp(3) + cfv(3)
 
                    case (cgnscmx)
                        monloc(mm) = monloc(mm) + cmp(1) + cmv(1)
 
                    case (cgnscmy)
                        monloc(mm) = monloc(mm) + cmp(2) + cmv(2)
 
                    case (cgnscmz)
                        monloc(mm) = monloc(mm) + cmp(3) + cmv(3)
 
                    case (cgnshdiffmax)
                        monloc(mm) = max(monloc(mm), hdiffmax)
 
                    case (cgnsmachmax)
                        monloc(mm) = max(monloc(mm), machmax)
 
                    case (cgnsyplusmax)
                        monloc(mm) = max(monloc(mm), yplusmax)
 
                    case (cgnseddymax)
                        monloc(mm) = max(monloc(mm), eddyvismax)
 
                    case (cgnssepsensor)
                        monloc(mm) = monloc(mm) + localvalues(isepsensor)
 
                    case (cgnssepsensorksarea)
                        monloc(mm) = monloc(mm) + localvalues(isepsensorksarea)
 
                    case (cgnscavitation)
                        monloc(mm) = monloc(mm) + localvalues(icavitation)
 
                    case (cgnsaxismoment)
                        monloc(mm) = monloc(mm) + localvalues(iaxismoment)
 
                    end select
 
                end do nmonitoringvar
            end do domains
 
            ! Add the corrections from zipper meshes from proc 0
            if (oversetpresent) then
                localvalues = zero
                call integratezippers(localvalues, fullfamlist, sps)
 
                fact = two / (gammainf * machcoef * machcoef &
                              * surfaceref * lref * lref * pref)
                cfp = localvalues(ifp:ifp + 2) * fact
                cfv = localvalues(ifv:ifv + 2) * fact
                fact = fact / (lengthref * lref)
                cmp = localvalues(imp:imp + 2) * fact
                cmv = localvalues(imv:imv + 2) * fact
 
                !Loop over the number of monitoring variables and just modify
                !the ones that need to be updated with the zipper forces we just
                !computed.
                nmonitoringvarzip: do mm = 1, nmon
 
                    ! Determine the monitoring variable and act accordingly.
 
                    select case (monnames(mm))
 
                    case (cgnscl)
                        monloc(mm) = monloc(mm) &
                                     + (cfp(1) + cfv(1)) * liftdirection(1) &
                                     + (cfp(2) + cfv(2)) * liftdirection(2) &
                                     + (cfp(3) + cfv(3)) * liftdirection(3)
 
                    case (cgnsclp)
                        monloc(mm) = monloc(mm) + cfp(1) * liftdirection(1) &
                                     + cfp(2) * liftdirection(2) &
                                     + cfp(3) * liftdirection(3)
 
                    case (cgnsclv)
                        monloc(mm) = monloc(mm) + cfv(1) * liftdirection(1) &
                                     + cfv(2) * liftdirection(2) &
                                     + cfv(3) * liftdirection(3)
 
                    case (cgnscd)
                        monloc(mm) = monloc(mm) &
                                     + (cfp(1) + cfv(1)) * dragdirection(1) &
                                     + (cfp(2) + cfv(2)) * dragdirection(2) &
                                     + (cfp(3) + cfv(3)) * dragdirection(3)
 
                    case (cgnscdp)
                        monloc(mm) = monloc(mm) + cfp(1) * dragdirection(1) &
                                     + cfp(2) * dragdirection(2) &
                                     + cfp(3) * dragdirection(3)
 
                    case (cgnscdv)
                        monloc(mm) = monloc(mm) + cfv(1) * dragdirection(1) &
                                     + cfv(2) * dragdirection(2) &
                                     + cfv(3) * dragdirection(3)
 
                    case (cgnscfx)
                        monloc(mm) = monloc(mm) + cfp(1) + cfv(1)
 
                    case (cgnscfy)
                        monloc(mm) = monloc(mm) + cfp(2) + cfv(2)
 
                    case (cgnscfz)
                        monloc(mm) = monloc(mm) + cfp(3) + cfv(3)
 
                    case (cgnscmx)
                        monloc(mm) = monloc(mm) + cmp(1) + cmv(1)
 
                    case (cgnscmy)
                        monloc(mm) = monloc(mm) + cmp(2) + cmv(2)
 
                    case (cgnscmz)
                        monloc(mm) = monloc(mm) + cmp(3) + cmv(3)
 
                    end select
 
                end do nmonitoringvarzip
            end if
            ! Determine the global sum of the summation monitoring
            ! variables. This is an all reduce since every processor needs to
            ! know the residual to make the same descisions.
 
            if (nmonsum > 0) &
                call mpi_allreduce(monloc, monglob, nmonsum, adflow_real, &
                                   mpi_sum, adflow_comm_world, ierr)
 
            ! Idem for the maximum monitoring variables.
#ifndef USE_COMPLEX
            if (nmonmax > 0) &
                call mpi_allreduce(monloc(nmonsum + 1), monglob(nmonsum + 1), &
                                   nmonmax, adflow_real, mpi_max, adflow_comm_world, ierr)
#else
            if (nmonmax < 0) &
                monglob(nmonsum + 1) = zero
#endif
 
            ! Write the convergence info; only processor 0 does this.
 
            testrootproc: if (myid == 0) then
 
                ! The variables which must always be written.
 
                if (printiterations) then
                    write (*, "(1x,i6,2x)", advance="no") groundlevel
 
                    if (equationmode == unsteady) then
 
                        write (*, "(i6,1x)", advance="no") timestepunsteady + &
                            ntimestepsrestart
                        write (*, "(es12.5,1x)", advance="no") timeunsteady + &
                            timeunsteadyrestart
 
                    else if (equationmode == timespectral) then
 
                        write (*, "(i8,3x)", advance="no") sps
 
                    end if
 
                    if (writeiterations) then
                        write (*, "(i6,1x)", advance="no") itertot
                        write (*, "(i6,1x)", advance="no") approxtotalits
                        write (*, "(a,1x)", advance="no") itertype
 
                        if (storeconvinneriter) then
                            solverdataarray(itertot, sps, 1) = approxtotalits
                            solvertypearray(itertot, sps) = itertype
                        end if
 
                        if (cflmonitor < zero) then
                            ! Print dashes if no cfl term is used, i.e. NK solver
                            write (*, "(a,1x)", advance="no") "    ----  "
                        else
#ifndef USE_COMPLEX
                            write (*, "(es10.2,1x)", advance="no") cflmonitor
                            if (storeconvinneriter) then
                                solverdataarray(itertot, sps, 2) = cflmonitor
                            end if
#else
                            write (*, "(es10.2,1x)", advance="no") real(cflmonitor)
#endif
                        end if
 
                        if (stepmonitor < zero) then
                            ! Print dashes in the first None iteration
                            write (*, "(a,1x)", advance="no") " ---- "
                        else
#ifndef USE_COMPLEX
                            write (*, "(f5.2,2x)", advance="no") stepmonitor
                            if (storeconvinneriter) then
                                solverdataarray(itertot, sps, 3) = stepmonitor
                            end if
#else
                            write (*, "(f5.2,2x)", advance="no") real(stepmonitor)
#endif
                        end if
 
                        if (linresmonitor < zero) then
                            ! For RK/DADI just print dashes
                            write (*, "(a,1x)", advance="no") " ----"
                        else
 
#ifndef USE_COMPLEX
                            write (*, "(f5.3,1x)", advance="no") linresmonitor
                            if (storeconvinneriter) then
                                solverdataarray(itertot, sps, 4) = linresmonitor
                            end if
#else
                            write (*, "(f5.3,1x)", advance="no") real(linresmonitor)
#endif
                        end if
 
                        if (showcpu) then
#ifndef USE_COMPLEX
                            write (*, "(es12.5,1x)", advance="no") mpi_wtime() - t0solver
                            if (storeconvinneriter) then
                                solverdataarray(itertot, sps, 5) = mpi_wtime() - t0solver
                            end if
#else
                            write (*, "(es12.5,1x)", advance="no") real(mpi_wtime() - t0solver)
#endif
 
                        end if
                    end if
                end if
            end if testrootproc
 
            ! Loop over the number of monitoring values.
            variableloop: do mm = 1, nmon
 
                ! The residual variables must be corrected.
 
                select case (monnames(mm))
 
                case (cgnsl2resrho, cgnsl2resmomx, &
                      cgnsl2resmomy, cgnsl2resmomz, &
                      cgnsl2resrhoe, cgnsl2resnu, &
                      cgnsl2resk, cgnsl2resomega, &
                      cgnsl2restau, cgnsl2resepsilon, &
                      cgnsl2resv2, cgnsl2resf)
#ifndef USE_COMPLEX
                    monglob(mm) = sqrt(monglob(mm) / ncellglobal(groundlevel))
#else
                    ! take the square roots separately in complex mode
                    monglob(mm) = cmplx(sqrt(real(monglob(mm) / ncellglobal(groundlevel))), &
                                        sqrt(aimag(monglob(mm) / ncellglobal(groundlevel))))
#endif
 
                    if (monnames(mm) == cgnsl2resrho) then
                        rhores = monglob(mm)
                    end if
                case ('totalR')
#ifndef USE_COMPLEX
                    monglob(mm) = sqrt(monglob(mm))
#else
                    ! take the square roots separately in complex mode
                    monglob(mm) = cmplx(sqrt(real(monglob(mm))), sqrt(aimag(monglob(mm))))
#endif
                    totalr = monglob(mm)
                end select
 
                if (myisnan(monglob(mm))) nanoccurred = .true.
 
                if (myid == 0 .and. printiterations) then
                    ! Write the convergence info to stdout.
#ifndef USE_COMPLEX
                    write (*, "(es24.16,1x)", advance="no") monglob(mm)
#else
                    select case (monnames(mm))
 
                    case (cgnsl2resrho, cgnsl2resmomx, &
                          cgnsl2resmomy, cgnsl2resmomz, &
                          cgnsl2resrhoe, cgnsl2resnu, &
                          cgnsl2resk, cgnsl2resomega, &
                          cgnsl2restau, cgnsl2resepsilon, &
                          cgnsl2resv2, cgnsl2resf, 'totalR')
 
                        ! we can do a shorter print for residuals because only the leading few digits
                        ! and the exponents are important anyways
                        write (*, '(es16.8,SP,es17.8E3,"i")', advance="no") monglob(mm)
 
                    case default
                        ! we need to do the regular full print for functionals because they are useful
                        write (*, '(es24.16,SP,es25.16E3,"i")', advance="no") monglob(mm)
                    end select
#endif
                end if
 
                ! Store the convergence info in convArray, if desired.
                if (storeconvinneriter) then
                    convarray(itertot, sps, mm) = monglob(mm)
                end if
            end do variableloop
 
            if (myid == 0 .and. printiterations) then
                ! Write the carriage return.
                print "(1x)"
            end if
 
            ! Determine whether or not the solution is converged.
            ! A distinction must be made between unsteady mode and the
            ! other modes, because in unsteady mode the convergence of
            ! the inner iterations is not necessarily stored.
 
            select case (equationmode)
            case (steady, timespectral)
 
                ! Steady or time spectral mode. The convergence histories
                ! are stored and this info can be used. The logical
                ! converged is set to .false. if the density residual
                ! has not converged yet.
 
                absconv = .false.
                relconv = .false.
 
                ! We make a split here based on if we are operating on a
                ! coarse grid or the finest. On the coarse grid levels, we
                ! l2convThisLevel refers to the relative convergence of
                ! rhoRes.
 
                if (currentlevel /= 1) then
                    if (rhores < l2convthislevel * rhoresstart) then
                        relconv = .true.
                    end if
                else
                    ! We are on the fine level. All checking is done using
                    ! the total residual.
 
                    ! Absolute convergence check
                    if (totalr < l2convthislevel * totalr0) then
                        absconv = .true.
                    end if
 
                    ! Relative check only done on finest level
                    if (totalr < l2convrel * totalrstart) then
                        relconv = .true.
                    end if
 
                    ! If the totla number of iterations is less than the
                    ! RKReset, don't check the residual
                    if (itertot < miniternum .and. rkreset) then
                        relconv = .false.
                        absconv = .false.
                    end if
                end if
 
                ! Combine the two flags.
                if (absconv .or. relconv) then
                    converged = .true.
                end if
 
                !===========================================================
 
            case (unsteady)
 
                ! Unsteady mode. The array convArray may not be present
                ! and therefore something else must be done.
                ! First determine the position in the array timeDataArray
                ! that can be used. For the coarser grids this is 1,
                ! because the time evolution is overwritten. For the fine
                ! mesh the actual position is determined.
 
                nn = 1
                if (groundlevel == 1) &
                    nn = timestepunsteady + ntimestepsrestart
                nn = max(nn, 1_inttype)
 
                ! Make a distinction between the time integration
                ! schemes.
 
                select case (timeintegrationscheme)
 
                case (explicitrk)
 
                    ! Explicit scheme. Simply store the data in the
                    ! convergence arrays.
 
                    timearray(nn) = timeunsteady + timeunsteadyrestart
 
                    ! For explicit schemes the residuals are not
                    ! monitored and therefore the monitoring variables
                    ! can simply be copied.
 
                    do mm = 1, nmon
                        timedataarray(nn, mm) = monglob(mm)
                    end do
 
                    !=======================================================
 
                case (bdf, implicitrk)
 
                    ! An implicit scheme is used and therefore an
                    ! iterative algorithm within every time step.
                    ! The array convArray may not be present and
                    ! therefore
                    ! something else must be done. First determine the
                    ! position in the array timeDataArray that can be
                    ! used.
                    ! For the coarser grids this is 1, because the time
                    ! evolution is overwritten. For the fine mesh the
                    ! actual position is determined.
 
                    ! Determine the situation we have here.
 
                    testinitunsteady: if (itertot == 0) then
 
                        ! This is the initialization phase for this time step.
                        ! Simply copy monGlob into monRef, store the value
                        ! of the physical time and set converged to .false..
 
                        do mm = 1, nmon
                            monref(mm) = monglob(mm)
                        end do
 
                        timearray(nn) = timeunsteady + timeunsteadyrestart
                        converged = .false.
 
                        else testinitunsteady
 
                        ! Iteration for this time step. Store the relative
                        ! convergence for the residual compared to the start
                        ! of this time step; an absolute norm does not give
                        ! any information here. For all other monitoring
                        ! variables the current value is stored.
 
                        do mm = 1, nmon
                            select case (monnames(mm))
 
                            case (cgnsl2resrho, cgnsl2resmomx, &
                                  cgnsl2resmomy, cgnsl2resmomz, &
                                  cgnsl2resrhoe, cgnsl2resnu, &
                                  cgnsl2resk, cgnsl2resomega, &
                                  cgnsl2restau, cgnsl2resepsilon, &
                                  cgnsl2resv2, cgnsl2resf)
 
                                timedataarray(nn, mm) = monglob(mm) &
                                                        / max(monref(mm), eps)
 
                                !===============================================
 
                            case default
 
                                timedataarray(nn, mm) = monglob(mm)
 
                            end select
                        end do
 
                        ! Set the logical converged to .false. if the density
                        ! residual has not converged yet.
 
                        if (timedataarray(nn, 1) < l2convthislevel) &
                            converged = .true.
 
                    end if testinitunsteady
                end select ! temporal integration scheme
 
            end select ! unsteady
 
        end do spectralloop
 
        if (nanoccurred) then
 
            if (myid == 0) then
                print *, 'Nan occured in Convergence Info on proc:', myid
            end if
 
            routinefailed = .true.
 
            call returnfail("convergenceInfo", &
                            "A NaN occurred during the computation.")
 
            ! in a normal computation, code will simply exit.
            ! in a python based computation, code will set
            ! routinedFailed to .True. and return to the
            ! python level...
            return
        end if
 
        ! ! If we are at the max iteration limit but the residual is
        ! ! *close*, ie within maxL2DeviationFactor we say that's fine
 
        if (frompython .and. groundlevel == 1 .and. approxtotalits >= ncycles) then
 
            !Check to see if residuals are diverging or stalled for python
            select case (equationmode)
 
            case (steady, timespectral)
 
                ! Steady or time spectral mode. The convergence histories
                ! are stored and this info can be used. If the residuals
                ! are diverging the, logical routineFailed in killSignals
                ! is set to true and the progress is halted.
                !only check on root porcessor
                if (myid == 0) then
 
                    ! If we made it to ncycles, check to see if we're
                    ! "close" to being converged.
                    do sps = 1, ntimeintervalsspectral
                        if (totalr > maxl2deviationfactor * totalr0 * l2convthislevel) then
                            routinefailed = .true.
                        end if
                    end do
                end if
                call mpi_bcast(routinefailed, 1, mpi_logical, 0, adflow_comm_world, ierr)
 
            end select
        end if
 
        ! Determine whether or not the solution is considered
        ! converged.  This info is only known at processor 0 and must
        ! therefore be broadcast to the other processors.
 
        call mpi_bcast(converged, 1, mpi_logical, 0, adflow_comm_world, ierr)
 
    end subroutine convergenceinfo
 
end module solvers