surfaceIntegrations.F90

Go to the documentation of this file.

module surfaceintegrations
 
contains
 
    subroutine getcostfunctions(globalVals, funcValues)
 
        use constants
        use inputtimespectral, only: ntimeintervalsspectral
        use flowvarrefstate, only: pref, rhoref, tref, lref, gammainf, pinf, uref, uinf
        use inputphysics, only: liftdirection, dragdirection, surfaceref, &
                                machcoef, lengthref, alpha, beta, liftindex, cpmin_family, &
                                cpmin_rho, sepsenmaxfamily, sepsenmaxrho
        use inputcostfunctions, only: computecavitation, computesepsensorks, sepsensorkssharpness, &
                                      sepsensorksoffset
        use inputtsstabderiv, only: tsstability
        use utils, only: computetsderivatives
        use flowutils, only: getdirvector
        implicit none
 
        ! Input/Output
        real(kind=realtype), intent(in), dimension(:, :) :: globalvals
        real(Kind=realtype), intent(out), dimension(:) :: funcvalues
 
        ! Working
        real(kind=realtype) :: fact, factmoment, ovrnts
        real(kind=realtype), dimension(3, nTimeIntervalsSpectral) :: force, forcep, forcev, forcem, &
                                                                     moment, cforce, cforcep, cforcev, cforcem, &
                                                                     cmoment, cofx, cofy, cofz
        real(kind=realtype), dimension(3) :: vcoordref, vfreestreamref
        real(kind=realtype) :: mavgptot, mavgttot, mavgrho, mavgps, mflow, mavgmn, mavga, &
                               mavgvx, mavgvy, mavgvz, garea, mavgvi, fxlift, fylift, fzlift
 
        real(kind=realtype) :: vdotn, mag, u, v, w, ks_comp
        integer(kind=intType) :: sps
        real(kind=realtype), dimension(8) :: dcdq, dcdqdot
        real(kind=realtype), dimension(8) :: dcdalpha, dcdalphadot
        real(kind=realtype), dimension(8) :: coef0
 
        ! Factor used for time-averaged quantities.
        ovrnts = one / ntimeintervalsspectral
 
        ! Sum pressure and viscous contributions
        force = globalvals(ifp:ifp + 2, :) + globalvals(ifv:ifv + 2, :) + globalvals(iflowfm:iflowfm + 2, :)
        forcep = globalvals(ifp:ifp + 2, :)
        forcev = globalvals(ifv:ifv + 2, :)
        forcem = globalvals(iflowfm:iflowfm + 2, :)
 
        cofx = globalvals(icoforcex:icoforcex + 2, :)
        cofy = globalvals(icoforcey:icoforcey + 2, :)
        cofz = globalvals(icoforcez:icoforcez + 2, :)
 
        moment = globalvals(imp:imp + 2, :) + globalvals(imv:imv + 2, :) + globalvals(iflowmm:iflowmm + 2, :)
 
        fact = two / (gammainf * machcoef * machcoef &
                      * surfaceref * lref * lref * pref)
        cforce = fact * force
        cforcep = fact * forcep
        cforcev = fact * forcev
        cforcem = fact * forcem
 
        ! Moment factor has an extra lengthRef
        fact = fact / (lengthref * lref)
        cmoment = fact * moment
 
        ! Zero values since we are summing.
        funcvalues = zero
 
        ! Here we finally assign the final function values
        !$AD II-LOOP
        do sps = 1, ntimeintervalsspectral
            funcvalues(costfuncforcex) = funcvalues(costfuncforcex) + ovrnts * force(1, sps)
            funcvalues(costfuncforcey) = funcvalues(costfuncforcey) + ovrnts * force(2, sps)
            funcvalues(costfuncforcez) = funcvalues(costfuncforcez) + ovrnts * force(3, sps)
 
            funcvalues(costfuncforcexpressure) = funcvalues(costfuncforcexpressure) + ovrnts * forcep(1, sps)
            funcvalues(costfuncforceypressure) = funcvalues(costfuncforceypressure) + ovrnts * forcep(2, sps)
            funcvalues(costfuncforcezpressure) = funcvalues(costfuncforcezpressure) + ovrnts * forcep(3, sps)
 
            funcvalues(costfuncforcexviscous) = funcvalues(costfuncforcexviscous) + ovrnts * forcev(1, sps)
            funcvalues(costfuncforceyviscous) = funcvalues(costfuncforceyviscous) + ovrnts * forcev(2, sps)
            funcvalues(costfuncforcezviscous) = funcvalues(costfuncforcezviscous) + ovrnts * forcev(3, sps)
 
            funcvalues(costfuncforcexmomentum) = funcvalues(costfuncforcexmomentum) + ovrnts * forcem(1, sps)
            funcvalues(costfuncforceymomentum) = funcvalues(costfuncforceymomentum) + ovrnts * forcem(2, sps)
            funcvalues(costfuncforcezmomentum) = funcvalues(costfuncforcezmomentum) + ovrnts * forcem(3, sps)
 
            ! ------------
 
            funcvalues(costfuncforcexcoef) = funcvalues(costfuncforcexcoef) + ovrnts * cforce(1, sps)
            funcvalues(costfuncforceycoef) = funcvalues(costfuncforceycoef) + ovrnts * cforce(2, sps)
            funcvalues(costfuncforcezcoef) = funcvalues(costfuncforcezcoef) + ovrnts * cforce(3, sps)
 
            funcvalues(costfuncforcexcoefpressure) = funcvalues(costfuncforcexcoefpressure) + ovrnts * cforcep(1, sps)
            funcvalues(costfuncforceycoefpressure) = funcvalues(costfuncforceycoefpressure) + ovrnts * cforcep(2, sps)
            funcvalues(costfuncforcezcoefpressure) = funcvalues(costfuncforcezcoefpressure) + ovrnts * cforcep(3, sps)
 
            funcvalues(costfuncforcexcoefviscous) = funcvalues(costfuncforcexcoefviscous) + ovrnts * cforcev(1, sps)
            funcvalues(costfuncforceycoefviscous) = funcvalues(costfuncforceycoefviscous) + ovrnts * cforcev(2, sps)
            funcvalues(costfuncforcezcoefviscous) = funcvalues(costfuncforcezcoefviscous) + ovrnts * cforcev(3, sps)
 
            funcvalues(costfuncforcexcoefmomentum) = funcvalues(costfuncforcexcoefmomentum) + ovrnts * cforcem(1, sps)
            funcvalues(costfuncforceycoefmomentum) = funcvalues(costfuncforceycoefmomentum) + ovrnts * cforcem(2, sps)
            funcvalues(costfuncforcezcoefmomentum) = funcvalues(costfuncforcezcoefmomentum) + ovrnts * cforcem(3, sps)
 
            ! ------------
 
            ! center of pressure (these are actually center of all forces)
            ! protect the divisions against zero, and divide the weighed sum by the force magnitude
            ! for this time spectral instance before we add it to the sum
            if (force(1, sps) /= zero) then
                cofx(:, sps) = cofx(:, sps) / force(1, sps)
            else
                cofx(:, sps) = zero
            end if
 
            if (force(2, sps) /= zero) then
                cofy(:, sps) = cofy(:, sps) / force(2, sps)
            else
                cofy(:, sps) = zero
            end if
 
            if (force(3, sps) /= zero) then
                cofz(:, sps) = cofz(:, sps) / force(3, sps)
            else
                cofz(:, sps) = zero
            end if
 
            ! Fx
            funcvalues(costfunccoforcexx) = funcvalues(costfunccoforcexx) + ovrnts * cofx(1, sps)
            funcvalues(costfunccoforcexy) = funcvalues(costfunccoforcexy) + ovrnts * cofx(2, sps)
            funcvalues(costfunccoforcexz) = funcvalues(costfunccoforcexz) + ovrnts * cofx(3, sps)
 
            ! Fy
            funcvalues(costfunccoforceyx) = funcvalues(costfunccoforceyx) + ovrnts * cofy(1, sps)
            funcvalues(costfunccoforceyy) = funcvalues(costfunccoforceyy) + ovrnts * cofy(2, sps)
            funcvalues(costfunccoforceyz) = funcvalues(costfunccoforceyz) + ovrnts * cofy(3, sps)
 
            ! Fz
            funcvalues(costfunccoforcezx) = funcvalues(costfunccoforcezx) + ovrnts * cofz(1, sps)
            funcvalues(costfunccoforcezy) = funcvalues(costfunccoforcezy) + ovrnts * cofz(2, sps)
            funcvalues(costfunccoforcezz) = funcvalues(costfunccoforcezz) + ovrnts * cofz(3, sps)
 
            ! ------------
 
            funcvalues(costfuncmomx) = funcvalues(costfuncmomx) + ovrnts * moment(1, sps)
            funcvalues(costfuncmomy) = funcvalues(costfuncmomy) + ovrnts * moment(2, sps)
            funcvalues(costfuncmomz) = funcvalues(costfuncmomz) + ovrnts * moment(3, sps)
 
            funcvalues(costfuncmomxcoef) = funcvalues(costfuncmomxcoef) + ovrnts * cmoment(1, sps)
            funcvalues(costfuncmomycoef) = funcvalues(costfuncmomycoef) + ovrnts * cmoment(2, sps)
            funcvalues(costfuncmomzcoef) = funcvalues(costfuncmomzcoef) + ovrnts * cmoment(3, sps)
 
            ! final part of the KS computation
            if (computesepsensorks) then
                ! only calculate the log part if we are actually computing for separation for KS method.
                ks_comp = ovrnts * (sepsenmaxfamily(sps) + &
                                    log(globalvals(isepsensorks, sps)) / sepsenmaxrho)
 
                funcvalues(costfuncsepsensorks) = funcvalues(costfuncsepsensorks) + ks_comp
 
                funcvalues(costfuncsepsensorksarea) = funcvalues(costfuncsepsensorksarea) + &
                                                    ovrnts * globalvals(isepsensorksarea, sps) * ks_comp * &
                                                    one / (one + exp(2 * sepsensorkssharpness &
                                                                     * (ks_comp + sepsensorksoffset))) + &
                                                    ovrnts * globalvals(isepsensorarea, sps)
 
            end if
 
            funcvalues(costfuncsepsensor) = funcvalues(costfuncsepsensor) + ovrnts * globalvals(isepsensor, sps)
 
            funcvalues(costfunccavitation) = funcvalues(costfunccavitation) + ovrnts * globalvals(icavitation, sps)
            ! final part of the KS computation
            if (computecavitation) then
                ! only calculate the log part if we are actually computing for cavitation.
                ! If we are not computing cavitation, the iCpMin in globalVals will be zero,
                ! which doesn't play well with log. we just want to return zero here.
                funcvalues(costfunccpmin) = funcvalues(costfunccpmin) + ovrnts * &
                                            (cpmin_family(sps) - log(globalvals(icpmin, sps)) / cpmin_rho)
            end if
 
            funcvalues(costfuncaxismoment) = funcvalues(costfuncaxismoment) + ovrnts * globalvals(iaxismoment, sps)
            funcvalues(costfuncsepsensoravgx) = funcvalues(costfuncsepsensoravgx) + ovrnts * globalvals(isepavg, sps)
            funcvalues(costfuncsepsensoravgy) = funcvalues(costfuncsepsensoravgy) + &
                                                ovrnts * globalvals(isepavg + 1, sps)
            funcvalues(costfuncsepsensoravgz) = funcvalues(costfuncsepsensoravgz) + &
                                                ovrnts * globalvals(isepavg + 2, sps)
            funcvalues(costfuncarea) = funcvalues(costfuncarea) + ovrnts * globalvals(iarea, sps)
            funcvalues(costfuncflowpower) = funcvalues(costfuncflowpower) + ovrnts * globalvals(ipower, sps)
 
            funcvalues(costfunccperror2) = funcvalues(costfunccperror2) + ovrnts * globalvals(icperror2, sps)
 
            ! Mass flow like objective
            mflow = globalvals(imassflow, sps)
            if (mflow /= zero) then
                mavgptot = globalvals(imassptot, sps) / mflow
                mavgttot = globalvals(imassttot, sps) / mflow
                mavgrho = globalvals(imassrho, sps) / mflow
                mavgps = globalvals(imassps, sps) / mflow
                mavgmn = globalvals(imassmn, sps) / mflow
                mavga = globalvals(imassa, sps) / mflow
 
                mavgvx = globalvals(imassvx, sps) / mflow
                mavgvy = globalvals(imassvy, sps) / mflow
                mavgvz = globalvals(imassvz, sps) / mflow
 
                mavgvi = (globalvals(imassvi, sps) / mflow)
 
                mag = sqrt(globalvals(imassnx, sps)**2 + &
                           globalvals(imassny, sps)**2 + &
                           globalvals(imassnz, sps)**2)
 
            else
                mavgptot = zero
                mavgttot = zero
                mavgrho = zero
                mavgps = zero
                mavgmn = zero
                mavga = zero
                mavgvx = zero
                mavgvy = zero
                mavgvz = zero
                mavgvi = zero
 
            end if
 
            ! area averaged objectives
            garea = globalvals(iarea, sps)
            if (garea /= zero) then
                ! area averaged pressure
                funcvalues(costfuncaavgptot) = funcvalues(costfuncaavgptot) + &
                                               ovrnts * globalvals(iareaptot, sps) / garea
                funcvalues(costfuncaavgps) = funcvalues(costfuncaavgps) + &
                                             ovrnts * globalvals(iareaps, sps) / garea
            end if
 
            funcvalues(costfuncmdot) = funcvalues(costfuncmdot) + ovrnts * mflow
            funcvalues(costfuncmavgptot) = funcvalues(costfuncmavgptot) + ovrnts * mavgptot
            funcvalues(costfuncmavgttot) = funcvalues(costfuncmavgttot) + ovrnts * mavgttot
            funcvalues(costfuncmavgrho) = funcvalues(costfuncmavgrho) + ovrnts * mavgrho
            funcvalues(costfuncmavgps) = funcvalues(costfuncmavgps) + ovrnts * mavgps
            funcvalues(costfuncmavgmn) = funcvalues(costfuncmavgmn) + ovrnts * mavgmn
            funcvalues(costfuncmavga) = funcvalues(costfuncmavga) + ovrnts * mavga
            funcvalues(costfuncmavgvx) = funcvalues(costfuncmavgvx) + ovrnts * mavgvx
            funcvalues(costfuncmavgvy) = funcvalues(costfuncmavgvy) + ovrnts * mavgvy
            funcvalues(costfuncmavgvz) = funcvalues(costfuncmavgvz) + ovrnts * mavgvz
            funcvalues(costfuncmavgvi) = funcvalues(costfuncmavgvi) + ovrnts * mavgvi
            ! Bending moment calc - also broken.
            ! call computeRootBendingMoment(cForce, cMoment, liftIndex, bendingMoment)
            ! funcValues(costFuncBendingCoef) = funcValues(costFuncBendingCoef) + ovrNTS*bendingMoment
 
        end do
 
        ! Lift and Drag (coefficients): Dot product with the lift/drag direction.
        funcvalues(costfunclift) = &
            funcvalues(costfuncforcex) * liftdirection(1) + &
            funcvalues(costfuncforcey) * liftdirection(2) + &
            funcvalues(costfuncforcez) * liftdirection(3)
 
        funcvalues(costfuncliftpressure) = &
            funcvalues(costfuncforcexpressure) * liftdirection(1) + &
            funcvalues(costfuncforceypressure) * liftdirection(2) + &
            funcvalues(costfuncforcezpressure) * liftdirection(3)
 
        funcvalues(costfuncliftviscous) = &
            funcvalues(costfuncforcexviscous) * liftdirection(1) + &
            funcvalues(costfuncforceyviscous) * liftdirection(2) + &
            funcvalues(costfuncforcezviscous) * liftdirection(3)
 
        funcvalues(costfuncliftmomentum) = &
            funcvalues(costfuncforcexmomentum) * liftdirection(1) + &
            funcvalues(costfuncforceymomentum) * liftdirection(2) + &
            funcvalues(costfuncforcezmomentum) * liftdirection(3)
 
        !-----
 
        funcvalues(costfuncdrag) = &
            funcvalues(costfuncforcex) * dragdirection(1) + &
            funcvalues(costfuncforcey) * dragdirection(2) + &
            funcvalues(costfuncforcez) * dragdirection(3)
 
        funcvalues(costfuncdragpressure) = &
            funcvalues(costfuncforcexpressure) * dragdirection(1) + &
            funcvalues(costfuncforceypressure) * dragdirection(2) + &
            funcvalues(costfuncforcezpressure) * dragdirection(3)
 
        funcvalues(costfuncdragviscous) = &
            funcvalues(costfuncforcexviscous) * dragdirection(1) + &
            funcvalues(costfuncforceyviscous) * dragdirection(2) + &
            funcvalues(costfuncforcezviscous) * dragdirection(3)
 
        funcvalues(costfuncdragmomentum) = &
            funcvalues(costfuncforcexmomentum) * dragdirection(1) + &
            funcvalues(costfuncforceymomentum) * dragdirection(2) + &
            funcvalues(costfuncforcezmomentum) * dragdirection(3)
 
        !-----
 
        funcvalues(costfuncliftcoef) = &
            funcvalues(costfuncforcexcoef) * liftdirection(1) + &
            funcvalues(costfuncforceycoef) * liftdirection(2) + &
            funcvalues(costfuncforcezcoef) * liftdirection(3)
 
        funcvalues(costfuncliftcoefpressure) = &
            funcvalues(costfuncforcexcoefpressure) * liftdirection(1) + &
            funcvalues(costfuncforceycoefpressure) * liftdirection(2) + &
            funcvalues(costfuncforcezcoefpressure) * liftdirection(3)
 
        funcvalues(costfuncliftcoefviscous) = &
            funcvalues(costfuncforcexcoefviscous) * liftdirection(1) + &
            funcvalues(costfuncforceycoefviscous) * liftdirection(2) + &
            funcvalues(costfuncforcezcoefviscous) * liftdirection(3)
 
        funcvalues(costfuncliftcoefmomentum) = &
            funcvalues(costfuncforcexcoefmomentum) * liftdirection(1) + &
            funcvalues(costfuncforceycoefmomentum) * liftdirection(2) + &
            funcvalues(costfuncforcezcoefmomentum) * liftdirection(3)
 
        !-----
 
        funcvalues(costfuncdragcoef) = &
            funcvalues(costfuncforcexcoef) * dragdirection(1) + &
            funcvalues(costfuncforceycoef) * dragdirection(2) + &
            funcvalues(costfuncforcezcoef) * dragdirection(3)
 
        funcvalues(costfuncdragcoefpressure) = &
            funcvalues(costfuncforcexcoefpressure) * dragdirection(1) + &
            funcvalues(costfuncforceycoefpressure) * dragdirection(2) + &
            funcvalues(costfuncforcezcoefpressure) * dragdirection(3)
 
        funcvalues(costfuncdragcoefviscous) = &
            funcvalues(costfuncforcexcoefviscous) * dragdirection(1) + &
            funcvalues(costfuncforceycoefviscous) * dragdirection(2) + &
            funcvalues(costfuncforcezcoefviscous) * dragdirection(3)
 
        funcvalues(costfuncdragcoefmomentum) = &
            funcvalues(costfuncforcexcoefmomentum) * dragdirection(1) + &
            funcvalues(costfuncforceycoefmomentum) * dragdirection(2) + &
            funcvalues(costfuncforcezcoefmomentum) * dragdirection(3)
 
        ! ----- Center of Lift
 
        ! dot product the 3 forces with the lift direction separately
        fxlift = funcvalues(costfuncforcex) * liftdirection(1)
        fylift = funcvalues(costfuncforcey) * liftdirection(2)
        fzlift = funcvalues(costfuncforcez) * liftdirection(3)
 
        ! run the weighed average for the 3 components of center of lift
        ! protect against division by zero
        if ((fxlift + fylift + fzlift) /= zero) then
            funcvalues(costfunccofliftx) = (fxlift * funcvalues(costfunccoforcexx) + &
                                            fylift * funcvalues(costfunccoforceyx) + &
                                            fzlift * funcvalues(costfunccoforcezx)) / &
                                           (fxlift + fylift + fzlift)
            funcvalues(costfunccoflifty) = (fxlift * funcvalues(costfunccoforcexy) + &
                                            fylift * funcvalues(costfunccoforceyy) + &
                                            fzlift * funcvalues(costfunccoforcezy)) / &
                                           (fxlift + fylift + fzlift)
            funcvalues(costfunccofliftz) = (fxlift * funcvalues(costfunccoforcexz) + &
                                            fylift * funcvalues(costfunccoforceyz) + &
                                            fzlift * funcvalues(costfunccoforcezz)) / &
                                           (fxlift + fylift + fzlift)
        else
            funcvalues(costfunccofliftx) = zero
            funcvalues(costfunccoflifty) = zero
            funcvalues(costfunccofliftz) = zero
        end if
 
        ! -------------------- Time Spectral Objectives ------------------
 
        if (tsstability) then
            print *, 'Error: TSStabilityDerivatives are *BROKEN*. They need to be '&
                 &'completely verifed from scratch'
            stop
 
            call computetsderivatives(force, moment, coef0, dcdalpha, &
                                      dcdalphadot, dcdq, dcdqdot)
 
            funcvalues(costfunccl0) = coef0(1)
            funcvalues(costfunccd0) = coef0(2)
            funcvalues(costfunccfy0) = coef0(4)
            funcvalues(costfunccm0) = coef0(8)
 
            funcvalues(costfuncclalpha) = dcdalpha(1)
            funcvalues(costfunccdalpha) = dcdalpha(2)
            funcvalues(costfunccfyalpha) = dcdalpha(4)
            funcvalues(costfunccmzalpha) = dcdalpha(8)
 
            funcvalues(costfuncclalphadot) = dcdalphadot(1)
            funcvalues(costfunccdalphadot) = dcdalphadot(2)
            funcvalues(costfunccfyalphadot) = dcdalphadot(4)
            funcvalues(costfunccmzalphadot) = dcdalphadot(8)
 
            funcvalues(costfuncclq) = dcdq(1)
            funcvalues(costfunccdq) = dcdq(2)
            funcvalues(costfunccfyq) = dcdq(4)
            funcvalues(costfunccmzq) = dcdq(8)
 
            funcvalues(costfuncclqdot) = dcdqdot(1)
            funcvalues(costfunccdqdot) = dcdqdot(2)
            funcvalues(costfunccfyqdot) = dcdqdot(4)
            funcvalues(costfunccmzqdot) = dcdqdot(8)
        end if
 
    end subroutine getcostfunctions
 
    subroutine wallintegrationface(localValues, mm)
        !
        !       wallIntegrations computes the contribution of the block
        !       given by the pointers in blockPointers to the force and
        !       moment of the geometry. A distinction is made
        !       between the inviscid and viscous parts. In case the maximum
        !       yplus value must be monitored (only possible for rans), this
        !       value is also computed. The separation sensor and the cavita-
        !       tion sensor is also computed
        !       here.
        !
        use constants
        use communication
        use blockpointers
        use flowvarrefstate
        use inputcostfunctions
        use inputphysics, only: machcoef, pointref, veldirfreestream, &
                                equations, momentaxis, cpmin_family, cpmin_rho, &
                                cavitationnumber, sepsenmaxfamily, sepsenmaxrho
        use bcpointers
        implicit none
 
        ! Input/output variables
        real(kind=realtype), dimension(nLocalValues), intent(inout) :: localvalues
        integer(kind=intType) :: mm
 
        ! Local variables.
        real(kind=realtype), dimension(3) :: fp, fv, mp, mv
        real(kind=realtype), dimension(3) :: cofsumfx, cofsumfy, cofsumfz
        real(kind=realtype) :: yplusmax, sepsensorks, sepsensor, sepsensoravg(3), &
                               sepsensorarea, cavitation, cpmin_ks_sum, sepsensorksarea
        integer(kind=intType) :: i, j, ii, blk
 
        real(kind=realtype) :: pm1, fx, fy, fz, fn
        real(kind=realtype) :: vecttangential(3)
        real(kind=realtype) :: vectdotproductfsnormal
        real(kind=realtype) :: xc, xco, yc, yco, zc, zco, qf(3), r(3), n(3), l
        real(kind=realtype) :: fact, rho, mul, yplus, dwall
        real(kind=realtype) :: v(3), sensor, sensor1, cp, tmp, plocal, ks_exponent
        real(kind=realtype) :: tauxx, tauyy, tauzz
        real(kind=realtype) :: tauxy, tauxz, tauyz
 
        real(kind=realtype), dimension(3) :: refpoint
        real(kind=realtype), dimension(3, 2) :: axispoints
        real(kind=realtype) :: mx, my, mz, cellarea, m0x, m0y, m0z, mvaxis, mpaxis
        real(kind=realtype) :: cperror, cperror2
 
        select case (bcfaceid(mm))
        case (imin, jmin, kmin)
            fact = -one
        case (imax, jmax, kmax)
            fact = one
        end select
 
        ! Determine the reference point for the moment computation in
        ! meters.
 
        refpoint(1) = lref * pointref(1)
        refpoint(2) = lref * pointref(2)
        refpoint(3) = lref * pointref(3)
 
        ! Determine the points defining the axis about which to compute a moment
        axispoints(1, 1) = lref * momentaxis(1, 1)
        axispoints(1, 2) = lref * momentaxis(1, 2)
        axispoints(2, 1) = lref * momentaxis(2, 1)
        axispoints(2, 2) = lref * momentaxis(2, 2)
        axispoints(3, 1) = lref * momentaxis(3, 1)
        axispoints(3, 2) = lref * momentaxis(3, 2)
 
        ! Initialize the force and moment coefficients to 0 as well as
        ! yplusMax.
 
        fp = zero; fv = zero;
        mp = zero; mv = zero;
        cofsumfx = zero; cofsumfy = zero; cofsumfz = zero
        yplusmax = zero
        sepsensor = zero
        sepsensorks = zero
        sepsensorarea = zero
        sepsensorksarea = zero
        cavitation = zero
        cpmin_ks_sum = zero
        sepsensoravg = zero
        mpaxis = zero; mvaxis = zero;
        cperror2 = zero;
        !
        !         Integrate the inviscid contribution over the solid walls,
        !         either inviscid or viscous. The integration is done with
        !         cp. For closed contours this is equal to the integration
        !         of p; for open contours this is not the case anymore.
        !         Question is whether a force for an open contour is
        !         meaningful anyway.
        !
 
        ! Loop over the quadrilateral faces of the subface. Note that
        ! the nodal range of BCData must be used and not the cell
        ! range, because the latter may include the halo's in i and
        ! j-direction. The offset +1 is there, because inBeg and jnBeg
        ! refer to nodal ranges and not to cell ranges. The loop
        ! (without the AD stuff) would look like:
        !
        ! do j=(BCData(mm)%jnBeg+1),BCData(mm)%jnEnd
        !    do i=(BCData(mm)%inBeg+1),BCData(mm)%inEnd
 
        !$AD II-LOOP
        do ii = 0, (bcdata(mm)%jnEnd - bcdata(mm)%jnBeg) * (bcdata(mm)%inEnd - bcdata(mm)%inBeg) - 1
            i = mod(ii, (bcdata(mm)%inEnd - bcdata(mm)%inBeg)) + bcdata(mm)%inBeg + 1
            j = ii / (bcdata(mm)%inEnd - bcdata(mm)%inBeg) + bcdata(mm)%jnBeg + 1
 
            ! Compute the average pressure minus 1 and the coordinates
            ! of the centroid of the face relative from from the
            ! moment reference point. Due to the usage of pointers for
            ! the coordinates, whose original array starts at 0, an
            ! offset of 1 must be used. The pressure is multipled by
            ! fact to account for the possibility of an inward or
            ! outward pointing normal.
 
            pm1 = fact * (half * (pp2(i, j) + pp1(i, j)) - pinf) * pref
 
            tmp = two / (gammainf * pinf * machcoef * machcoef)
            cp = (half * (pp2(i, j) + pp1(i, j)) - pinf) * tmp
            cperror = (cp - bcdata(mm)%CpTarget(i, j))
            cperror2 = cperror2 + cperror * cperror
 
            xc = fourth * (xx(i, j, 1) + xx(i + 1, j, 1) &
                           + xx(i, j + 1, 1) + xx(i + 1, j + 1, 1)) - refpoint(1)
            yc = fourth * (xx(i, j, 2) + xx(i + 1, j, 2) &
                           + xx(i, j + 1, 2) + xx(i + 1, j + 1, 2)) - refpoint(2)
            zc = fourth * (xx(i, j, 3) + xx(i + 1, j, 3) &
                           + xx(i, j + 1, 3) + xx(i + 1, j + 1, 3)) - refpoint(3)
 
            ! Compute the force components.
            blk = max(bcdata(mm)%iblank(i, j), 0)
            fx = pm1 * ssi(i, j, 1)
            fy = pm1 * ssi(i, j, 2)
            fz = pm1 * ssi(i, j, 3)
 
            ! Note from AY: Technically, we can just compute the moments using the center of force
            ! terms. However, the moment computations coded here distinguish pressure,
            ! viscous, and momentum contributions to moment. Even though these individual
            ! contributions are not exposed to python, I still wanted to keep how it's done in the
            ! code in case its useful in the future. This is also true for the face integrations
 
            ! Update the inviscid force and moment coefficients. Iblank as we sum
            fp(1) = fp(1) + fx * blk
            fp(2) = fp(2) + fy * blk
            fp(3) = fp(3) + fz * blk
 
            mx = yc * fz - zc * fy
            my = zc * fx - xc * fz
            mz = xc * fy - yc * fx
 
            mp(1) = mp(1) + mx * blk
            mp(2) = mp(2) + my * blk
            mp(3) = mp(3) + mz * blk
 
            ! the force integral for the center of pressure computation.
            ! We need the cell centers wrt origin
            xco = fourth * (xx(i, j, 1) + xx(i + 1, j, 1) &
                            + xx(i, j + 1, 1) + xx(i + 1, j + 1, 1))
            yco = fourth * (xx(i, j, 2) + xx(i + 1, j, 2) &
                            + xx(i, j + 1, 2) + xx(i + 1, j + 1, 2))
            zco = fourth * (xx(i, j, 3) + xx(i + 1, j, 3) &
                            + xx(i, j + 1, 3) + xx(i + 1, j + 1, 3))
 
            ! accumulate in the sums. each force component is tracked separately
 
            ! Force-X
            cofsumfx(1) = cofsumfx(1) + xco * fx * blk
            cofsumfx(2) = cofsumfx(2) + yco * fx * blk
            cofsumfx(3) = cofsumfx(3) + zco * fx * blk
 
            ! Force-Y
            cofsumfy(1) = cofsumfy(1) + xco * fy * blk
            cofsumfy(2) = cofsumfy(2) + yco * fy * blk
            cofsumfy(3) = cofsumfy(3) + zco * fy * blk
 
            ! Force-Z
            cofsumfz(1) = cofsumfz(1) + xco * fz * blk
            cofsumfz(2) = cofsumfz(2) + yco * fz * blk
            cofsumfz(3) = cofsumfz(3) + zco * fz * blk
 
            ! Compute the r and n vectores for the moment around an
            ! axis computation where r is the distance from the
            ! force to the first point on the axis and n is a unit
            ! normal in the direction of the axis
            r(1) = fourth * (xx(i, j, 1) + xx(i + 1, j, 1) &
                             + xx(i, j + 1, 1) + xx(i + 1, j + 1, 1)) - axispoints(1, 1)
            r(2) = fourth * (xx(i, j, 2) + xx(i + 1, j, 2) &
                             + xx(i, j + 1, 2) + xx(i + 1, j + 1, 2)) - axispoints(2, 1)
            r(3) = fourth * (xx(i, j, 3) + xx(i + 1, j, 3) &
                             + xx(i, j + 1, 3) + xx(i + 1, j + 1, 3)) - axispoints(3, 1)
 
            l = sqrt((axispoints(1, 2) - axispoints(1, 1))**2 &
                     + (axispoints(2, 2) - axispoints(2, 1))**2 &
                     + (axispoints(3, 2) - axispoints(3, 1))**2)
 
            n(1) = (axispoints(1, 2) - axispoints(1, 1)) / l
            n(2) = (axispoints(2, 2) - axispoints(2, 1)) / l
            n(3) = (axispoints(3, 2) - axispoints(3, 1)) / l
 
            ! Compute the moment of the force about the first point
            ! used to define the axis, and the project that axis in
            ! the n direction
            m0x = r(2) * fz - r(3) * fy
            m0y = r(3) * fx - r(1) * fz
            m0z = r(1) * fy - r(2) * fx
            mpaxis = mpaxis + (m0x * n(1) + m0y * n(2) + m0z * n(3)) * blk
 
            ! Save the face-based forces and area
            bcdata(mm)%Fp(i, j, 1) = fx
            bcdata(mm)%Fp(i, j, 2) = fy
            bcdata(mm)%Fp(i, j, 3) = fz
            cellarea = sqrt(ssi(i, j, 1)**2 + ssi(i, j, 2)**2 + ssi(i, j, 3)**2)
 
            bcdata(mm)%area(i, j) = cellarea
 
            ! Get normalized surface velocity:
            v(1) = ww2(i, j, ivx)
            v(2) = ww2(i, j, ivy)
            v(3) = ww2(i, j, ivz)
            v = v / (sqrt(v(1)**2 + v(2)**2 + v(3)**2) + 1e-16)
 
            if (computesepsensorks) then
                ! Freestream projection over the surface.
                vectdotproductfsnormal = veldirfreestream(1) * bcdata(mm)%norm(i, j, 1) + &
                                         veldirfreestream(2) * bcdata(mm)%norm(i, j, 2) + &
                                         veldirfreestream(3) * bcdata(mm)%norm(i, j, 3)
                ! Tangential Vector on the surface, which is the freestream projected vector
                vecttangential(1) = veldirfreestream(1) - vectdotproductfsnormal * bcdata(mm)%norm(i, j, 1)
                vecttangential(2) = veldirfreestream(2) - vectdotproductfsnormal * bcdata(mm)%norm(i, j, 2)
                vecttangential(3) = veldirfreestream(3) - vectdotproductfsnormal * bcdata(mm)%norm(i, j, 3)
 
                vecttangential = vecttangential / (sqrt(vecttangential(1)**2 + vecttangential(2)**2 + &
                                                        vecttangential(3)**2) + 1e-16)
 
                ! computing separation sensor
                ! velocity dot products
                sensor = (v(1) * vecttangential(1) + v(2) * vecttangential(2) + &
                          v(3) * vecttangential(3))
 
                ! sepsensor value
                sensor = (cos(degtorad * sepsensorksphi) - sensor) / &
                         (-cos(degtorad * sepsensorksphi) + cos(zero) + 1e-16)
 
                ! also do the ks-based spensenor max computation
                call ksaggregationfunction(sensor, sepsenmaxfamily(spectralsol), sepsenmaxrho, ks_exponent)
 
                sepsensorksarea = sepsensorksarea + blk * cellarea
                sepsensorarea = cellarea * blk * one / (one + exp(-2 * sepsensorkssharpness &
                                                                  * (sensor + sepsensorksoffset))) + sepsensorarea
 
                sepsensorks = sepsensorks + ks_exponent * blk
 
            end if
 
            ! Dot product with free stream
            sensor = -(v(1) * veldirfreestream(1) + v(2) * veldirfreestream(2) + &
                       v(3) * veldirfreestream(3))
 
            !Now run through a smooth heaviside function:
            sensor = one / (one + exp(-2 * sepsensorsharpness * (sensor - sepsensoroffset)))
 
            ! And integrate over the area of this cell and save, blanking as we go.
            sensor = sensor * cellarea * blk
            sepsensor = sepsensor + sensor
 
            ! Also accumulate into the sepSensorAvg
            ! x-y-zco are computed above for center of force computations
            sepsensoravg(1) = sepsensoravg(1) + sensor * xco
            sepsensoravg(2) = sepsensoravg(2) + sensor * yco
            sepsensoravg(3) = sepsensoravg(3) + sensor * zco
 
            if (computecavitation) then
                plocal = pp2(i, j)
                tmp = two / (gammainf * machcoef * machcoef)
                cp = tmp * (plocal - pinf)
                sensor1 = -cp - cavitationnumber
                sensor1 = (sensor1**cavexponent) / (one + exp(2 * cavsensorsharpness * (-sensor1 + cavsensoroffset)))
                sensor1 = sensor1 * cellarea * blk
                cavitation = cavitation + sensor1
 
                ! also do the ks-based cpmin computation
                ks_exponent = exp(cpmin_rho * (-cp + cpmin_family(spectralsol)))
                cpmin_ks_sum = cpmin_ks_sum + ks_exponent * blk
            end if
        end do
 
        !
        ! Integration of the viscous forces.
        ! Only for viscous boundaries.
        !
        visforce: if (bctype(mm) == nswalladiabatic .or. &
                      bctype(mm) == nswallisothermal) then
 
            ! Initialize dwall for the laminar case and set the pointer
            ! for the unit normals.
 
            dwall = zero
 
            ! Loop over the quadrilateral faces of the subface and
            ! compute the viscous contribution to the force and
            ! moment and update the maximum value of y+.
 
            !$AD II-LOOP
            do ii = 0, (bcdata(mm)%jnEnd - bcdata(mm)%jnBeg) * (bcdata(mm)%inEnd - bcdata(mm)%inBeg) - 1
                i = mod(ii, (bcdata(mm)%inEnd - bcdata(mm)%inBeg)) + bcdata(mm)%inBeg + 1
                j = ii / (bcdata(mm)%inEnd - bcdata(mm)%inBeg) + bcdata(mm)%jnBeg + 1
 
                ! Store the viscous stress tensor a bit easier.
                blk = max(bcdata(mm)%iblank(i, j), 0)
 
                tauxx = viscsubface(mm)%tau(i, j, 1)
                tauyy = viscsubface(mm)%tau(i, j, 2)
                tauzz = viscsubface(mm)%tau(i, j, 3)
                tauxy = viscsubface(mm)%tau(i, j, 4)
                tauxz = viscsubface(mm)%tau(i, j, 5)
                tauyz = viscsubface(mm)%tau(i, j, 6)
 
                ! Compute the viscous force on the face. A minus sign
                ! is now present, due to the definition of this force.
 
                fx = -fact * (tauxx * ssi(i, j, 1) + tauxy * ssi(i, j, 2) &
                              + tauxz * ssi(i, j, 3)) * pref
                fy = -fact * (tauxy * ssi(i, j, 1) + tauyy * ssi(i, j, 2) &
                              + tauyz * ssi(i, j, 3)) * pref
                fz = -fact * (tauxz * ssi(i, j, 1) + tauyz * ssi(i, j, 2) &
                              + tauzz * ssi(i, j, 3)) * pref
 
                ! Compute the coordinates of the centroid of the face
                ! relative from the moment reference point. Due to the
                ! usage of pointers for xx and offset of 1 is present,
                ! because x originally starts at 0.
 
                xc = fourth * (xx(i, j, 1) + xx(i + 1, j, 1) &
                               + xx(i, j + 1, 1) + xx(i + 1, j + 1, 1)) - refpoint(1)
                yc = fourth * (xx(i, j, 2) + xx(i + 1, j, 2) &
                               + xx(i, j + 1, 2) + xx(i + 1, j + 1, 2)) - refpoint(2)
                zc = fourth * (xx(i, j, 3) + xx(i + 1, j, 3) &
                               + xx(i, j + 1, 3) + xx(i + 1, j + 1, 3)) - refpoint(3)
 
                ! Update the viscous force and moment coefficients, blanking as we go.
 
                fv(1) = fv(1) + fx * blk
                fv(2) = fv(2) + fy * blk
                fv(3) = fv(3) + fz * blk
 
                mx = yc * fz - zc * fy
                my = zc * fx - xc * fz
                mz = xc * fy - yc * fx
 
                mv(1) = mv(1) + mx * blk
                mv(2) = mv(2) + my * blk
                mv(3) = mv(3) + mz * blk
 
                ! the force integral for the center of pressure computation.
                ! We need the cell centers wrt origin
                xco = fourth * (xx(i, j, 1) + xx(i + 1, j, 1) &
                                + xx(i, j + 1, 1) + xx(i + 1, j + 1, 1))
                yco = fourth * (xx(i, j, 2) + xx(i + 1, j, 2) &
                                + xx(i, j + 1, 2) + xx(i + 1, j + 1, 2))
                zco = fourth * (xx(i, j, 3) + xx(i + 1, j, 3) &
                                + xx(i, j + 1, 3) + xx(i + 1, j + 1, 3))
 
                ! accumulate in the sums. each force component is tracked separately
 
                ! Force-X
                cofsumfx(1) = cofsumfx(1) + xco * fx * blk
                cofsumfx(2) = cofsumfx(2) + yco * fx * blk
                cofsumfx(3) = cofsumfx(3) + zco * fx * blk
 
                ! Force-Y
                cofsumfy(1) = cofsumfy(1) + xco * fy * blk
                cofsumfy(2) = cofsumfy(2) + yco * fy * blk
                cofsumfy(3) = cofsumfy(3) + zco * fy * blk
 
                ! Force-Z
                cofsumfz(1) = cofsumfz(1) + xco * fz * blk
                cofsumfz(2) = cofsumfz(2) + yco * fz * blk
                cofsumfz(3) = cofsumfz(3) + zco * fz * blk
 
                ! Compute the r and n vectors for the moment around an
                ! axis computation where r is the distance from the
                ! force to the first point on the axis and n is a unit
                ! normal in the direction of the axis
                r(1) = fourth * (xx(i, j, 1) + xx(i + 1, j, 1) &
                                 + xx(i, j + 1, 1) + xx(i + 1, j + 1, 1)) - axispoints(1, 1)
                r(2) = fourth * (xx(i, j, 2) + xx(i + 1, j, 2) &
                                 + xx(i, j + 1, 2) + xx(i + 1, j + 1, 2)) - axispoints(2, 1)
                r(3) = fourth * (xx(i, j, 3) + xx(i + 1, j, 3) &
                                 + xx(i, j + 1, 3) + xx(i + 1, j + 1, 3)) - axispoints(3, 1)
 
                l = sqrt((axispoints(1, 2) - axispoints(1, 1))**2 &
                         + (axispoints(2, 2) - axispoints(2, 1))**2 &
                         + (axispoints(3, 2) - axispoints(3, 1))**2)
 
                n(1) = (axispoints(1, 2) - axispoints(1, 1)) / l
                n(2) = (axispoints(2, 2) - axispoints(2, 1)) / l
                n(3) = (axispoints(3, 2) - axispoints(3, 1)) / l
 
                ! Compute the moment of the force about the first point
                ! used to define the axis, and then project that axis in
                ! the n direction
                m0x = r(2) * fz - r(3) * fy
                m0y = r(3) * fx - r(1) * fz
                m0z = r(1) * fy - r(2) * fx
                mvaxis = mvaxis + (m0x * n(1) + m0y * n(2) + m0z * n(3)) * blk
 
                ! Save the face based forces for the slice operations
                bcdata(mm)%Fv(i, j, 1) = fx
                bcdata(mm)%Fv(i, j, 2) = fy
                bcdata(mm)%Fv(i, j, 3) = fz
 
                ! Compute the tangential component of the stress tensor,
                ! which is needed to monitor y+. The result is stored
                ! in fx, fy, fz, although it is not really a force.
                ! As later on only the magnitude of the tangential
                ! component is important, there is no need to take the
                ! sign into account (it should be a minus sign).
 
                fx = tauxx * bcdata(mm)%norm(i, j, 1) + tauxy * bcdata(mm)%norm(i, j, 2) &
                     + tauxz * bcdata(mm)%norm(i, j, 3)
                fy = tauxy * bcdata(mm)%norm(i, j, 1) + tauyy * bcdata(mm)%norm(i, j, 2) &
                     + tauyz * bcdata(mm)%norm(i, j, 3)
                fz = tauxz * bcdata(mm)%norm(i, j, 1) + tauyz * bcdata(mm)%norm(i, j, 2) &
                     + tauzz * bcdata(mm)%norm(i, j, 3)
 
                fn = fx * bcdata(mm)%norm(i, j, 1) + fy * bcdata(mm)%norm(i, j, 2) + fz * bcdata(mm)%norm(i, j, 3)
 
                fx = fx - fn * bcdata(mm)%norm(i, j, 1)
                fy = fy - fn * bcdata(mm)%norm(i, j, 2)
                fz = fz - fn * bcdata(mm)%norm(i, j, 3)
 
                ! Compute the local value of y+. Due to the usage
                ! of pointers there is on offset of -1 in dd2Wall..
