masterRoutines.F90

Go to the documentation of this file.

module masterroutines
contains
    subroutine master(useSpatial, famLists, funcValues, forces, &
                      bcDataNames, bcDataValues, bcDataFamLists)
 
        use constants
        use communication, only: adflow_comm_world
        use bcroutines, only: applyallbc_block
        use bcdata, only: setbcdata, setbcdatafinegrid
        use turbbcroutines, only: applyallturbbcthisblock, bcturbtreatment
        use iteration, only: currentlevel
        use inputadjoint, only: viscpc
        use flowvarrefstate, only: nwf, nw
        use blockpointers, only: ndom, il, jl, kl
        use flowvarrefstate, only: viscous
        use inputphysics, only: turbprod, equationmode, equations, turbmodel
        use inputdiscretization, only: lowspeedpreconditioner, lumpeddiss, spacediscr, useapproxwalldistance
        use inputtimespectral, only: ntimeintervalsspectral
        use initializeflow, only: referencestate
        use section, only: sections, nsections
        use monitor, only: timeunsteadyrestart
        use sa, only: sasource, saviscous, saresscale, qq
        use haloexchange, only: exchangecoor, whalo2
        use walldistance, only: updatewalldistancesquickly
        use solverutils, only: timestep_block
        use flowutils, only: allnodalgradients, computelamviscosity, computepressuresimple, &
                             computespeedofsoundsquared, adjustinflowangle
        use fluxes, only: invisciddissfluxscalarapprox, invisciddissfluxmatrixapprox, &
                          inviscidupwindflux, invisciddissfluxscalar, invisciddissfluxmatrix, &
                          viscousflux, viscousfluxapprox, inviscidcentralflux
        use utils, only: setpointers, echk
        use turbutils, only: turbadvection, computeeddyviscosity
        use residuals, only: initres_block, sourceterms_block
        use surfaceintegrations, only: getsolution
        use adjointextra, only: volume_block, metric_block, boundarynormals, &
                                xhalo_block, sumdwandfw, resscale
        use oversetdata, only: oversetpresent
        use inputoverset, only: oversetupdatemode
        use oversetcommutilities, only: updateoversetconnectivity
        use actuatorregiondata, only: nactuatorregions, actuatorregions
        use walldistancedata, only: xsurfvec, xsurf
        use actuatorregion, only: computeactuatorregionvolume
 
        implicit none
 
        ! Input Arguments:
        logical, intent(in) :: useSpatial
        integer(kind=intType), optional, dimension(:, :), intent(in) :: famLists
        real(kind=realtype), optional, dimension(:, :), intent(out) :: funcvalues
        character, optional, dimension(:, :), intent(in) :: bcDataNames
        real(kind=realtype), optional, dimension(:), intent(in) :: bcdatavalues
        integer(kind=intType), optional, dimension(:, :) :: bcDataFamLists
 
        ! Output Variables
        real(kind=realtype), intent(out), optional, dimension(:, :, :) :: forces
 
        ! Working Variables
        integer(kind=intType) :: ierr, nn, sps, fSize, iRegion
        real(kind=realtype), dimension(nSections) :: time
        real(kind=realtype) :: dummyreal
 
        if (usespatial) then
            call adjustinflowangle()
 
            ! Update all the BCData
            call referencestate
            if (present(bcdatanames)) then
                do sps = 1, ntimeintervalsspectral
                    call setbcdata(bcdatanames, bcdatavalues, bcdatafamlists, sps, &
                                   size(bcdatavalues), size(bcdatafamlists, 2))
                end do
                call setbcdatafinegrid(.true.)
            end if
 
            do sps = 1, ntimeintervalsspectral
                do nn = 1, ndom
                    call setpointers(nn, 1, sps)
                    call xhalo_block()
                end do
            end do
 
            ! Now exchange the coordinates (fine level only)
            call exchangecoor(1)
 
            do sps = 1, ntimeintervalsspectral
                ! Update overset connectivity if necessary
                if (oversetpresent .and. oversetupdatemode == updatefast) then
                    call updateoversetconnectivity(1_inttype, sps)
                end if
            end do
        end if
 
        if (usespatial) then
            ! Zero out the local volume pointers for the actuator zone
            do iregion = 1, nactuatorregions
                actuatorregions(iregion)%volLocal = zero
            end do
        end if
 
        do sps = 1, ntimeintervalsspectral
            do nn = 1, ndom
                call setpointers(nn, 1, sps)
 
                if (usespatial) then
 
                    call vecgetarrayf90(xsurfvec(1, sps), xsurf, ierr)
                    call echk(ierr, __file__, __line__)
 
                    call volume_block
 
                    ! Compute the volume of each actuator region
                    do iregion = 1, nactuatorregions
                        call computeactuatorregionvolume(nn, iregion)
                    end do
 
                    call metric_block
                    call boundarynormals
 
                    if (equations == ransequations .and. useapproxwalldistance) then
                        call updatewalldistancesquickly(nn, 1, sps)
                    end if
 
                    ! These arrays need to be restored before we can move to the next spectral instance.
                    call vecrestorearrayf90(xsurfvec(1, sps), xsurf, ierr)
                    call echk(ierr, __file__, __line__)
 
                end if
 
                ! Compute the pressures/viscositites
                call computepressuresimple(.false.)
 
                ! Compute Laminar/eddy viscosity if required
                call computelamviscosity(.false.)
                call computeeddyviscosity(.false.)
 
                ! Make sure to call the turb BC's first incase we need to
                ! correct for K
                if (equations == ransequations) then
                    call bcturbtreatment
                    call applyallturbbcthisblock(.true.)
                end if
                call applyallbc_block(.true.)
            end do
        end do
 
        ! Sum the local actuator zone volumes into the global actuator
        ! zone volumes.
        if (usespatial) then
            do iregion = 1, nactuatorregions
                call mpi_allreduce(actuatorregions(iregion)%volLocal, actuatorregions(iregion)%volume, 1, &
                                   adflow_real, mpi_sum, adflow_comm_world, ierr)
                call echk(ierr, __file__, __line__)
            end do
        end if
 
        ! Exchange values
        call whalo2(currentlevel, 1_inttype, nw, .true., .true., .true.)
 
        ! Need to re-apply the BCs. The reason is that BC halos behind
        ! interpolated cells need to be recomputed with their new
        ! interpolated values from actual compute cells. Only needed for
        ! overset.
        if (oversetpresent) then
            do sps = 1, ntimeintervalsspectral
                do nn = 1, ndom
                    call setpointers(nn, 1, sps)
                    if (equations == ransequations) then
                        call bcturbtreatment
                        call applyallturbbcthisblock(.true.)
                    end if
                    call applyallbc_block(.true.)
                end do
            end do
        end if
 
        ! Main loop for the residual
        do sps = 1, ntimeintervalsspectral
            do nn = 1, ndom
                call setpointers(nn, 1, sps)
                call initres_block(1, nw, nn, sps)
                do iregion = 1, nactuatorregions
                    call sourceterms_block(nn, .true., iregion, dummyreal)
                end do
 
                call timestep_block(.false.)
 
                ! Compute turbulence residual for RANS equations
                if (equations == ransequations) then
 
                    ! Initialize only the Turblent Variables
                    !call unsteadyTurbSpectral_block(itu1, itu1, nn, sps)
 
                    select case (turbmodel)
 
                    case (spalartallmaras)
                        allocate (qq(2:il, 2:jl, 2:kl))
                        call sasource
                        call turbadvection(1_inttype, 1_inttype, itu1 - 1, qq)
                        !call unsteadyTurbTerm(1_intType, 1_intType, itu1-1, qq)
                        call saviscous
                        call saresscale
                        deallocate (qq)
                    end select
                end if
 
                ! Compute the mean flow residuals
                call inviscidcentralflux
                if (lumpeddiss) then
                    select case (spacediscr)
                    case (dissscalar)
                        call invisciddissfluxscalarapprox
                    case (dissmatrix)
                        call invisciddissfluxmatrixapprox
                    case (upwind)
                        call inviscidupwindflux(.true.)
                    end select
                else
                    select case (spacediscr)
                    case (dissscalar)
                        call invisciddissfluxscalar
                    case (dissmatrix)
                        call invisciddissfluxmatrix
                    case (upwind)
                        call inviscidupwindflux(.true.)
                    end select
                end if
 
                if (viscous) then
                    call computespeedofsoundsquared
                    if (.not. lumpeddiss .or. viscpc) then
                        call allnodalgradients
                        call viscousflux
                    else
                        call viscousfluxapprox
                    end if
                end if
                call sumdwandfw
                ! if (lowSpeedPreconditioner) then
                !    call applyLowSpeedPreconditioner
                ! end if
                call resscale
            end do
        end do
 
        ! Compute the final solution values
        if (present(famlists)) then
            call getsolution(famlists, funcvalues)
        end if
 
        do sps = 1, ntimeintervalsspectral
            if (present(forces)) then
                ! Now we can retrieve the forces/tractions for this spectral instance
                fsize = size(forces, 2)
                call getforces(forces(:, :, sps), fsize, sps)
            end if
        end do
 
    end subroutine master