#ifndef USE_TAPENADE
                if (equations == ransequations) then
                    dwall = dd2wall(i - 1, j - 1)
                    rho = half * (ww2(i, j, irho) + ww1(i, j, irho))
                    mul = half * (rlv2(i, j) + rlv1(i, j))
                    yplus = sqrt(rho * sqrt(fx * fx + fy * fy + fz * fz)) * dwall / mul
 
                    ! Store this value if this value is larger than the
                    ! currently stored value. Blank non-active cells.
 
                    yplusmax = max(yplusmax, yplus * blk)
                end if
#endif
            end do
        else
            ! If we had no viscous force, set the viscous component to zero
            bcdata(mm)%Fv = zero
        end if visforce
 
        ! Increment the local values array with the values we computed here.
        localvalues(ifp:ifp + 2) = localvalues(ifp:ifp + 2) + fp
        localvalues(ifv:ifv + 2) = localvalues(ifv:ifv + 2) + fv
        localvalues(imp:imp + 2) = localvalues(imp:imp + 2) + mp
        localvalues(imv:imv + 2) = localvalues(imv:imv + 2) + mv
        localvalues(icoforcex:icoforcex + 2) = localvalues(icoforcex:icoforcex + 2) + cofsumfx
        localvalues(icoforcey:icoforcey + 2) = localvalues(icoforcey:icoforcey + 2) + cofsumfy
        localvalues(icoforcez:icoforcez + 2) = localvalues(icoforcez:icoforcez + 2) + cofsumfz
        localvalues(isepsensor) = localvalues(isepsensor) + sepsensor
        localvalues(isepsensorks) = localvalues(isepsensorks) + sepsensorks
        localvalues(isepsensorksarea) = localvalues(isepsensorksarea) + sepsensorksarea
        localvalues(isepsensorarea) = localvalues(isepsensorarea) + sepsensorarea
        localvalues(icavitation) = localvalues(icavitation) + cavitation
        localvalues(icpmin) = localvalues(icpmin) + cpmin_ks_sum
        localvalues(isepavg:isepavg + 2) = localvalues(isepavg:isepavg + 2) + sepsensoravg
        localvalues(iaxismoment) = localvalues(iaxismoment) + mpaxis + mvaxis
        localvalues(icperror2) = localvalues(icperror2) + cperror2
 
#ifndef USE_TAPENADE
        localvalues(iyplus) = max(localvalues(iyplus), yplusmax)
#endif
    end subroutine wallintegrationface
 
    subroutine ksaggregationfunction(g, max_g, g_rho, ks_g)
 
        use precision
 
        real(kind=realtype), intent(inout) :: ks_g
        real(kind=realtype), intent(in) :: g, max_g, g_rho
 
        ks_g = exp(g_rho * (g - max_g))
 
    end subroutine ksaggregationfunction
 
    subroutine flowintegrationface(isInflow, localValues, mm)
 
        use constants
        use blockpointers, only: bctype, bcfaceid, bcdata, addgridvelocities
        use flowvarrefstate, only: pref, pinf, rhoref, timeref, lref, tref, rgas, uref, uinf, rhoinf, gammainf
        use inputphysics, only: pointref, flowtype, rgasdim
        use flowutils, only: computeptot, computettot
        use bcpointers, only: ssi, sface, ww1, ww2, pp1, pp2, xx, gamma1, gamma2
        use utils, only: mynorm2
        implicit none
 
        ! Input/Output variables
        logical, intent(in) :: isInflow
        real(kind=realtype), dimension(nLocalValues), intent(inout) :: localvalues
        integer(kind=intType), intent(in) :: mm
 
        ! Local variables
        real(kind=realtype) :: massflowrate, mass_ptot, mass_ttot, mass_ps, mass_mn, mass_a, mass_rho, &
                               mass_vx, mass_vy, mass_vz, mass_nx, mass_ny, mass_nz, mass_vi
        real(kind=realtype) :: area_ptot, area_ps
        real(kind=realtype) :: govgm1, gm1ovg, viconst, vilocal, pratio
        real(kind=realtype) :: mredim
        integer(kind=intType) :: i, j, ii, blk
        real(kind=realtype) :: internalflowfact, inflowfact, fact, xc, xco, yc, yco, zc, zco, mx, my, mz
        real(kind=realtype) :: sf, vmag, vnm, vnmfreestreamref, vxm, vym, vzm, fx, fy, fz, u, v, w
        real(kind=realtype) :: pm, ptot, ttot, rhom, gammam, am
        real(kind=realtype) :: area, cellarea, overcellarea
        real(kind=realtype), dimension(3) :: fp, mp, fmom, mmom, refpoint, sfacecoordref
        real(kind=realtype), dimension(3) :: cofsumfx, cofsumfy, cofsumfz
        real(kind=realtype) :: mnm, massflowratelocal
 
        refpoint(1) = lref * pointref(1)
        refpoint(2) = lref * pointref(2)
        refpoint(3) = lref * pointref(3)
 
        ! Note that these are *opposite* of force integrations. The reason
        ! is that we want positive mass flow into the domain and negative
        ! mass flow out of the domain. Since the low faces have ssi
        ! vectors pointining into the domain, this is correct. The high
        ! end faces need to flip this.
        select case (bcfaceid(mm))
        case (imin, jmin, kmin)
            fact = one
        case (imax, jmax, kmax)
            fact = -one
        end select
 
        ! the sign of momentum forces are flipped for internal flows
        internalflowfact = one
        if (flowtype == internalflow) then
            internalflowfact = -one
        end if
 
        inflowfact = one
        if (isinflow) then
            inflowfact = -one
        end if
 
        ! Loop over the quadrilateral faces of the subface. Note that
        ! the nodal range of BCData must be used and not the cell
        ! range, because the latter may include the halo's in i and
        ! j-direction. The offset +1 is there, because inBeg and jnBeg
        ! refer to nodal ranges and not to cell ranges. The loop
        ! (without the AD stuff) would look like:
        !
        ! do j=(BCData(mm)%jnBeg+1),BCData(mm)%jnEnd
        !    do i=(BCData(mm)%inBeg+1),BCData(mm)%inEnd
 
        mredim = sqrt(pref * rhoref)
        fp = zero
        mp = zero
        fmom = zero
        mmom = zero
 
        cofsumfx = zero
        cofsumfy = zero
        cofsumfz = zero
 
        massflowrate = zero
        area = zero
        mass_ptot = zero
        mass_ttot = zero
        mass_ps = zero
        mass_mn = zero
        mass_a = zero
        mass_rho = zero
 
        mass_vx = zero
        mass_vy = zero
        mass_vz = zero
        mass_nx = zero
        mass_ny = zero
        mass_nz = zero
        mass_vi = zero
 
        area_ptot = zero
        area_ps = zero
 
        !$AD II-LOOP
        do ii = 0, (bcdata(mm)%jnEnd - bcdata(mm)%jnBeg) * (bcdata(mm)%inEnd - bcdata(mm)%inBeg) - 1
            i = mod(ii, (bcdata(mm)%inEnd - bcdata(mm)%inBeg)) + bcdata(mm)%inBeg + 1
            j = ii / (bcdata(mm)%inEnd - bcdata(mm)%inBeg) + bcdata(mm)%jnBeg + 1
 
            if (addgridvelocities) then
                sf = sface(i, j)
            else
                sf = zero
            end if
 
            ! Compute the force components.
            blk = max(bcdata(mm)%iblank(i, j), 0) ! iBlank forces for overset stuff
 
            vxm = half * (ww1(i, j, ivx) + ww2(i, j, ivx))
            vym = half * (ww1(i, j, ivy) + ww2(i, j, ivy))
            vzm = half * (ww1(i, j, ivz) + ww2(i, j, ivz))
            rhom = half * (ww1(i, j, irho) + ww2(i, j, irho))
            pm = half * (pp1(i, j) + pp2(i, j))
            gammam = half * (gamma1(i, j) + gamma2(i, j))
 
            vnm = vxm * ssi(i, j, 1) + vym * ssi(i, j, 2) + vzm * ssi(i, j, 3) - sf
            vmag = sqrt((vxm**2 + vym**2 + vzm**2)) - sf
            am = sqrt(gammam * pm / rhom)
            mnm = vmag / am
 
            cellarea = sqrt(ssi(i, j, 1)**2 + ssi(i, j, 2)**2 + ssi(i, j, 3)**2)
            area = area + cellarea * blk
            overcellarea = 1 / cellarea
 
            call computeptot(rhom, vxm, vym, vzm, pm, ptot)
            call computettot(rhom, vxm, vym, vzm, pm, ttot)
 
            massflowratelocal = rhom * vnm * blk * fact * mredim
 
            massflowrate = massflowrate + massflowratelocal
 
            ! re-dimentionalize quantities
            pm = pm * pref
 
            mass_ptot = mass_ptot + ptot * massflowratelocal * pref
            mass_ttot = mass_ttot + ttot * massflowratelocal * tref
            mass_rho = mass_rho + rhom * massflowratelocal * rhoref
            mass_a = mass_a + am * massflowratelocal * uref
 
            mass_ps = mass_ps + pm * massflowratelocal
            mass_mn = mass_mn + mnm * massflowratelocal
 
            area_ptot = area_ptot + ptot * pref * cellarea * blk
            area_ps = area_ps + pm * cellarea * blk
 
            sfacecoordref(1) = sf * ssi(i, j, 1) * overcellarea
            sfacecoordref(2) = sf * ssi(i, j, 2) * overcellarea
            sfacecoordref(3) = sf * ssi(i, j, 3) * overcellarea
 
            mass_vx = mass_vx + (vxm * uref - sfacecoordref(1)) * massflowratelocal
            mass_vy = mass_vy + (vym * uref - sfacecoordref(2)) * massflowratelocal
            mass_vz = mass_vz + (vzm * uref - sfacecoordref(3)) * massflowratelocal
 
            govgm1 = gammainf / (gammainf - one)
            gm1ovg = one / govgm1
            viconst = two * govgm1 * rgasdim
            ! the prefs in psinf / ptot cancel out so we can just take the ratio
            ! we need to clip the ratio to stay under one. right next to the wall,
            ! the pTot can go below the static free stream pressure. To prevent
            ! nans from the sqrt, we just clip this. This does not affect the computation
            ! because when pTot is this small, the velocities are also small, and the
            ! mdot is almost zero, so the cells in this area don't contribute much
            ! to the mass weighed sum.
            pratio = min(one, one / ptot)
            vilocal = sqrt(viconst * (one - (pratio)**gm1ovg) * ttot * tref)
            mass_vi = mass_vi + vilocal * massflowratelocal
 
            mass_nx = mass_nx + ssi(i, j, 1) * overcellarea * massflowratelocal
            mass_ny = mass_ny + ssi(i, j, 2) * overcellarea * massflowratelocal
            mass_nz = mass_nz + ssi(i, j, 3) * overcellarea * massflowratelocal
 
            xc = fourth * (xx(i, j, 1) + xx(i + 1, j, 1) &
                           + xx(i, j + 1, 1) + xx(i + 1, j + 1, 1)) - refpoint(1)
            yc = fourth * (xx(i, j, 2) + xx(i + 1, j, 2) &
                           + xx(i, j + 1, 2) + xx(i + 1, j + 1, 2)) - refpoint(2)
            zc = fourth * (xx(i, j, 3) + xx(i + 1, j, 3) &
                           + xx(i, j + 1, 3) + xx(i + 1, j + 1, 3)) - refpoint(3)
 
            ! Pressure forces. Note that these need a *negative* and to subtract
            ! the reference pressure sign to be consistent with the force
            ! computation on the walls.
            pm = -(pm - pinf * pref) * fact * blk
 
            fx = pm * ssi(i, j, 1)
            fy = pm * ssi(i, j, 2)
            fz = pm * ssi(i, j, 3)
 
            ! Update the pressure force and moment coefficients.
            fp(1) = fp(1) + fx
            fp(2) = fp(2) + fy
            fp(3) = fp(3) + fz
 
            mx = yc * fz - zc * fy
            my = zc * fx - xc * fz
            mz = xc * fy - yc * fx
 
            mp(1) = mp(1) + mx
            mp(2) = mp(2) + my
            mp(3) = mp(3) + mz
 
            ! the force integral for the center of pressure computation.
            ! We need the cell centers wrt origin
            xco = fourth * (xx(i, j, 1) + xx(i + 1, j, 1) &
                            + xx(i, j + 1, 1) + xx(i + 1, j + 1, 1))
            yco = fourth * (xx(i, j, 2) + xx(i + 1, j, 2) &
                            + xx(i, j + 1, 2) + xx(i + 1, j + 1, 2))
            zco = fourth * (xx(i, j, 3) + xx(i + 1, j, 3) &
                            + xx(i, j + 1, 3) + xx(i + 1, j + 1, 3))
 
            ! Center of force computations. Here we accumulate in the sums.
            ! accumulate in the sums. each force component is tracked separately
            ! blanking is included in the mdot multiplier for the force.
 
            ! Force-X
            cofsumfx(1) = cofsumfx(1) + xco * fx
            cofsumfx(2) = cofsumfx(2) + yco * fx
            cofsumfx(3) = cofsumfx(3) + zco * fx
 
            ! Force-Y
            cofsumfy(1) = cofsumfy(1) + xco * fy
            cofsumfy(2) = cofsumfy(2) + yco * fy
            cofsumfy(3) = cofsumfy(3) + zco * fy
 
            ! Force-Z
            cofsumfz(1) = cofsumfz(1) + xco * fz
            cofsumfz(2) = cofsumfz(2) + yco * fz
            cofsumfz(3) = cofsumfz(3) + zco * fz
 
            ! Momentum forces are a little tricky.  We negate because
            ! have to re-apply fact to massFlowRateLocal to undoo it, because
            ! we need the signed behavior of ssi to get the momentum forces correct.
            ! Also, the sign is flipped between inflow and outflow types
 
            massflowratelocal = massflowratelocal * fact / timeref * blk / cellarea * internalflowfact * inflowfact
 
            fx = massflowratelocal * ssi(i, j, 1) * vxm
            fy = massflowratelocal * ssi(i, j, 2) * vym
            fz = massflowratelocal * ssi(i, j, 3) * vzm
 
            fmom(1) = fmom(1) + fx
            fmom(2) = fmom(2) + fy
            fmom(3) = fmom(3) + fz
 
            mx = yc * fz - zc * fy
            my = zc * fx - xc * fz
            mz = xc * fy - yc * fx
 
            mmom(1) = mmom(1) + mx
            mmom(2) = mmom(2) + my
            mmom(3) = mmom(3) + mz
 
            ! Center of force computations. Here we accumulate in the sums.
            ! each force component is tracked separately
            ! blanking is included in the mdot multiplier for the force.
 
            ! Force-X
            cofsumfx(1) = cofsumfx(1) + xco * fx
            cofsumfx(2) = cofsumfx(2) + yco * fx
            cofsumfx(3) = cofsumfx(3) + zco * fx
 
            ! Force-Y
            cofsumfy(1) = cofsumfy(1) + xco * fy
            cofsumfy(2) = cofsumfy(2) + yco * fy
            cofsumfy(3) = cofsumfy(3) + zco * fy
 
            ! Force-Z
            cofsumfz(1) = cofsumfz(1) + xco * fz
            cofsumfz(2) = cofsumfz(2) + yco * fz
            cofsumfz(3) = cofsumfz(3) + zco * fz
 
        end do
 
        ! Increment the local values array with what we computed here
        localvalues(imassflow) = localvalues(imassflow) + massflowrate
        localvalues(iarea) = localvalues(iarea) + area
        localvalues(imassrho) = localvalues(imassrho) + mass_rho
        localvalues(imassa) = localvalues(imassa) + mass_a
        localvalues(imassptot) = localvalues(imassptot) + mass_ptot
        localvalues(imassttot) = localvalues(imassttot) + mass_ttot
        localvalues(imassps) = localvalues(imassps) + mass_ps
        localvalues(imassmn) = localvalues(imassmn) + mass_mn
        localvalues(ifp:ifp + 2) = localvalues(ifp:ifp + 2) + fp
        localvalues(iflowfm:iflowfm + 2) = localvalues(iflowfm:iflowfm + 2) + fmom
        localvalues(iflowmp:iflowmp + 2) = localvalues(iflowmp:iflowmp + 2) + mp
        localvalues(iflowmm:iflowmm + 2) = localvalues(iflowmm:iflowmm + 2) + mmom
 
        localvalues(icoforcex:icoforcex + 2) = localvalues(icoforcex:icoforcex + 2) + cofsumfx
        localvalues(icoforcey:icoforcey + 2) = localvalues(icoforcey:icoforcey + 2) + cofsumfy
        localvalues(icoforcez:icoforcez + 2) = localvalues(icoforcez:icoforcez + 2) + cofsumfz
 
        localvalues(iareaptot) = localvalues(iareaptot) + area_ptot
        localvalues(iareaps) = localvalues(iareaps) + area_ps
 
        localvalues(imassvx) = localvalues(imassvx) + mass_vx
        localvalues(imassvy) = localvalues(imassvy) + mass_vy
        localvalues(imassvz) = localvalues(imassvz) + mass_vz
        localvalues(imassnx) = localvalues(imassnx) + mass_nx
        localvalues(imassny) = localvalues(imassny) + mass_ny
        localvalues(imassnz) = localvalues(imassnz) + mass_nz
 
        localvalues(imassvi) = localvalues(imassvi) + mass_vi
 
    end subroutine flowintegrationface
 
    ! ----------------------------------------------------------------------
    !                                                                      |
    !                    No Tapenade Routine below this line               |
    !                                                                      |
    ! ----------------------------------------------------------------------
 