#ifndef USE_COMPLEX
    subroutine master_d(wdot, xdot, forcesDot, dwDot, famLists, funcValues, funcValuesd, &
                        bcDataNames, bcDataValues, bcDataValuesd, bcDataFamLists)
 
        use constants
        use diffsizes, only: isize1ofdrfbcdata, isize1ofdrfviscsubface
        use communication, only: adflow_comm_world, myid
        use iteration, only: currentlevel
        use bcextra_d, only: applyallbc_block_d
        use inputadjoint, only: viscpc
        use flowvarrefstate, only: nw, nwf
        use blockpointers, only: ndom, il, jl, kl, wd, xd, dw, dwd, nbocos, nviscbocos
        use flowvarrefstate, only: viscous, timerefd
        use inputphysics, only: turbprod, equationmode, equations, turbmodel, walldistanceneeded
        use inputdiscretization, only: lowspeedpreconditioner, lumpeddiss, spacediscr, useapproxwalldistance
        use inputtimespectral, only: ntimeintervalsspectral
        use section, only: sections, nsections
        use monitor, only: timeunsteadyrestart
        use utils, only: iswalltype, setpointers, setpointers_d, echk
        use sa_d, only: sasource_d, saviscous_d, saresscale_d, qq
        use turbutils_d, only: turbadvection_d, computeeddyviscosity_d
        use fluxes_d, only: invisciddissfluxscalarapprox_d, invisciddissfluxmatrixapprox_d, &
                            inviscidupwindflux_d, invisciddissfluxscalar_d, invisciddissfluxmatrix_d, &
                            inviscidupwindflux_d, viscousflux_d, viscousfluxapprox_d, inviscidcentralflux_d
        use residuals_d, only: sourceterms_block_d, initres_block_d
        use adjointpetsc, only: x_like
        use haloexchange, only: whalo2_d, exchangecoor_d, exchangecoor, whalo2
        use walldistance_d, only: updatewalldistancesquickly_d
        use walldistancedata, only: xsurfvec, xsurfvecd, xsurf, xsurfd, wallscatter
        use flowutils_d, only: computepressuresimple_d, computelamviscosity_d, &
                               computespeedofsoundsquared_d, allnodalgradients_d, adjustinflowangle_d
        use solverutils_d, only: timestep_block_d, gridvelocitiesfinelevel_block_d, slipvelocitiesfinelevel_block_d, &
                                 normalvelocities_block_d
        use turbbcroutines_d, only: applyallturbbcthisblock_d, bcturbtreatment_d
        use initializeflow_d, only: referencestate_d
        use surfaceintegrations, only: getsolution_d
        use adjointextra_d, only: xhalo_block_d, volume_block_d, metric_block_d, boundarynormals_d
        use adjointextra_d, only: resscale_d, sumdwandfw_d
        use bcdata, only: setbcdata_d, setbcdatafinegrid_d
        use oversetdata, only: oversetpresent
        use inputoverset, only: oversetupdatemode
        use oversetcommutilities, only: updateoversetconnectivity_d
        use actuatorregiondata, only: nactuatorregions, actuatorregionsd
        use actuatorregion_d, only: computeactuatorregionvolume_d