#ifndef USE_TAPENADE
 
    subroutine getsolutionwrap(famLists, funcValues, nCost, nGroups, nFamMax)
 
        use constants
        use inputtimespectral, only: ntimeintervalsspectral
        implicit none
        ! Input/output Variables
        integer(kind=intType) :: nGroups, nCost, nFamMax
        integer(kind=intType), dimension(nGroups, nFamMax) :: famLists
        real(kind=realtype), dimension(nCost, nGroups), intent(out) :: funcvalues
 
        ! Local variable
 
        call getsolution(famlists, funcvalues)
    end subroutine getsolutionwrap
 
    subroutine getsolution(famLists, funcValues, globalValues)
        !--------------------------------------------------------------
        ! Manual Differentiation Warning: Modifying this routine requires
        ! modifying the hand-written forward and reverse routines.
        ! --------------------------------------------------------------
 
        use constants
        use inputtimespectral, only: ntimeintervalsspectral
        use communication, only: adflow_comm_world
        use blockpointers, only: ndom
        use utils, only: setpointers, echk
        use zipperintegrations, only: integratezippers
        use usersurfaceintegrations, only: integrateusersurfaces
        use actuatorregion, only: integrateactuatorregions
        use inputcostfunctions, only: computecavitation, computesepsensorks
        implicit none
 
        ! Input/Output Variables
        integer(kind=intType), dimension(:, :), intent(in), target :: famLists
        real(kind=realtype), dimension(:, :), intent(out) :: funcvalues
        real(kind=realtype), optional, intent(out), dimension(:, :, :) :: globalvalues
 
        ! Working
        real(kind=realtype), dimension(nLocalValues, nTimeIntervalsSpectral) :: localval, globalval
        integer(kind=intType) :: nn, sps, ierr, iGroup, nFam
        integer(kind=intType), dimension(:), pointer :: famList
        ! Master loop over the each of the groups we have
 
        grouploop: do igroup = 1, size(famlists, 1)
 
            ! Extract the current family list
            nfam = famlists(igroup, 1)
            famlist => famlists(igroup, 2:2 + nfam - 1)
            funcvalues(:, igroup) = zero
            localval = zero
 
            ! compute the current cp min value for the cavitation computation with KS aggregation
            if (computecavitation) then
                call computecpminfamily(famlist)
            end if
 
            ! compute the current sepsensor max value for the separation computation with KS aggregation
            if (computesepsensorks) then
                call computesepsenmaxfamily(famlist)
            end if
 
            do sps = 1, ntimeintervalsspectral
                ! Integrate the normal block surfaces.
                do nn = 1, ndom
                    call setpointers(nn, 1, sps)
                    call integratesurfaces(localval(:, sps), famlist)
                end do
 
                ! Integrate any zippers we have
                call integratezippers(localval(:, sps), famlist, sps)
 
                ! Integrate any user-supplied surfaces as have as well.
                call integrateusersurfaces(localval(:, sps), famlist, sps)
 
                ! Integrate any actuator regions we have
                call integrateactuatorregions(localval(:, sps), famlist, sps)
            end do
 
            ! Now we need to reduce all the cost functions
            call mpi_allreduce(localval, globalval, nlocalvalues * ntimeintervalsspectral, adflow_real, &
                               mpi_sum, adflow_comm_world, ierr)
            call echk(ierr, __file__, __line__)
 
            ! Call the final routine that will comptue all of our functions of
            ! interest.
            call getcostfunctions(globalval, funcvalues(:, igroup))
 
            if (present(globalvalues)) then
                globalvalues(:, :, igroup) = globalval
            end if
        end do grouploop
    end subroutine getsolution
 
    subroutine integratesurfaces(localValues, famList)
        !--------------------------------------------------------------
        ! Manual Differentiation Warning: Modifying this routine requires
        ! modifying the hand-written forward and reverse routines.
        ! --------------------------------------------------------------
        !
        ! This is a shell routine that calls the specific surface
        ! integration routines. Currently we have have the forceAndMoment
        ! routine as well as the flow properties routine. This routine
        ! takes care of setting pointers, while the actual computational
        ! routine just acts on a specific fast pointed to by pointers.
 
        use constants
        use blockpointers, only: nbocos, bcdata, bctype, sk, sj, si, x, rlv, &
                                 sfacei, sfacej, sfacek, gamma, rev, p, viscsubface
        use utils, only: setbcpointers, iswalltype
        use sorting, only: faminlist
        ! Tapenade needs to see these modules that the callees use.
        use bcpointers
        use flowvarrefstate
        use inputphysics
 
        implicit none
 
        ! Input/output Variables
        real(kind=realtype), dimension(nLocalValues), intent(inout) :: localvalues
        integer(kind=intType), dimension(:), intent(in) :: famList
 
        ! Working variables
        integer(kind=intType) :: mm
 
        ! Loop over all possible boundary conditions
        bocos: do mm = 1, nbocos
 
            ! Determine if this boundary condition is to be incldued in the
            ! currently active group
            faminclude: if (faminlist(bcdata(mm)%famID, famlist)) then
 
                ! Set a bunch of pointers depending on the face id to make
                ! a generic treatment possible.
                call setbcpointers(mm, .true.)
 
                ! no post gathered integrations currently
                iswall: if (iswalltype(bctype(mm))) then
                    call wallintegrationface(localvalues, mm)
                end if iswall
 
                isinflowoutflow: if (bctype(mm) == subsonicinflow .or. &
                                     bctype(mm) == supersonicinflow) then
                    call flowintegrationface(.true., localvalues, mm)
                else if (bctype(mm) == subsonicoutflow .or. &
                         bctype(mm) == supersonicoutflow) then
                    call flowintegrationface(.false., localvalues, mm)
                end if isinflowoutflow
 
            end if faminclude
        end do bocos
 
    end subroutine integratesurfaces
 
    subroutine computecpminfamily(famList)
 
        use constants
        use inputtimespectral, only: ntimeintervalsspectral
        use communication, only: adflow_comm_world, myid
        use blockpointers, only: ndom
        use inputphysics, only: cpmin_family, machcoef
        use blockpointers
        use flowvarrefstate
        use bcpointers
        use utils, only: setpointers, setbcpointers, iswalltype, echk
        use sorting, only: faminlist
 
        implicit none
 
        integer(kind=intType), dimension(:), intent(in) :: famList
        integer(kind=intType) :: mm, nn, sps
        integer(kind=intType) :: i, j, ii, blk, ierr
        real(kind=realtype) :: cp, plocal, tmp, cpmin_local
 
        ! this routine loops over the surface cells in the given family
        ! and computes the true minimum Cp value.
        ! this is then used in the surface integration routine to compute
        ! the cpmin using KS aggregation.
        ! the goal is to get a differentiable cpmin output.
 
        ! loop over the TS instances and compute cpmin_family for each TS instance
        do sps = 1, ntimeintervalsspectral
            ! set the local cp min to a large value so that we get the actual min
            cpmin_local = 10000.0_realtype
            do nn = 1, ndom
                call setpointers(nn, 1, sps)
 
                ! Loop over all possible boundary conditions
                bocos: do mm = 1, nbocos
                    ! Determine if this boundary condition is to be incldued in the
                    ! currently active group
                    faminclude: if (faminlist(bcdata(mm)%famID, famlist)) then
 
                        ! Set a bunch of pointers depending on the face id to make
                        ! a generic treatment possible.
                        call setbcpointers(mm, .true.)
 
                        ! no post gathered integrations currently
                        iswall: if (iswalltype(bctype(mm))) then
 
                            !$AD II-LOOP
                            do ii = 0, (bcdata(mm)%jnEnd - bcdata(mm)%jnBeg) * (bcdata(mm)%inEnd - bcdata(mm)%inBeg) - 1
                                i = mod(ii, (bcdata(mm)%inEnd - bcdata(mm)%inBeg)) + bcdata(mm)%inBeg + 1
                                j = ii / (bcdata(mm)%inEnd - bcdata(mm)%inBeg) + bcdata(mm)%jnBeg + 1
 
                                ! only take this if its a compute cell
                                if (bcdata(mm)%iblank(i, j) .eq. one) then
 
                                    ! compute local CP
                                    plocal = pp2(i, j)
                                    tmp = two / (gammainf * machcoef * machcoef)
                                    cp = tmp * (plocal - pinf)
 
                                    ! compare it against the current value on this proc
                                    cpmin_local = min(cpmin_local, cp)
                                end if
                            end do
                        end if iswall
 
                    end if faminclude
                end do bocos
            end do
            ! finally communicate across all processors for this time spectral instance
            call mpi_allreduce(cpmin_local, cpmin_family(sps), 1, mpi_double, &
                               mpi_min, adflow_comm_world, ierr)
            call echk(ierr, __file__, __line__)
        end do
 
    end subroutine computecpminfamily
 
    subroutine computesepsenmaxfamily(famList)
 
        use constants
        use inputtimespectral, only: ntimeintervalsspectral
        use communication, only: adflow_comm_world, myid
        use blockpointers, only: ndom
        use inputphysics, only: sepsenmaxfamily, veldirfreestream
        use blockpointers
        use flowvarrefstate
        use inputcostfunctions, only: sepsensorksphi
        use bcpointers
        use utils, only: setpointers, setbcpointers, iswalltype, echk
        use sorting, only: faminlist
 
        implicit none
 
        integer(kind=intType), dimension(:), intent(in) :: famList
        integer(kind=intType) :: mm, nn, sps
        integer(kind=intType) :: i, j, ii, blk, ierr
        real(kind=realtype) :: sepsensor_local
        real(kind=realtype) :: vecttangential(3), v(3)
        real(kind=realtype) :: vectdotproductfsnormal, sensor
 
        ! this routine loops over the surface cells in the given family
        ! and computes the true maximum sepsensor value.
        ! this is then used in the surface integration routine to compute
        ! the maximum sepsensor value using KS aggregation.
        ! the goal is to get a differentiable maximum sepsensor output.
 
        ! loop over the TS instances and compute sepSenMaxFamily for each TS instance
        do sps = 1, ntimeintervalsspectral
            ! set the local sepsensor to a smaller value so that we get the actual max
            sepsensor_local = -10000.0_realtype
            do nn = 1, ndom
                call setpointers(nn, 1, sps)
 
                ! Loop over all possible boundary conditions
                bocos: do mm = 1, nbocos
                    ! Determine if this boundary condition is to be incldued in the
                    ! currently active group
                    faminclude: if (faminlist(bcdata(mm)%famID, famlist)) then
 
                        ! Set a bunch of pointers depending on the face id to make
                        ! a generic treatment possible.
                        call setbcpointers(mm, .true.)
 
                        ! no post gathered integrations currently
                        iswall: if (iswalltype(bctype(mm))) then
 
                            !$AD II-LOOP
                            do ii = 0, (bcdata(mm)%jnEnd - bcdata(mm)%jnBeg) * (bcdata(mm)%inEnd - bcdata(mm)%inBeg) - 1
                                i = mod(ii, (bcdata(mm)%inEnd - bcdata(mm)%inBeg)) + bcdata(mm)%inBeg + 1
                                j = ii / (bcdata(mm)%inEnd - bcdata(mm)%inBeg) + bcdata(mm)%jnBeg + 1
 
                                ! only take this if its a compute cell
                                if (bcdata(mm)%iblank(i, j) .eq. one) then
                                    ! Get normalized surface velocity:
                                    v(1) = ww2(i, j, ivx)
                                    v(2) = ww2(i, j, ivy)
                                    v(3) = ww2(i, j, ivz)
                                    v = v / (sqrt(v(1)**2 + v(2)**2 + v(3)**2) + 1e-16)
 
                                    vectdotproductfsnormal = veldirfreestream(1) * bcdata(mm)%norm(i, j, 1) + &
                                                             veldirfreestream(2) * bcdata(mm)%norm(i, j, 2) + &
                                                             veldirfreestream(3) * bcdata(mm)%norm(i, j, 3)
                                    ! Tangential Vector on the surface, which is the freestream projected vector
                                    vecttangential(1) = veldirfreestream(1) - &
                                                        vectdotproductfsnormal * bcdata(mm)%norm(i, j, 1)
                                    vecttangential(2) = veldirfreestream(2) - &
                                                        vectdotproductfsnormal * bcdata(mm)%norm(i, j, 2)
                                    vecttangential(3) = veldirfreestream(3) - &
                                                        vectdotproductfsnormal * bcdata(mm)%norm(i, j, 3)
 
                                    vecttangential = vecttangential / &
                                                     (sqrt(vecttangential(1)**2 + vecttangential(2)**2 + &
                                                           vecttangential(3)**2) + 1e-16)
 
                                    ! computing separation sensor
                                    ! velocity dot products
                                    sensor = (v(1) * vecttangential(1) + v(2) * vecttangential(2) + &
                                              v(3) * vecttangential(3))
 
                                    ! sepsensor value
                                    sensor = (cos(degtorad * sepsensorksphi) - sensor) / &
                                             (-cos(degtorad * sepsensorksphi) + cos(zero) + 1e-16)
 
                                    ! compare it against the current value on this proc
                                    sepsensor_local = max(sepsensor_local, sensor)
                                end if
 
                                ! end if
                            end do
                        end if iswall
 
                    end if faminclude
                end do bocos
            end do
            ! finally communicate across all processors for this time spectral instance
            call mpi_allreduce(sepsensor_local, sepsenmaxfamily(sps), 1, mpi_double, &
                               mpi_max, adflow_comm_world, ierr)
            call echk(ierr, __file__, __line__)
        end do
 
    end subroutine computesepsenmaxfamily
 
#ifndef USE_COMPLEX
    subroutine integratesurfaces_d(localValues, localValuesd, famList)
        !------------------------------------------------------------------------
        ! Manual Differentiation Warning: This routine is differentiated by hand.
        ! -----------------------------------------------------------------------
 
        ! Forward mode linearization of integrateSurfaces
        use constants
        use blockpointers, only: nbocos, bcdata, bctype
        use utils, only: setbcpointers_d, iswalltype
        use sorting, only: faminlist
        use surfaceintegrations_d, only: wallintegrationface_d, flowintegrationface_d
        implicit none
 
        ! Input/output Variables
        real(kind=realtype), dimension(nLocalValues), intent(inout) :: localvalues, localvaluesd
        integer(kind=intType), dimension(:), intent(in) :: famList
 
        ! Working variables
        integer(kind=intType) :: mm
 
        ! Loop over all possible boundary conditions
        do mm = 1, nbocos
            ! Determine if this boundary condition is to be incldued in the
            ! currently active group
            faminclude: if (faminlist(bcdata(mm)%famID, famlist)) then
 
                ! Set a bunch of pointers depending on the face id to make
                ! a generic treatment possible.
                call setbcpointers_d(mm, .true.)
 
                ! not post gathered integrations currently
                iswall: if (iswalltype(bctype(mm))) then
                    call wallintegrationface_d(localvalues, localvaluesd, mm)
                end if iswall
 
                isinflowoutflow: if (bctype(mm) == subsonicinflow .or. &
                                     bctype(mm) == supersonicinflow) then
                    call flowintegrationface_d(.true., localvalues, localvaluesd, mm)
                else if (bctype(mm) == subsonicoutflow .or. &
                         bctype(mm) == supersonicoutflow) then
 
                    call flowintegrationface_d(.false., localvalues, localvaluesd, mm)
 
                end if isinflowoutflow
            end if faminclude
        end do
    end subroutine integratesurfaces_d
 
    subroutine integratesurfaces_b(localValues, localValuesd, famList)
        !------------------------------------------------------------------------
        ! Manual Differentiation Warning: This routine is differentiated by hand.
        ! -----------------------------------------------------------------------
 
        ! Reverse mode linearization of integrateSurfaces
        use constants
        use blockpointers, only: nbocos, bcdata, bctype, bcdatad
        use utils, only: setbcpointers_d, iswalltype
        use sorting, only: faminlist
        use surfaceintegrations_b, only: wallintegrationface_b, flowintegrationface_b
        implicit none
 
        ! Input/output Variables
        real(kind=realtype), dimension(nLocalValues), intent(inout) :: localvalues, localvaluesd
        integer(kind=intType), dimension(:), intent(in) :: famList
        ! Working variables
        integer(kind=intType) :: mm
 
        ! Call the individual integration routines.
        do mm = 1, nbocos
            ! Determine if this boundary condition is to be incldued in the
            ! currently active group
            faminclude: if (faminlist(bcdata(mm)%famID, famlist)) then
 
                ! Set a bunch of pointers depending on the face id to make
                ! a generic treatment possible.
                call setbcpointers_d(mm, .true.)
 
                ! not post gathered integrations currently
                iswall: if (iswalltype(bctype(mm))) then
                    call wallintegrationface_b(localvalues, localvaluesd, mm)
                end if iswall
 
                isinflowoutflow: if (bctype(mm) == subsonicinflow .or. &
                                     bctype(mm) == supersonicinflow) then
                    call flowintegrationface_b(.true., localvalues, localvaluesd, mm)
                else if (bctype(mm) == subsonicoutflow .or. &
                         bctype(mm) == supersonicoutflow) then
                    call flowintegrationface_b(.false., localvalues, localvaluesd, mm)
                end if isinflowoutflow
            end if faminclude
        end do
    end subroutine integratesurfaces_b
 
    subroutine getsolution_d(famLists, funcValues, funcValuesd)
        !------------------------------------------------------------------------
        ! Manual Differentiation Warning: This routine is differentiated by hand.
        ! -----------------------------------------------------------------------
 
        use constants
        use inputtsstabderiv, only: tsstability
        use inputtimespectral, only: ntimeintervalsspectral
        use communication, only: adflow_comm_world
        use blockpointers, only: ndom
        use utils, only: setpointers_d, echk
        use surfaceintegrations_d, only: getcostfunctions_d
        use zipperintegrations, only: integratezippers_d
        use usersurfaceintegrations, only: integrateusersurfaces_d
        use actuatorregion, only: integrateactuatorregions_d
        use inputcostfunctions, only: computecavitation, computesepsensorks
        implicit none
 
        ! Input/Output Variables
        integer(kind=intType), dimension(:, :), target, intent(in) :: famLists
        real(kind=realtype), dimension(:, :), intent(out) :: funcvalues, funcvaluesd
 
        ! Working
        real(kind=realtype), dimension(nLocalValues, nTimeIntervalsSpectral) :: localval, globalval
        real(kind=realtype), dimension(nLocalValues, nTimeIntervalsSpectral) :: localvald, globalvald
        integer(kind=intType) :: nn, sps, ierr, iGroup, nFam
        integer(kind=intType), dimension(:), pointer :: famList
 
        grouploop: do igroup = 1, size(famlists, 1)
 
            ! Extract the current family list
            nfam = famlists(igroup, 1)
            famlist => famlists(igroup, 2:2 + nfam - 1)
 
            ! compute the current cp min value for the cavitation computation with KS aggregation
            if (computecavitation) then
                call computecpminfamily(famlist)
            end if
 
            ! compute the current sepsensor max value for the separation computation with KS aggregation
            if (computesepsensorks) then
                call computesepsenmaxfamily(famlist)
            end if
 
            localval = zero
            localvald = zero
            do sps = 1, ntimeintervalsspectral
                ! Integrate the normal block surfaces.
                do nn = 1, ndom
                    call setpointers_d(nn, 1, sps)
                    call integratesurfaces_d(localval(:, sps), localvald(:, sps), famlist)
                end do
 
                ! Integrate any zippers we have
                call integratezippers_d(localval(:, sps), localvald(:, sps), famlist, sps)
 
                ! Integrate any user-supplied surface as have as well.
                call integrateusersurfaces_d(localval(:, sps), localvald(:, sps), famlist, sps)
 
                ! Integrate any actuator regions we have
                call integrateactuatorregions_d(localval(:, sps), localvald(:, sps), famlist, sps)
            end do
 
            ! Now we need to reduce all the cost functions
            call mpi_allreduce(localval, globalval, nlocalvalues * ntimeintervalsspectral, adflow_real, &
                               mpi_sum, adflow_comm_world, ierr)
            call echk(ierr, __file__, __line__)
 
            ! Now we need to reduce all the cost functions
            call mpi_allreduce(localvald, globalvald, nlocalvalues * ntimeintervalsspectral, adflow_real, &
                               mpi_sum, adflow_comm_world, ierr)
            call echk(ierr, __file__, __line__)
 
            ! Call the final routine that will comptue all of our functions of
            ! interest.
 
            call getcostfunctions_d(globalval, globalvald, funcvalues(:, igroup), funcvaluesd(:, igroup))
 
            ! if (present(globalValues)) then
            !    globalValues = globalVal
            !    globalValuesd = globalVald
            ! end if
        end do grouploop
 
    end subroutine getsolution_d
 
    subroutine getsolution_b(famLists, funcValues, funcValuesd)
        ! -----------------------------------------------------------------------
        ! Manual Differentiation Warning: This routine is differentiated by hand.
        ! -----------------------------------------------------------------------
 
        use constants
        use communication, only: myid
        use inputtsstabderiv, only: tsstability
        use inputtimespectral, only: ntimeintervalsspectral
        use communication, only: adflow_comm_world
        use blockpointers, only: ndom
        use utils, only: setpointers_b, echk, setpointers
        use surfaceintegrations_b, only: getcostfunctions_b
        use zipperintegrations, only: integratezippers_b
        use usersurfaceintegrations, only: integrateusersurfaces_b
        use actuatorregion, only: integrateactuatorregions_b
        use inputcostfunctions, only: computecavitation, computesepsensorks
        implicit none
 
        ! Input/Output Variables
        integer(kind=intType), dimension(:, :), target, intent(in) :: famLists
        real(kind=realtype), dimension(:, :) :: funcvalues, funcvaluesd
 
        ! Working
        real(kind=realtype), dimension(nLocalValues, nTimeIntervalsSpectral) :: localval, globalval
        real(kind=realtype), dimension(nLocalValues, nTimeIntervalsSpectral) :: localvald, globalvald
        real(kind=realtype), dimension(nLocalValues, nTimeIntervalsSpectral, size(famLists, 1)) :: globalvalues
        integer(kind=intType) :: nn, sps, ierr, iGroup, nFam
        integer(kind=intType), dimension(:), pointer :: famList
 
        call getsolution(famlists, funcvalues, globalvalues)
 
        grouploop: do igroup = 1, size(famlists, 1)
 
            ! Extract the current family list
            nfam = famlists(igroup, 1)
            famlist => famlists(igroup, 2:2 + nfam - 1)
 
            ! compute the current cp min value for the cavitation computation with KS aggregation
            if (computecavitation) then
                call computecpminfamily(famlist)
            end if
 
            ! compute the current sepsensor max value for the separation computation with KS aggregation
            if (computesepsensorks) then
                call computesepsenmaxfamily(famlist)
            end if
 
            localval = zero
            localvald = zero
 
            ! Retrive the forward pass values from getSolution
            globalval = globalvalues(:, :, igroup)
 
            if (myid == 0) then
                call getcostfunctions_b(globalval, globalvald, funcvalues(:, igroup), funcvaluesd(:, igroup))
                localvald = globalvald
            end if
 
            ! Now we need to bcast out the localValues to all procs.
            call mpi_bcast(localvald, nlocalvalues * ntimeintervalsspectral, &
                           adflow_real, 0, adflow_comm_world, ierr)
            call echk(ierr, __file__, __line__)
 
            do sps = 1, ntimeintervalsspectral
 
                ! Integrate any actuator regions we have:
                call integrateactuatorregions_b(localval(:, sps), localvald(:, sps), famlist, sps)
 
                ! Integrate any user-supplied planes as have as well.
                call integrateusersurfaces_b(localval(:, sps), localvald(:, sps), famlist, sps)
 
                ! Integrate any zippers we have
                call integratezippers_b(localval(:, sps), localvald(:, sps), famlist, sps)
 
                ! Integrate the normal block surfaces.
                do nn = 1, ndom
                    call setpointers_b(nn, 1, sps)
                    call integratesurfaces_b(localval(:, sps), localvald(:, sps), famlist)
                end do
 
            end do
        end do grouploop
    end subroutine getsolution_b
#endif
#endif
end module surfaceintegrations