#include <petsc/finclude/petsc.h>
        use petsc
        implicit none
 
        ! Input Arguments:
        real(kind=realtype), intent(in), dimension(:) :: wdot, xdot
        integer(kind=intType), optional, dimension(:, :), intent(in) :: famLists
        real(kind=realtype), optional, dimension(:, :), intent(out) :: funcvalues, funcvaluesd
        character, optional, dimension(:, :), intent(in) :: bcDataNames
        real(kind=realtype), optional, dimension(:), intent(in) :: bcdatavalues, bcdatavaluesd
        integer(kind=intType), optional, dimension(:, :) :: bcDataFamLists
 
        ! Output variables:
        real(kind=realtype), intent(out), dimension(:) :: dwdot
        real(kind=realtype), intent(out), dimension(:, :, :) :: forcesdot
 
        ! Working Variables
        real(kind=realtype), dimension(:, :, :), allocatable :: forces
        integer(kind=intType) :: ierr, nn, sps, mm, i, j, k, l, fSize, ii, jj, iRegion
        real(kind=realtype), dimension(nSections) :: time
        real(kind=realtype) :: dummyreal, dummyreald
 
        logical :: useOldCoor ! for adjointextra_d() functions
        useoldcoor = .false.
 
        fsize = size(forcesdot, 2)
        allocate (forces(3, fsize, ntimeintervalsspectral))
 
        call vecplacearray(x_like, xdot, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Set the provided w and x seeds:
        ii = 0
        jj = 0
        domainloop1: do nn = 1, ndom
            spectalloop1: do sps = 1, ntimeintervalsspectral
                call setpointers_d(nn, 1, sps)
                do k = 1, kl
                    do j = 1, jl
                        do i = 1, il
                            do l = 1, 3
                                ii = ii + 1
                                xd(i, j, k, l) = xdot(ii)
                            end do
                        end do
                    end do
                end do
                do k = 2, kl
                    do j = 2, jl
                        do i = 2, il
                            do l = 1, nw
                                jj = jj + 1
                                wd(i, j, k, l) = wdot(jj)
                            end do
                        end do
                    end do
                end do
            end do spectalloop1
        end do domainloop1
 
        do sps = 1, ntimeintervalsspectral
            do nn = 1, ndom
                call setpointers_d(nn, 1, sps)
                call xhalo_block_d()
            end do
        end do
 
        ! Now exchange the coordinates. Note that we *must* exhchange the
        ! actual coordinates as well becuase xhao_block overwrites all
        ! halo nodes and exchange coor corrects them.
        call exchangecoor_d(1)
        call exchangecoor(1)
 
        do sps = 1, ntimeintervalsspectral
            ! Update overset connectivity if necessary
            if (oversetpresent) then
                if (oversetupdatemode == updatefast) then
                    call updateoversetconnectivity_d(1_inttype, sps)
                else if (oversetupdatemode == updatefull) then
                    if (myid == 0) then
                        print *, 'Full overset update derivatives not implemented'
                    end if
                end if
            end if
        end do
 
        ! Now set the xsurfd contribution from the full x perturbation.
        ! scatter from the global seed (in x_like) to xSurfVecd...but only
        ! if wallDistances were used
        if (walldistanceneeded .and. useapproxwalldistance) then
            do sps = 1, ntimeintervalsspectral
                call vecscatterbegin(wallscatter(1, sps), x_like, xsurfvecd(sps), insert_values, scatter_forward, ierr)
                call echk(ierr, __file__, __line__)
 
                call vecscatterend(wallscatter(1, sps), x_like, xsurfvecd(sps), insert_values, scatter_forward, ierr)
                call echk(ierr, __file__, __line__)
            end do
        end if
 
        call adjustinflowangle_d
        call referencestate_d
        if (present(bcdatanames)) then
            do sps = 1, ntimeintervalsspectral
                call setbcdata_d(bcdatanames, bcdatavalues, bcdatavaluesd, &
                                 bcdatafamlists, sps, size(bcdatavalues), size(bcdatafamlists, 2))
            end do
            call setbcdatafinegrid_d(.true.)
        end if
 
        ! Zero out the local volume seeds for the actuator zones
        do iregion = 1, nactuatorregions
            actuatorregionsd(iregion)%volLocal = zero
        end do
 
        do sps = 1, ntimeintervalsspectral
            do nn = 1, ndom
 
                call setpointers_d(nn, 1, sps)
                isize1ofdrfbcdata = nbocos
                isize1ofdrfviscsubface = nviscbocos
 
                ! Get the pointers from the petsc vector for the surface
                ! perturbation for wall distance.
                call vecgetarrayf90(xsurfvec(1, sps), xsurf, ierr)
                call echk(ierr, __file__, __line__)
 
                ! And it's derivative
                call vecgetarrayf90(xsurfvecd(sps), xsurfd, ierr)
                call echk(ierr, __file__, __line__)
 
                call volume_block_d()
                call metric_block_d()
 
                ! Loop over the actuator regions to compute the local
                ! volume seeds
                do iregion = 1, nactuatorregions
                    call computeactuatorregionvolume_d(nn, iregion)
                end do
                call boundarynormals_d()
 
                time = timeunsteadyrestart
                if (equationmode .eq. timespectral) then
                    do mm = 1, nsections
                        time(mm) = time(mm) + (sps - 1) * sections(mm)%timeperiod / real(&
                             &         ntimeintervalsspectral, realtype)
                    end do
                end if
 
                call gridvelocitiesfinelevel_block_d(useoldcoor, time, sps, nn)
                ! required for ts
                call normalvelocities_block_d(sps)
                ! required for ts
                call slipvelocitiesfinelevel_block_d(useoldcoor, time, sps, nn)
 
                if (equations == ransequations .and. useapproxwalldistance) then
                    call updatewalldistancesquickly_d(nn, 1, sps)
                end if
 
                call computepressuresimple_d(.false.)
                call computelamviscosity_d(.false.)
                call computeeddyviscosity_d(.false.)
 
                ! Make sure to call the turb BC's first incase we need to
                ! correct for K
                if (equations == ransequations) then
                    call bcturbtreatment_d
                    call applyallturbbcthisblock_d(.true.)
                end if
 
                call applyallbc_block_d(.true.)
 
                ! These arrays need to be restored before we can move to the next spectral instance.
                call vecrestorearrayf90(xsurfvec(1, sps), xsurf, ierr)
                call echk(ierr, __file__, __line__)
 
                ! And it's derivative
                call vecrestorearrayf90(xsurfvecd(sps), xsurfd, ierr)
                call echk(ierr, __file__, __line__)
            end do
        end do
 
        ! Loop over the actuator regions again to sum the local volume
        ! seeds into the global volume seeds
        do iregion = 1, nactuatorregions
            call mpi_allreduce(actuatorregionsd(iregion)%volLocal, actuatorregionsd(iregion)%volume, 1, &
                               adflow_real, mpi_sum, adflow_comm_world, ierr)
            call echk(ierr, __file__, __line__)
        end do
 
        ! Just exchange the derivative values.
        call whalo2_d(1, 1, nw, .true., .true., .true.)
 
        ! Need to re-apply the BCs. The reason is that BC halos behind
        ! interpolated cells need to be recomputed with their new
        ! interpolated values from actual compute cells. Only needed for
        ! overset.
        if (oversetpresent) then
            do sps = 1, ntimeintervalsspectral
                do nn = 1, ndom
                    call setpointers_d(nn, 1, sps)
                    if (equations == ransequations) then
                        call bcturbtreatment_d
                        call applyallturbbcthisblock_d(.true.)
                    end if
                    call applyallbc_block_d(.true.)
                end do
            end do
        end if
 
        do sps = 1, ntimeintervalsspectral
            do nn = 1, ndom
                call setpointers_d(nn, 1, sps)
                isize1ofdrfbcdata = nbocos
                isize1ofdrfviscsubface = nviscbocos
 
                ! initalize the residuals for this block
                call initres_block_d(1, nw, nn, sps)
 
                ! Compute any source terms
                do iregion = 1, nactuatorregions
                    call sourceterms_block_d(nn, .true., iregion, dummyreal, dummyreald)
                end do
 
                call timestep_block_d(.false.)
 
                !Compute turbulence residual for RANS equations
                if (equations == ransequations) then
                    !call unsteadyTurbSpectral_block(itu1, itu1, nn, sps)
 
                    select case (turbmodel)
                    case (spalartallmaras)
                        call sasource_d
                        call turbadvection_d(1_inttype, 1_inttype, itu1 - 1, qq)
                !!call unsteadyTurbTerm_d(1_intType, 1_intType, itu1-1, qq)
                        call saviscous_d
                        call saresscale_d
                    end select
                end if
 
                ! compute the mean flow residual
                call inviscidcentralflux_d
 
                if (lumpeddiss) then
                    select case (spacediscr)
                    case (dissscalar)
                        call invisciddissfluxscalarapprox_d
                    case (dissmatrix)
                        call invisciddissfluxmatrixapprox_d
                    case (upwind)
                        call inviscidupwindflux_d(.true.)
                    end select
                else
                    select case (spacediscr)
                    case (dissscalar)
                        call invisciddissfluxscalar_d
                    case (dissmatrix)
                        call invisciddissfluxmatrix_d
                    case (upwind)
                        call inviscidupwindflux_d(.true.)
                    end select
                end if
 
                if (viscous) then
                    call computespeedofsoundsquared_d
                    if (.not. lumpeddiss .or. viscpc) then
                        call allnodalgradients_d
                        call viscousflux_d
                    else
                        call viscousfluxapprox_d
                    end if
                end if
 
                ! if (lowSpeedPreconditioner) then
                !    call applyLowSpeedPreconditioner_d
                ! end if
                call sumdwandfw_d
                call resscale_d
            end do
        end do
 
        ! Compute final solution values
        if (present(famlists)) then
            call getsolution_d(famlists, funcvalues, funcvaluesd)
        end if
 
        do sps = 1, ntimeintervalsspectral
            call getforces_d(forces(:, :, sps), forcesdot(:, :, sps), fsize, sps)
        end do
 
        ! Copy out the residual derivative into the provided dwDot
        ii = 0
        do nn = 1, ndom
            do sps = 1, ntimeintervalsspectral
                call setpointers_d(nn, 1, sps)
                do k = 2, kl
                    do j = 2, jl
                        do i = 2, il
                            do l = 1, nw
                                ii = ii + 1
                                dwdot(ii) = dwd(i, j, k, l)
                            end do
                        end do
                    end do
                end do
            end do
        end do
        call vecresetarray(x_like, ierr)
        call echk(ierr, __file__, __line__)
        deallocate (forces)
    end subroutine master_d
 
    subroutine master_b(wbar, xbar, extraBar, forcesBar, dwBar, nState, famLists, &
                        funcValues, funcValuesd, bcDataNames, bcDataValues, bcDataValuesd, bcDataFamLists)
 
        ! This is the main reverse mode differentiaion of master. It
        ! compute the reverse mode sensitivity of *all* outputs with
        ! respect to *all* inputs. Anything that needs to be
        ! differentiated for the adjoint method should be included in this
        ! function. This routine is written by hand, assembling the various
        ! individually differentiated tapenade routines.
 
        use constants
        use adjointvars, only: ialpha, ibeta, imach, imachgrid, itemperature, idensity, &
                               ipointrefx, ipointrefy, ipointrefz, ipressure
        use communication, only: adflow_comm_world, myid
        use iteration, only: currentlevel
        use inputadjoint, only: viscpc
        use fluxes, only: viscousflux
        use flowvarrefstate, only: nw, nwf, viscous, pinfdimd, rhoinfdimd, tinfdimd
        use blockpointers, only: ndom, il, jl, kl, wd, xd, dw, dwd
        use inputphysics, only: pointrefd, alphad, betad, equations, machcoefd, &
                                machd, machgridd, rgasdimd, equationmode, turbmodel, walldistanceneeded
        use inputdiscretization, only: lowspeedpreconditioner, lumpeddiss, spacediscr, useapproxwalldistance
        use inputtimespectral, only: ntimeintervalsspectral
        use inputadjoint, only: frozenturbulence
        use utils, only: iswalltype, setpointers_b, echk
        use adjointpetsc, only: x_like
        use haloexchange, only: whalo2_b, exchangecoor_b, exchangecoor, whalo2
        use walldistancedata, only: xsurfvec, xsurfvecd, xsurf, xsurfd, wallscatter
        use surfaceintegrations, only: getsolution_b
        use flowutils, only: fixallnodalgradientsfromad
        use adjointextra_b, only: resscale_b, sumdwandfw_b
        use adjointextra_b, only: xhalo_block_b, volume_block_b, metric_block_b, boundarynormals_b
        use flowutils_b, only: computepressuresimple_b, computelamviscosity_b, &
                               computespeedofsoundsquared_b, allnodalgradients_b, adjustinflowangle_b
        use solverutils_b, only: timestep_block_b, gridvelocitiesfinelevel_block_b, slipvelocitiesfinelevel_block_b, &
                                 normalvelocities_block_b
        use turbbcroutines_b, only: applyallturbbcthisblock_b, bcturbtreatment_b
        use initializeflow_b, only: referencestate_b
        use walldistance_b, only: updatewalldistancesquickly_b
        use sa_b, only: sasource_b, saviscous_b, saresscale_b, qq
        use turbutils_b, only: turbadvection_b, computeeddyviscosity_b
        use residuals_b, only: sourceterms_block_b, initres_block_b
        use fluxes_b, only: inviscidupwindflux_b, invisciddissfluxscalar_b, &
                            invisciddissfluxmatrix_b, viscousflux_b, inviscidcentralflux_b
        use bcextra_b, only: applyallbc_block_b
        use bcdata, only: setbcdata_b, setbcdatafinegrid_b
        use oversetdata, only: oversetpresent
        use inputoverset, only: oversetupdatemode
        use oversetcommutilities, only: updateoversetconnectivity_b
        use bcroutines, only: applyallbc_block
        use actuatorregiondata, only: nactuatorregions, actuatorregionsd
        use actuatorregion_b, only: computeactuatorregionvolume_b
        use monitor, only: timeunsteadyrestart
        use section, only: sections, nsections ! used in time-declaration
 
#include <petsc/finclude/petsc.h>
        use petsc
        implicit none
 
        ! Input variables:
        real(kind=realtype), intent(in), dimension(:) :: dwbar
        real(kind=realtype), intent(in), dimension(:, :, :) :: forcesbar
        integer(kind=intType), intent(in) :: nState
        integer(kind=intType), optional, dimension(:, :), intent(in) :: famLists
        real(kind=realtype), optional, dimension(:, :) :: funcvalues, funcvaluesd
        character, optional, dimension(:, :), intent(in) :: bcDataNames
        real(kind=realtype), optional, dimension(:), intent(in) :: bcdatavalues
        integer(kind=intType), optional, dimension(:, :) :: bcDataFamLists
 
        ! Output Arguments:
        real(kind=realtype), optional, intent(out), dimension(:) :: wbar, xbar, extrabar
        real(kind=realtype), optional, dimension(:), intent(out) :: bcdatavaluesd
 
        ! Working Variables
        integer(kind=intType) :: ierr, nn, sps, mm, i, j, k, l, fSize, ii, jj, level, iRegion
        real(kind=realtype), dimension(:), allocatable :: extralocalbar, bcdatavaluesdlocal
        real(kind=realtype) :: dummyreal, dummyreald
        logical :: resetToRans
        real(kind=realtype), dimension(nSections) :: time
        logical :: useOldCoor ! for solverutils_b() functions
        useoldcoor = .false.
 
        ! extraLocalBar accumulates the seeds onto the extra variables
        allocate (extralocalbar(size(extrabar)))
        extralocalbar = zero
 
        ! Place the output spatial seed into the temporary petsc x-like
        ! vector.
        xbar = zero
        call vecplacearray(x_like, xbar, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Set the residual seeds.
        ii = 0
        do nn = 1, ndom
            do sps = 1, ntimeintervalsspectral
 
                ! Set pointers and derivative pointers
                call setpointers_b(nn, 1, sps)
 
                ! Set the dw seeds
                do k = 2, kl
                    do j = 2, jl
                        do i = 2, il
                            do l = 1, nstate
                                ii = ii + 1
                                dwd(i, j, k, l) = dwbar(ii)
                            end do
                        end do
                    end do
                end do
            end do
        end do
 
        forcespsloop: do sps = 1, ntimeintervalsspectral
            fsize = size(forcesbar, 2)
            call getforces_b(forcesbar(:, :, sps), fsize, sps)
        end do forcespsloop
 
        ! Call the final getSolution_b routine
        if (present(famlists)) then
            call getsolution_b(famlists, funcvalues, funcvaluesd)
        end if
 
        spsloop1: do sps = 1, ntimeintervalsspectral
 
            domainloop1: do nn = 1, ndom
                call setpointers_b(nn, 1, sps)
 
                ! Now we start running back through the main residual code:
                call resscale_b
                call sumdwandfw_b
 
                ! if (lowSpeedPreconditioner) then
                !    call applyLowSpeedPreconditioner_b
                ! end if
 
                ! Note that master_b does not include the first order flux
                ! approxation codes as those are never needed in reverse.
                if (viscous) then
                    call viscousflux_b
                    call allnodalgradients_b
                    call computespeedofsoundsquared_b
                end if
 
                ! So the all nodal gradients doesnt' perform the final
                ! scaling by the volume since it isn't necessary for the
                ! derivative. We have a special routine to fix that.
                if (viscous) then
                    call fixallnodalgradientsfromad()
                    call viscousflux
                end if
 
                select case (spacediscr)
                case (dissscalar)
                    call invisciddissfluxscalar_b
                case (dissmatrix)
                    call invisciddissfluxmatrix_b
                case (upwind)
                    call inviscidupwindflux_b(.true.)
                end select
 
                call inviscidcentralflux_b
                ! Compute turbulence residual for RANS equations
                if (equations == ransequations) then
                    select case (turbmodel)
                    case (spalartallmaras)
                        call saresscale_b
                        call saviscous_b
                        !call unsteadyTurbTerm_b(1_intType, 1_intType, itu1-1, qq)
                        call turbadvection_b(1_inttype, 1_inttype, itu1 - 1, qq)
                        ! turbAdvection_b zeros the faceid. This should be ok since
                        ! it presumably is the last call in master using faceid and
                        ! therefore should be the first call in master_b to use faceid
                        call sasource_b
                    end select
 
                    !call unsteadyTurbSpectral_block_b(itu1, itu1, nn, sps)
                end if
 
                call timestep_block_b(.false.)
 
                ! Just to be safe, zero the pLocald value...should not matter
                dummyreald = zero
                do iregion = 1, nactuatorregions
                    call sourceterms_block_b(nn, .true., iregion, dummyreal, dummyreald)
                end do
 
                call initres_block_b(1, nw, nn, sps)
            end do domainloop1
        end do spsloop1
 
        ! All reduce the global AZ volume seeds and store them in the
        ! local AZ volume seeds.
        ! This is the inverse of the all reduce that is done in forward
        ! mode to sum the local seeds into the global seeds.
        do iregion = 1, nactuatorregions
            call mpi_allreduce(actuatorregionsd(iregion)%volume, actuatorregionsd(iregion)%volLocal, 1, &
                               adflow_real, mpi_sum, adflow_comm_world, ierr)
            call echk(ierr, __file__, __line__)
        end do
 
        ! Need to re-apply the BCs. The reason is that BC halos behind
        ! interpolated cells need to be recomputed with their new
        ! interpolated values from actual compute cells. Only needed for
        ! overset.
        if (oversetpresent) then
            do sps = 1, ntimeintervalsspectral
                do nn = 1, ndom
                    call setpointers_b(nn, 1, sps)
                    call applyallbc_block_b(.true.)
                    call applyallbc_block(.true.)
 
                    if (equations == ransequations) then
                        call applyallturbbcthisblock_b(.true.)
                        call bcturbtreatment_b
                    end if
                end do
            end do
        end if
 
        ! Exchange the adjoint values.
        call whalo2_b(currentlevel, 1_inttype, nw, .true., .true., .true.)
 
        spsloop2: do sps = 1, ntimeintervalsspectral
 
            ! Get the pointers from the petsc vector for the wall
            ! surface and it's accumulation. Only necessary for wall
            ! distance.
            call vecgetarrayf90(xsurfvec(1, sps), xsurf, ierr)
            call echk(ierr, __file__, __line__)
 
            ! And it's derivative
            call vecgetarrayf90(xsurfvecd(sps), xsurfd, ierr)
            call echk(ierr, __file__, __line__)
 
            !Zero the accumulation vector on a per-time-spectral instance basis
            xsurfd = zero
 
            domainloop2: do nn = 1, ndom
                call setpointers_b(nn, 1, sps)
                call applyallbc_block_b(.true.)
 
                ! Run the forward application of the BCs. The reason is that
                ! the reverse application of the BCs can result in
                ! inconsisent values in the halos. This has only been
                ! observed to caluse problems with hot subsonic flow in an
                ! engine core. This is the same reason we need the applyBCs
                ! after the _b version above.
                call applyallbc_block(.true.)
 
                if (equations == ransequations) then
                    call applyallturbbcthisblock_b(.true.)
                    call bcturbtreatment_b
                end if
                call computeeddyviscosity_b(.false.)
                call computelamviscosity_b(.false.)
                call computepressuresimple_b(.false.)
 
                if (equations == ransequations .and. useapproxwalldistance) then
                    call updatewalldistancesquickly_b(nn, 1, sps)
                end if
 
                ! Here we insert the functions related to
                ! rotational (mesh movement) setup
                time = timeunsteadyrestart
                if (equationmode .eq. timespectral) then
                    do mm = 1, nsections
                        time(mm) = time(mm) + (sps - 1) * sections(mm)%timeperiod / real(&
                             &         ntimeintervalsspectral, realtype)
                    end do
                end if
 
                call slipvelocitiesfinelevel_block_b(useoldcoor, time, sps, nn)
                call normalvelocities_block_b(sps)
                call gridvelocitiesfinelevel_block_b(useoldcoor, time, sps, nn)
 
                call boundarynormals_b
                do iregion = 1, nactuatorregions
                    call computeactuatorregionvolume_b(nn, iregion)
                end do
                call metric_block_b
                call volume_block_b
 
            end do domainloop2
 
            ! Restore the petsc pointers.
            call vecgetarrayf90(xsurfvec(1, sps), xsurf, ierr)
            call echk(ierr, __file__, __line__)
 
            ! And it's derivative
            call vecgetarrayf90(xsurfvecd(sps), xsurfd, ierr)
            call echk(ierr, __file__, __line__)
 
            ! Now accumulate the xsurfd accumulation by using the wall scatter
            ! in reverse.
            if (walldistanceneeded .and. useapproxwalldistance) then
 
                call vecscatterbegin(wallscatter(1, sps), xsurfvecd(sps), x_like, add_values, scatter_reverse, ierr)
                call echk(ierr, __file__, __line__)
 
                call vecscatterend(wallscatter(1, sps), xsurfvecd(sps), x_like, add_values, scatter_reverse, ierr)
                call echk(ierr, __file__, __line__)
            end if
        end do spsloop2
 
        ! Zero out the local volume seeds of the actuator zone
        do iregion = 1, nactuatorregions
            actuatorregionsd(iregion)%volLocal = zero
        end do
 
        if (present(bcdatanames)) then
            allocate (bcdatavaluesdlocal(size(bcdatavaluesd)))
            bcdatavaluesdlocal = zero
            call setbcdatafinegrid_b(.true.)
            do sps = 1, ntimeintervalsspectral
                call setbcdata_b(bcdatanames, bcdatavalues, bcdatavaluesdlocal, bcdatafamlists, &
                                 sps, size(bcdatavalues), size(bcdatafamlists, 2))
            end do
            ! Reverse seeds need to accumulated across all processors:
            call mpi_allreduce(bcdatavaluesdlocal, bcdatavaluesd, size(bcdatavaluesd), adflow_real, &
                               mpi_sum, adflow_comm_world, ierr)
            deallocate (bcdatavaluesdlocal)
        end if
        call referencestate_b
        call adjustinflowangle_b
 
        do sps = 1, ntimeintervalsspectral
            ! Update overset connectivity if necessary
            if (oversetpresent) then
                if (oversetupdatemode == updatefast) then
                    call updateoversetconnectivity_b(1_inttype, sps)
                else if (oversetupdatemode == updatefull) then
                    if (myid == 0) then
                        print *, 'Full overset update derivatives not implemented'
                    end if
                end if
            end if
        end do
        ! Now the adjoint of the coordinate exhcange
        call exchangecoor_b(1)
        do nn = 1, ndom
            do sps = 1, ntimeintervalsspectral
                call setpointers_b(nn, 1, sps)
                call xhalo_block_b()
            end do
        end do
 
        call vecresetarray(x_like, ierr)
        call echk(ierr, __file__, __line__)
 
        ! =========================================
        !            End of reverse pass
        ! =========================================
 
        ! Store the extra derivatives
        extralocalbar(ialpha) = alphad
        extralocalbar(ibeta) = betad
        extralocalbar(imach) = machd + machcoefd
        extralocalbar(imachgrid) = machgridd
        extralocalbar(ipressure) = pinfdimd
        extralocalbar(itemperature) = tinfdimd
        extralocalbar(idensity) = rhoinfdimd
        extralocalbar(ipointrefx) = pointrefd(1)
        extralocalbar(ipointrefy) = pointrefd(2)
        extralocalbar(ipointrefz) = pointrefd(3)
 
        ! Finally put the output seeds into the provided vectors.
        ii = 0
        jj = 0
        do nn = 1, ndom
            do sps = 1, ntimeintervalsspectral
 
                ! Set pointers and derivative pointers
                call setpointers_b(nn, 1, sps)
                ! Set the wbar accumulation
                do k = 2, kl
                    do j = 2, jl
                        do i = 2, il
                            do l = 1, nstate
                                ii = ii + 1
                                wbar(ii) = wd(i, j, k, l)
                            end do
                        end do
                    end do
                end do
 
                ! Set the xvbar accumulation. Note that this must be a sum,
                ! becuase we may already have wall distance accumulation
                ! from the wallScatter directly into xbar (through x_like).
                do k = 1, kl
                    do j = 1, jl
                        do i = 1, il
                            do l = 1, 3
                                jj = jj + 1
                                xbar(jj) = xbar(jj) + xd(i, j, k, l)
                            end do
                        end do
                    end do
                end do
            end do
        end do
 
        ! Finally get the full contribution of the extra variables by
        ! summing all local contributions.
        extrabar = zero
        call mpi_allreduce(extralocalbar, extrabar, size(extrabar), adflow_real, &
                           mpi_sum, adflow_comm_world, ierr)
        call echk(ierr, __file__, __line__)
 
    end subroutine master_b
 
    subroutine master_state_b(wbar, dwBar, nState)
 
        ! This is specialized form of master that *ONLY* computes drdw
        ! products. It uses a few specialzed routines that are
        ! differentiated without including spatial dependencies. This
        ! results in slightly faster code. This specialization is
        ! justififed since this routine is needed for the transpose
        ! matrix-vector products in solving the adjoint system and thus
        ! this routine is called several orders of magniutde more than
        ! master_b. This routine has to be fast!
 
        use constants
        use iteration, only: currentlevel
        use flowvarrefstate, only: nw, viscous
        use blockpointers, only: ndom, il, jl, kl, wd, dwd, iblank
        use inputphysics, only: equationmode, turbmodel, equations
        use inputdiscretization, only: lowspeedpreconditioner, spacediscr
        use inputtimespectral, only: ntimeintervalsspectral
        use utils, only: setpointers_d
        use haloexchange, only: whalo2_b
        use flowutils, only: fixallnodalgradientsfromad
        use adjointextra_b, only: resscale_b, sumdwandfw_b
        use flowutils_b, only: computepressuresimple_b, computelamviscosity_b, &
                               computespeedofsoundsquared_b
        use turbbcroutines_b, only: applyallturbbcthisblock_b, bcturbtreatment_b
        use turbutils_b, only: computeeddyviscosity_b
        use bcextra_b, only: applyallbc_block_b
 
        use sa_fast_b, only: saresscale_fast_b, saviscous_fast_b, &
                             sasource_fast_b, qq
        use turbutils_fast_b, only: turbadvection_fast_b
        use fluxes_fast_b, only: inviscidupwindflux_fast_b, invisciddissfluxscalar_fast_b, &
                                 invisciddissfluxmatrix_fast_b, viscousflux_fast_b, inviscidcentralflux_fast_b
        use solverutils_fast_b, only: timestep_block_fast_b
        use flowutils_fast_b, only: allnodalgradients_fast_b
        use residuals_fast_b, only: sourceterms_block_fast_b, initres_block_fast_b
        use oversetdata, only: oversetpresent
        use bcroutines, only: applyallbc_block
        use actuatorregiondata, only: nactuatorregions
        implicit none
 
        ! Input variables:
        real(kind=realtype), intent(in), dimension(:) :: dwbar
        integer(kind=intType), intent(in) :: nState
 
        ! Input Arguments:
        real(kind=realtype), intent(out), dimension(:) :: wbar
 
        ! Working Variables
        integer(kind=intType) :: ierr, nn, sps, mm, i, j, k, l, fSize, ii, jj, level, iRegion
        real(kind=realtype) :: dummyreal
 
        ! Set the residual seeds.
        ii = 0
        do nn = 1, ndom
            do sps = 1, ntimeintervalsspectral
                call setpointers_d(nn, 1, sps)
                do k = 2, kl
                    do j = 2, jl
                        do i = 2, il
                            do l = 1, nstate
                                ii = ii + 1
                                dwd(i, j, k, l) = dwbar(ii)
                            end do
                        end do
                    end do
                end do
            end do
        end do
 
        ! ============================================
        !  reverse the order of calls from master
        ! ============================================
 
        spsloop1: do sps = 1, ntimeintervalsspectral
            domainloop1: do nn = 1, ndom
                call setpointers_d(nn, 1, sps)
 
                ! Now we start running back through the main residual code:
                call resscale_b
                call sumdwandfw_b
 
                ! if (lowSpeedPreconditioner) then
                !    call applyLowSpeedPreconditioner_b
                ! end if
 
                ! Note that master_b does not include the approximation codes
                ! as those are never needed in reverse.
                if (viscous) then
                    call viscousflux_fast_b
                    call allnodalgradients_fast_b
                    call computespeedofsoundsquared_b
                end if
 
                select case (spacediscr)
                case (dissscalar)
                    call invisciddissfluxscalar_fast_b
                case (dissmatrix)
                    call invisciddissfluxmatrix_fast_b
                case (upwind)
                    call inviscidupwindflux_fast_b(.true.)
                end select
 
                call inviscidcentralflux_fast_b
 
                ! Compute turbulence residual for RANS equations
                if (equations == ransequations) then
                    select case (turbmodel)
                    case (spalartallmaras)
                        call saresscale_fast_b
                        call saviscous_fast_b
                        !call unsteadyTurbTerm_b(1_intType, 1_intType, itu1-1, qq)
                        call turbadvection_fast_b(1_inttype, 1_inttype, itu1 - 1, qq)
                        call sasource_fast_b
                    end select
 
                    !call unsteadyTurbSpectral_block_b(itu1, itu1, nn, sps)
                end if
 
                call timestep_block_fast_b(.false.)
                do iregion = 1, nactuatorregions
                    call sourceterms_block_fast_b(nn, .true., iregion, dummyreal)
                end do
 
                call initres_block_fast_b(1, nw, nn, sps)
            end do domainloop1
        end do spsloop1
 
        ! Need to re-apply the BCs. The reason is that BC halos behind
        ! interpolated cells need to be recomputed with their new
        ! interpolated values from actual compute cells. Only needed for
        ! overset.
        if (oversetpresent) then
            do sps = 1, ntimeintervalsspectral
                do nn = 1, ndom
                    call setpointers_d(nn, 1, sps)
                    call applyallbc_block_b(.true.)
                    call applyallbc_block(.true.)
 
                    if (equations == ransequations) then
                        call applyallturbbcthisblock_b(.true.)
                        call bcturbtreatment_b
                    end if
                end do
            end do
        end if
 
        ! Exchange the adjoint values.
        call whalo2_b(currentlevel, 1_inttype, nw, .true., .true., .true.)
 
        spsloop2: do sps = 1, ntimeintervalsspectral
            domainloop2: do nn = 1, ndom
                call setpointers_d(nn, 1, sps)
 
                call applyallbc_block_b(.true.)
                call applyallbc_block(.true.)
 
                if (equations == ransequations) then
                    call applyallturbbcthisblock_b(.true.)
                    call bcturbtreatment_b
                end if
 
                call computeeddyviscosity_b(.false.)
                call computelamviscosity_b(.false.)
                call computepressuresimple_b(.false.)
 
            end do domainloop2
        end do spsloop2
 
        ! Finally put the output seeds into wbar
        ii = 0
        do nn = 1, ndom
            do sps = 1, ntimeintervalsspectral
                call setpointers_d(nn, 1, sps)
                do k = 2, kl
                    do j = 2, jl
                        do i = 2, il
                            do l = 1, nstate
                                ii = ii + 1
                                wbar(ii) = wd(i, j, k, l)!*max(real(iblank(i,j,k)), zero)
                            end do
                        end do
                    end do
                end do
            end do
        end do
    end subroutine master_state_b
#endif
    subroutine block_res_state(nn, sps, useFlowRes, useTurbRes)
 
        ! This is a special state-only routine used only for finite
        ! differce computations of the jacobian
        use constants
        use bcroutines, only: applyallbc_block
        use inputadjoint, only: viscpc
        use blockpointers, only: ndom, wd, xd, dw, il, jl, kl
        use flowvarrefstate, only: viscous
        use inputphysics, only: equations, turbmodel
        use inputdiscretization, only: lowspeedpreconditioner, lumpeddiss, spacediscr
        use utils, only: setpointers, echk
        use residuals, only: sourceterms_block
        use turbutils, only: computeeddyviscosity
        use turbbcroutines, only: applyallturbbcthisblock, bcturbtreatment
        use adjointextra, only: sumdwandfw, resscale
        use inputdiscretization, only: useblockettes
        use blockette, only: blocketterescore, blockrescore
        use flowutils, only: computelamviscosity, computepressuresimple
        use actuatorregiondata, only: nactuatorregions
        implicit none
 
        ! Input Arguments:
        integer(kind=intType), intent(in) :: nn, sps
        logical, optional, intent(in) :: useFlowRes, useTurbRes
 
        ! Working Variables
        integer(kind=intType) :: ierr, mm, i, j, k, l, fSize, ii, jj, iRegion
        real(kind=realtype) :: plocal
        logical :: dissApprox, viscApprox, updateIntermed, flowRes, turbRes, storeWall
 
        flowres = .true.
        if (present(useflowres)) then
            flowres = useflowres
        end if
 
        turbres = .true.
        if (present(useturbres)) then
            turbres = useturbres
        end if
 
        call computepressuresimple(.true.)
        call computelamviscosity(.true.)
        call computeeddyviscosity(.true.)
 
        ! Make sure to call the turb BC's first incase we need to
        ! correct for K
        if (equations == ransequations) then
            call bcturbtreatment
            call applyallturbbcthisblock(.true.)
        end if
 
        call applyallbc_block(.true.)
 
        ! Set the flags:
        dissapprox = lumpeddiss
        viscapprox = lumpeddiss
        updateintermed = .false.
        storewall = .true.
        blockettes: if (useblockettes) then
            call blocketterescore(dissapprox, viscapprox, updateintermed, flowres, turbres, storewall)
        else
            call blockrescore(dissapprox, viscapprox, updateintermed, flowres, turbres, storewall, nn, sps)
        end if blockettes
        do iregion = 1, nactuatorregions
            call sourceterms_block(nn, .true., iregion, plocal)
        end do
        call resscale
 
    end subroutine block_res_state
#ifndef USE_COMPLEX
    subroutine block_res_state_d(nn, sps)
 
        ! This is a special state-only forward mode linearization
        ! computation used to assemble the jacobian.
        use constants
        use bcextra_d, only: applyallbc_block_d
        use inputadjoint, only: viscpc
        use blockpointers, only: ndom, wd, xd, dw, dwd
        use flowvarrefstate, only: viscous
        use inputphysics, only: equations, turbmodel
        use inputdiscretization, only: lowspeedpreconditioner, lumpeddiss, spacediscr
        use inputtimespectral, only: ntimeintervalsspectral
        use utils, only: setpointers_d, echk
        use sa_d, only: sasource_d, saviscous_d, saresscale_d, qq
        use turbutils_d, only: turbadvection_d, computeeddyviscosity_d
        use fluxes_d, only: invisciddissfluxscalarapprox_d, invisciddissfluxmatrixapprox_d, &
                            inviscidupwindflux_d, invisciddissfluxscalar_d, invisciddissfluxmatrix_d, &
                            inviscidupwindflux_d, viscousflux_d, viscousfluxapprox_d, inviscidcentralflux_d
        use flowutils_d, only: computepressuresimple_d, computelamviscosity_d, &
                               computespeedofsoundsquared_d, allnodalgradients_d
        use solverutils_d, only: timestep_block_d
        use turbbcroutines_d, only: applyallturbbcthisblock_d, bcturbtreatment_d
        use adjointextra_d, only: resscale_d, sumdwandfw_d
        use residuals_d, only: sourceterms_block_d
        use actuatorregiondata, only: nactuatorregions
        implicit none
 
        ! Input Arguments:
        integer(kind=intType), intent(in) :: nn, sps
 
        ! Working Variables
        integer(kind=intType) :: ierr, mm, i, j, k, l, fSize, ii, jj, iRegion
        real(kind=realtype) :: dummyreal, dummyreald
 
        call computepressuresimple_d(.true.)
        call computelamviscosity_d(.true.)
        call computeeddyviscosity_d(.true.)
 
        ! Make sure to call the turb BC's first incase we need to
        ! correct for K
        if (equations == ransequations) then
            call bcturbtreatment_d
            call applyallturbbcthisblock_d(.true.)
        end if
 
        call applyallbc_block_d(.true.)
        call timestep_block_d(.false.)
 
        dw = zero ! These two lines are init_res
        dwd = zero
 
        ! Compute any source terms
        do iregion = 1, nactuatorregions
            call sourceterms_block_d(nn, .true., iregion, dummyreal, dummyreald)
        end do
 
        !Compute turbulence residual for RANS equations
        if (equations == ransequations) then
            !call unsteadyTurbSpectral_block(itu1, itu1, nn, sps)
 
            select case (turbmodel)
            case (spalartallmaras)
                call sasource_d
                call turbadvection_d(1_inttype, 1_inttype, itu1 - 1, qq)
          !!call unsteadyTurbTerm_d(1_intType, 1_intType, itu1-1, qq)
                call saviscous_d
                call saresscale_d
            end select
        end if
 
        ! compute the mean flow residual
        call inviscidcentralflux_d
 
        if (lumpeddiss) then
            select case (spacediscr)
            case (dissscalar)
                call invisciddissfluxscalarapprox_d
            case (dissmatrix)
                call invisciddissfluxmatrixapprox_d
            case (upwind)
                call inviscidupwindflux_d(.true.)
            end select
        else
            select case (spacediscr)
            case (dissscalar)
                call invisciddissfluxscalar_d
            case (dissmatrix)
                call invisciddissfluxmatrix_d
            case (upwind)
                call inviscidupwindflux_d(.true.)
            end select
        end if
 
        if (viscous) then
            call computespeedofsoundsquared_d
            if (.not. lumpeddiss .or. viscpc) then
                call allnodalgradients_d
                call viscousflux_d
            else
                call viscousfluxapprox_d
            end if
        end if
 
        ! if (lowSpeedPreconditioner) then
        !    call applyLowSpeedPreconditioner_d
        ! end if
        call sumdwandfw_d
        call resscale_d
    end subroutine block_res_state_d
#endif
 
end module masterroutines