surfaceIntegrations_b.f90

Go to the documentation of this file.

!        generated by tapenade     (inria, ecuador team)
!  tapenade 3.16 (develop) - 22 aug 2023 15:51
!
module surfaceintegrations_b
  implicit none
 
contains
!  differentiation of getcostfunctions in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: pref machcoef dragdirection
!                liftdirection funcvalues
!   with respect to varying inputs: pref machcoef dragdirection
!                liftdirection globalvals funcvalues
!   rw status of diff variables: pref:incr machcoef:incr dragdirection:incr
!                liftdirection:incr globalvals:out funcvalues:in-zero
  subroutine getcostfunctions_b(globalvals, globalvalsd, funcvalues, &
&   funcvaluesd)
    use constants
    use inputtimespectral, only : ntimeintervalsspectral
    use flowvarrefstate, only : pref, prefd, rhoref, rhorefd, tref, &
&   trefd, lref, gammainf, pinf, pinfd, uref, urefd, uinf, uinfd
    use inputphysics, only : liftdirection, liftdirectiond, &
&   dragdirection, dragdirectiond, surfaceref, machcoef, machcoefd, &
&   lengthref, alpha, alphad, beta, betad, liftindex, cpmin_family, &
&   cpmin_rho, sepsenmaxfamily, sepsenmaxrho
    use inputcostfunctions, only : computecavitation, computesepsensorks&
&   , sepsensorkssharpness, sepsensorksoffset
    use inputtsstabderiv, only : tsstability
    use utils_b, only : computetsderivatives
    use flowutils_b, only : getdirvector
    implicit none
! input/output
    real(kind=realtype), dimension(:, :), intent(in) :: globalvals
    real(kind=realtype), dimension(:, :) :: globalvalsd
    real(kind=realtype), dimension(:) :: funcvalues
    real(kind=realtype), dimension(:) :: funcvaluesd
! working
    real(kind=realtype) :: fact, factmoment, ovrnts
    real(kind=realtype) :: factd
    real(kind=realtype), dimension(3, ntimeintervalsspectral) :: force, &
&   forcep, forcev, forcem, moment, cforce, cforcep, cforcev, cforcem, &
&   cmoment, cofx, cofy, cofz
    real(kind=realtype), dimension(3, ntimeintervalsspectral) :: forced&
&   , forcepd, forcevd, forcemd, momentd, cforced, cforcepd, cforcevd, &
&   cforcemd, cmomentd, cofxd, cofyd, cofzd
    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) :: mavgptotd, mavgttotd, mavgrhod, mavgpsd, &
&   mflowd, mavgmnd, mavgad, mavgvxd, mavgvyd, mavgvzd, garead, mavgvid&
&   , fxliftd, fyliftd, fzliftd
    real(kind=realtype) :: vdotn, mag, u, v, w, ks_comp
    real(kind=realtype) :: ks_compd
    integer(kind=inttype) :: sps
    real(kind=realtype), dimension(8) :: dcdq, dcdqdot
    real(kind=realtype), dimension(8) :: dcdalpha, dcdalphadot
    real(kind=realtype), dimension(8) :: coef0
    intrinsic log
    intrinsic exp
    intrinsic sqrt
    real(kind=realtype) :: temp
    real(kind=realtype) :: temp0
    real(kind=realtype) :: temp1
    real(kind=realtype) :: tempd
    real(kind=realtype) :: tmp
    real(kind=realtype) :: tmpd
    real(kind=realtype) :: tmp0
    real(kind=realtype) :: tmpd0
    real(kind=realtype) :: tmp1
    real(kind=realtype) :: tmpd1
    real(kind=realtype) :: tmp2
    real(kind=realtype) :: tmpd2
    real(kind=realtype) :: tmp3
    real(kind=realtype) :: tmpd3
    real(kind=realtype) :: tmp4
    real(kind=realtype) :: tmpd4
    real(kind=realtype) :: tmp5
    real(kind=realtype) :: tmpd5
    real(kind=realtype) :: tmp6
    real(kind=realtype) :: tmpd6
    real(kind=realtype) :: tmp7
    real(kind=realtype) :: tmpd7
    real(kind=realtype) :: tmp8
    real(kind=realtype) :: tmpd8
    real(kind=realtype) :: tmp9
    real(kind=realtype) :: tmpd9
    real(kind=realtype) :: tmp10
    real(kind=realtype) :: tmpd10
    real(kind=realtype) :: tmp11
    real(kind=realtype) :: tmpd11
    real(kind=realtype) :: tmp12
    real(kind=realtype) :: tmpd12
    real(kind=realtype) :: tmp13
    real(kind=realtype) :: tmpd13
    real(kind=realtype) :: tmp14
    real(kind=realtype) :: tmpd14
    real(kind=realtype) :: tmp15
    real(kind=realtype) :: tmpd15
    real(kind=realtype) :: tempd0
    real(kind=realtype) :: tmp16
    real(kind=realtype) :: tmpd16
    real(kind=realtype) :: tmp17
    real(kind=realtype) :: tmpd17
    real(kind=realtype) :: tempd1
    integer :: branch
! 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
    call pushreal8(fact)
    fact = fact/(lengthref*lref)
    cmoment = fact*moment
! zero values since we are summing.
    funcvalues = zero
    call pushreal8array(cofx, 3*ntimeintervalsspectral)
    call pushreal8array(cofy, 3*ntimeintervalsspectral)
    call pushreal8array(cofz, 3*ntimeintervalsspectral)
!$fwd-of ii-loop 
! here we finally assign the final function values
    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) .ne. zero) then
        cofx(:, sps) = cofx(:, sps)/force(1, sps)
      else
        cofx(:, sps) = zero
      end if
      if (force(2, sps) .ne. zero) then
        cofy(:, sps) = cofy(:, sps)/force(2, sps)
      else
        cofy(:, sps) = zero
      end if
      if (force(3, sps) .ne. 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 .ne. 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 .ne. 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.
    tmp = funcvalues(costfuncforcex)*liftdirection(1) + funcvalues(&
&     costfuncforcey)*liftdirection(2) + funcvalues(costfuncforcez)*&
&     liftdirection(3)
    call pushreal8(funcvalues(costfunclift))
    funcvalues(costfunclift) = tmp
    tmp0 = funcvalues(costfuncforcexpressure)*liftdirection(1) + &
&     funcvalues(costfuncforceypressure)*liftdirection(2) + funcvalues(&
&     costfuncforcezpressure)*liftdirection(3)
    call pushreal8(funcvalues(costfuncliftpressure))
    funcvalues(costfuncliftpressure) = tmp0
    tmp1 = funcvalues(costfuncforcexviscous)*liftdirection(1) + &
&     funcvalues(costfuncforceyviscous)*liftdirection(2) + funcvalues(&
&     costfuncforcezviscous)*liftdirection(3)
    call pushreal8(funcvalues(costfuncliftviscous))
    funcvalues(costfuncliftviscous) = tmp1
    tmp2 = funcvalues(costfuncforcexmomentum)*liftdirection(1) + &
&     funcvalues(costfuncforceymomentum)*liftdirection(2) + funcvalues(&
&     costfuncforcezmomentum)*liftdirection(3)
    call pushreal8(funcvalues(costfuncliftmomentum))
    funcvalues(costfuncliftmomentum) = tmp2
!-----
    tmp3 = funcvalues(costfuncforcex)*dragdirection(1) + funcvalues(&
&     costfuncforcey)*dragdirection(2) + funcvalues(costfuncforcez)*&
&     dragdirection(3)
    call pushreal8(funcvalues(costfuncdrag))
    funcvalues(costfuncdrag) = tmp3
    tmp4 = funcvalues(costfuncforcexpressure)*dragdirection(1) + &
&     funcvalues(costfuncforceypressure)*dragdirection(2) + funcvalues(&
&     costfuncforcezpressure)*dragdirection(3)
    call pushreal8(funcvalues(costfuncdragpressure))
    funcvalues(costfuncdragpressure) = tmp4
    tmp5 = funcvalues(costfuncforcexviscous)*dragdirection(1) + &
&     funcvalues(costfuncforceyviscous)*dragdirection(2) + funcvalues(&
&     costfuncforcezviscous)*dragdirection(3)
    call pushreal8(funcvalues(costfuncdragviscous))
    funcvalues(costfuncdragviscous) = tmp5
    tmp6 = funcvalues(costfuncforcexmomentum)*dragdirection(1) + &
&     funcvalues(costfuncforceymomentum)*dragdirection(2) + funcvalues(&
&     costfuncforcezmomentum)*dragdirection(3)
    call pushreal8(funcvalues(costfuncdragmomentum))
    funcvalues(costfuncdragmomentum) = tmp6
!-----
    tmp7 = funcvalues(costfuncforcexcoef)*liftdirection(1) + funcvalues(&
&     costfuncforceycoef)*liftdirection(2) + funcvalues(&
&     costfuncforcezcoef)*liftdirection(3)
    call pushreal8(funcvalues(costfuncliftcoef))
    funcvalues(costfuncliftcoef) = tmp7
    tmp8 = funcvalues(costfuncforcexcoefpressure)*liftdirection(1) + &
&     funcvalues(costfuncforceycoefpressure)*liftdirection(2) + &
&     funcvalues(costfuncforcezcoefpressure)*liftdirection(3)
    call pushreal8(funcvalues(costfuncliftcoefpressure))
    funcvalues(costfuncliftcoefpressure) = tmp8
    tmp9 = funcvalues(costfuncforcexcoefviscous)*liftdirection(1) + &
&     funcvalues(costfuncforceycoefviscous)*liftdirection(2) + &
&     funcvalues(costfuncforcezcoefviscous)*liftdirection(3)
    call pushreal8(funcvalues(costfuncliftcoefviscous))
    funcvalues(costfuncliftcoefviscous) = tmp9
    tmp10 = funcvalues(costfuncforcexcoefmomentum)*liftdirection(1) + &
&     funcvalues(costfuncforceycoefmomentum)*liftdirection(2) + &
&     funcvalues(costfuncforcezcoefmomentum)*liftdirection(3)
    call pushreal8(funcvalues(costfuncliftcoefmomentum))
    funcvalues(costfuncliftcoefmomentum) = tmp10
!-----
    tmp11 = funcvalues(costfuncforcexcoef)*dragdirection(1) + funcvalues&
&     (costfuncforceycoef)*dragdirection(2) + funcvalues(&
&     costfuncforcezcoef)*dragdirection(3)
    call pushreal8(funcvalues(costfuncdragcoef))
    funcvalues(costfuncdragcoef) = tmp11
    tmp12 = funcvalues(costfuncforcexcoefpressure)*dragdirection(1) + &
&     funcvalues(costfuncforceycoefpressure)*dragdirection(2) + &
&     funcvalues(costfuncforcezcoefpressure)*dragdirection(3)
    call pushreal8(funcvalues(costfuncdragcoefpressure))
    funcvalues(costfuncdragcoefpressure) = tmp12
    tmp13 = funcvalues(costfuncforcexcoefviscous)*dragdirection(1) + &
&     funcvalues(costfuncforceycoefviscous)*dragdirection(2) + &
&     funcvalues(costfuncforcezcoefviscous)*dragdirection(3)
    call pushreal8(funcvalues(costfuncdragcoefviscous))
    funcvalues(costfuncdragcoefviscous) = tmp13
    tmp14 = funcvalues(costfuncforcexcoefmomentum)*dragdirection(1) + &
&     funcvalues(costfuncforceycoefmomentum)*dragdirection(2) + &
&     funcvalues(costfuncforcezcoefmomentum)*dragdirection(3)
    call pushreal8(funcvalues(costfuncdragcoefmomentum))
    funcvalues(costfuncdragcoefmomentum) = tmp14
! ----- 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 .ne. zero) then
      tmp15 = (fxlift*funcvalues(costfunccoforcexx)+fylift*funcvalues(&
&       costfunccoforceyx)+fzlift*funcvalues(costfunccoforcezx))/(fxlift&
&       +fylift+fzlift)
      call pushreal8(funcvalues(costfunccofliftx))
      funcvalues(costfunccofliftx) = tmp15
      tmp16 = (fxlift*funcvalues(costfunccoforcexy)+fylift*funcvalues(&
&       costfunccoforceyy)+fzlift*funcvalues(costfunccoforcezy))/(fxlift&
&       +fylift+fzlift)
      call pushreal8(funcvalues(costfunccoflifty))
      funcvalues(costfunccoflifty) = tmp16
      call pushcontrol1b(0)
    else
      call pushcontrol1b(1)
    end if
! -------------------- time spectral objectives ------------------
    if (tsstability) then
      stop
    else
      call popcontrol1b(branch)
      if (branch .eq. 0) then
        tmpd17 = funcvaluesd(costfunccofliftz)
        funcvaluesd(costfunccofliftz) = 0.0_8
        tempd = tmpd17/(fxlift+fylift+fzlift)
        funcvaluesd(costfunccoforcexz) = funcvaluesd(costfunccoforcexz) &
&         + fxlift*tempd
        funcvaluesd(costfunccoforceyz) = funcvaluesd(costfunccoforceyz) &
&         + fylift*tempd
        funcvaluesd(costfunccoforcezz) = funcvaluesd(costfunccoforcezz) &
&         + fzlift*tempd
        tempd0 = -((fxlift*funcvalues(costfunccoforcexz)+fylift*&
&         funcvalues(costfunccoforceyz)+fzlift*funcvalues(&
&         costfunccoforcezz))*tempd/(fxlift+fylift+fzlift))
        fxliftd = funcvalues(costfunccoforcexz)*tempd + tempd0
        fyliftd = funcvalues(costfunccoforceyz)*tempd + tempd0
        fzliftd = funcvalues(costfunccoforcezz)*tempd + tempd0
        call popreal8(funcvalues(costfunccoflifty))
        tmpd16 = funcvaluesd(costfunccoflifty)
        funcvaluesd(costfunccoflifty) = 0.0_8
        tempd = tmpd16/(fxlift+fylift+fzlift)
        funcvaluesd(costfunccoforcexy) = funcvaluesd(costfunccoforcexy) &
&         + fxlift*tempd
        funcvaluesd(costfunccoforceyy) = funcvaluesd(costfunccoforceyy) &
&         + fylift*tempd
        funcvaluesd(costfunccoforcezy) = funcvaluesd(costfunccoforcezy) &
&         + fzlift*tempd
        tempd0 = -((fxlift*funcvalues(costfunccoforcexy)+fylift*&
&         funcvalues(costfunccoforceyy)+fzlift*funcvalues(&
&         costfunccoforcezy))*tempd/(fxlift+fylift+fzlift))
        fxliftd = fxliftd + funcvalues(costfunccoforcexy)*tempd + tempd0
        fyliftd = fyliftd + funcvalues(costfunccoforceyy)*tempd + tempd0
        fzliftd = fzliftd + funcvalues(costfunccoforcezy)*tempd + tempd0
        call popreal8(funcvalues(costfunccofliftx))
        tmpd15 = funcvaluesd(costfunccofliftx)
        funcvaluesd(costfunccofliftx) = 0.0_8
        tempd = tmpd15/(fxlift+fylift+fzlift)
        funcvaluesd(costfunccoforcexx) = funcvaluesd(costfunccoforcexx) &
&         + fxlift*tempd
        funcvaluesd(costfunccoforceyx) = funcvaluesd(costfunccoforceyx) &
&         + fylift*tempd
        funcvaluesd(costfunccoforcezx) = funcvaluesd(costfunccoforcezx) &
&         + fzlift*tempd
        tempd0 = -((fxlift*funcvalues(costfunccoforcexx)+fylift*&
&         funcvalues(costfunccoforceyx)+fzlift*funcvalues(&
&         costfunccoforcezx))*tempd/(fxlift+fylift+fzlift))
        fxliftd = fxliftd + funcvalues(costfunccoforcexx)*tempd + tempd0
        fyliftd = fyliftd + funcvalues(costfunccoforceyx)*tempd + tempd0
        fzliftd = fzliftd + funcvalues(costfunccoforcezx)*tempd + tempd0
      else
        funcvaluesd(costfunccofliftz) = 0.0_8
        funcvaluesd(costfunccoflifty) = 0.0_8
        funcvaluesd(costfunccofliftx) = 0.0_8
        fxliftd = 0.0_8
        fyliftd = 0.0_8
        fzliftd = 0.0_8
      end if
      funcvaluesd(costfuncforcez) = funcvaluesd(costfuncforcez) + &
&       liftdirection(3)*fzliftd
      liftdirectiond(3) = liftdirectiond(3) + funcvalues(costfuncforcez)&
&       *fzliftd
      funcvaluesd(costfuncforcey) = funcvaluesd(costfuncforcey) + &
&       liftdirection(2)*fyliftd
      liftdirectiond(2) = liftdirectiond(2) + funcvalues(costfuncforcey)&
&       *fyliftd
      funcvaluesd(costfuncforcex) = funcvaluesd(costfuncforcex) + &
&       liftdirection(1)*fxliftd
      liftdirectiond(1) = liftdirectiond(1) + funcvalues(costfuncforcex)&
&       *fxliftd
      call popreal8(funcvalues(costfuncdragcoefmomentum))
      tmpd14 = funcvaluesd(costfuncdragcoefmomentum)
      funcvaluesd(costfuncdragcoefmomentum) = 0.0_8
      funcvaluesd(costfuncforcexcoefmomentum) = funcvaluesd(&
&       costfuncforcexcoefmomentum) + dragdirection(1)*tmpd14
      dragdirectiond(1) = dragdirectiond(1) + funcvalues(&
&       costfuncforcexcoefmomentum)*tmpd14
      funcvaluesd(costfuncforceycoefmomentum) = funcvaluesd(&
&       costfuncforceycoefmomentum) + dragdirection(2)*tmpd14
      dragdirectiond(2) = dragdirectiond(2) + funcvalues(&
&       costfuncforceycoefmomentum)*tmpd14
      funcvaluesd(costfuncforcezcoefmomentum) = funcvaluesd(&
&       costfuncforcezcoefmomentum) + dragdirection(3)*tmpd14
      dragdirectiond(3) = dragdirectiond(3) + funcvalues(&
&       costfuncforcezcoefmomentum)*tmpd14
      call popreal8(funcvalues(costfuncdragcoefviscous))
      tmpd13 = funcvaluesd(costfuncdragcoefviscous)
      funcvaluesd(costfuncdragcoefviscous) = 0.0_8
      funcvaluesd(costfuncforcexcoefviscous) = funcvaluesd(&
&       costfuncforcexcoefviscous) + dragdirection(1)*tmpd13
      dragdirectiond(1) = dragdirectiond(1) + funcvalues(&
&       costfuncforcexcoefviscous)*tmpd13
      funcvaluesd(costfuncforceycoefviscous) = funcvaluesd(&
&       costfuncforceycoefviscous) + dragdirection(2)*tmpd13
      dragdirectiond(2) = dragdirectiond(2) + funcvalues(&
&       costfuncforceycoefviscous)*tmpd13
      funcvaluesd(costfuncforcezcoefviscous) = funcvaluesd(&
&       costfuncforcezcoefviscous) + dragdirection(3)*tmpd13
      dragdirectiond(3) = dragdirectiond(3) + funcvalues(&
&       costfuncforcezcoefviscous)*tmpd13
      call popreal8(funcvalues(costfuncdragcoefpressure))
      tmpd12 = funcvaluesd(costfuncdragcoefpressure)
      funcvaluesd(costfuncdragcoefpressure) = 0.0_8
      funcvaluesd(costfuncforcexcoefpressure) = funcvaluesd(&
&       costfuncforcexcoefpressure) + dragdirection(1)*tmpd12
      dragdirectiond(1) = dragdirectiond(1) + funcvalues(&
&       costfuncforcexcoefpressure)*tmpd12
      funcvaluesd(costfuncforceycoefpressure) = funcvaluesd(&
&       costfuncforceycoefpressure) + dragdirection(2)*tmpd12
      dragdirectiond(2) = dragdirectiond(2) + funcvalues(&
&       costfuncforceycoefpressure)*tmpd12
      funcvaluesd(costfuncforcezcoefpressure) = funcvaluesd(&
&       costfuncforcezcoefpressure) + dragdirection(3)*tmpd12
      dragdirectiond(3) = dragdirectiond(3) + funcvalues(&
&       costfuncforcezcoefpressure)*tmpd12
      call popreal8(funcvalues(costfuncdragcoef))
      tmpd11 = funcvaluesd(costfuncdragcoef)
      funcvaluesd(costfuncdragcoef) = 0.0_8
      funcvaluesd(costfuncforcexcoef) = funcvaluesd(costfuncforcexcoef) &
&       + dragdirection(1)*tmpd11
      dragdirectiond(1) = dragdirectiond(1) + funcvalues(&
&       costfuncforcexcoef)*tmpd11
      funcvaluesd(costfuncforceycoef) = funcvaluesd(costfuncforceycoef) &
&       + dragdirection(2)*tmpd11
      dragdirectiond(2) = dragdirectiond(2) + funcvalues(&
&       costfuncforceycoef)*tmpd11
      funcvaluesd(costfuncforcezcoef) = funcvaluesd(costfuncforcezcoef) &
&       + dragdirection(3)*tmpd11
      dragdirectiond(3) = dragdirectiond(3) + funcvalues(&
&       costfuncforcezcoef)*tmpd11
      call popreal8(funcvalues(costfuncliftcoefmomentum))
      tmpd10 = funcvaluesd(costfuncliftcoefmomentum)
      funcvaluesd(costfuncliftcoefmomentum) = 0.0_8
      funcvaluesd(costfuncforcexcoefmomentum) = funcvaluesd(&
&       costfuncforcexcoefmomentum) + liftdirection(1)*tmpd10
      liftdirectiond(1) = liftdirectiond(1) + funcvalues(&
&       costfuncforcexcoefmomentum)*tmpd10
      funcvaluesd(costfuncforceycoefmomentum) = funcvaluesd(&
&       costfuncforceycoefmomentum) + liftdirection(2)*tmpd10
      liftdirectiond(2) = liftdirectiond(2) + funcvalues(&
&       costfuncforceycoefmomentum)*tmpd10
      funcvaluesd(costfuncforcezcoefmomentum) = funcvaluesd(&
&       costfuncforcezcoefmomentum) + liftdirection(3)*tmpd10
      liftdirectiond(3) = liftdirectiond(3) + funcvalues(&
&       costfuncforcezcoefmomentum)*tmpd10
      call popreal8(funcvalues(costfuncliftcoefviscous))
      tmpd9 = funcvaluesd(costfuncliftcoefviscous)
      funcvaluesd(costfuncliftcoefviscous) = 0.0_8
      funcvaluesd(costfuncforcexcoefviscous) = funcvaluesd(&
&       costfuncforcexcoefviscous) + liftdirection(1)*tmpd9
      liftdirectiond(1) = liftdirectiond(1) + funcvalues(&
&       costfuncforcexcoefviscous)*tmpd9
      funcvaluesd(costfuncforceycoefviscous) = funcvaluesd(&
&       costfuncforceycoefviscous) + liftdirection(2)*tmpd9
      liftdirectiond(2) = liftdirectiond(2) + funcvalues(&
&       costfuncforceycoefviscous)*tmpd9
      funcvaluesd(costfuncforcezcoefviscous) = funcvaluesd(&
&       costfuncforcezcoefviscous) + liftdirection(3)*tmpd9
      liftdirectiond(3) = liftdirectiond(3) + funcvalues(&
&       costfuncforcezcoefviscous)*tmpd9
      call popreal8(funcvalues(costfuncliftcoefpressure))
      tmpd8 = funcvaluesd(costfuncliftcoefpressure)
      funcvaluesd(costfuncliftcoefpressure) = 0.0_8
      funcvaluesd(costfuncforcexcoefpressure) = funcvaluesd(&
&       costfuncforcexcoefpressure) + liftdirection(1)*tmpd8
      liftdirectiond(1) = liftdirectiond(1) + funcvalues(&
&       costfuncforcexcoefpressure)*tmpd8
      funcvaluesd(costfuncforceycoefpressure) = funcvaluesd(&
&       costfuncforceycoefpressure) + liftdirection(2)*tmpd8
      liftdirectiond(2) = liftdirectiond(2) + funcvalues(&
&       costfuncforceycoefpressure)*tmpd8
      funcvaluesd(costfuncforcezcoefpressure) = funcvaluesd(&
&       costfuncforcezcoefpressure) + liftdirection(3)*tmpd8
      liftdirectiond(3) = liftdirectiond(3) + funcvalues(&
&       costfuncforcezcoefpressure)*tmpd8
      call popreal8(funcvalues(costfuncliftcoef))
      tmpd7 = funcvaluesd(costfuncliftcoef)
      funcvaluesd(costfuncliftcoef) = 0.0_8
      funcvaluesd(costfuncforcexcoef) = funcvaluesd(costfuncforcexcoef) &
&       + liftdirection(1)*tmpd7
      liftdirectiond(1) = liftdirectiond(1) + funcvalues(&
&       costfuncforcexcoef)*tmpd7
      funcvaluesd(costfuncforceycoef) = funcvaluesd(costfuncforceycoef) &
&       + liftdirection(2)*tmpd7
      liftdirectiond(2) = liftdirectiond(2) + funcvalues(&
&       costfuncforceycoef)*tmpd7
      funcvaluesd(costfuncforcezcoef) = funcvaluesd(costfuncforcezcoef) &
&       + liftdirection(3)*tmpd7
      liftdirectiond(3) = liftdirectiond(3) + funcvalues(&
&       costfuncforcezcoef)*tmpd7
      call popreal8(funcvalues(costfuncdragmomentum))
      tmpd6 = funcvaluesd(costfuncdragmomentum)
      funcvaluesd(costfuncdragmomentum) = 0.0_8
      funcvaluesd(costfuncforcexmomentum) = funcvaluesd(&
&       costfuncforcexmomentum) + dragdirection(1)*tmpd6
      dragdirectiond(1) = dragdirectiond(1) + funcvalues(&
&       costfuncforcexmomentum)*tmpd6
      funcvaluesd(costfuncforceymomentum) = funcvaluesd(&
&       costfuncforceymomentum) + dragdirection(2)*tmpd6
      dragdirectiond(2) = dragdirectiond(2) + funcvalues(&
&       costfuncforceymomentum)*tmpd6
      funcvaluesd(costfuncforcezmomentum) = funcvaluesd(&
&       costfuncforcezmomentum) + dragdirection(3)*tmpd6
      dragdirectiond(3) = dragdirectiond(3) + funcvalues(&
&       costfuncforcezmomentum)*tmpd6
      call popreal8(funcvalues(costfuncdragviscous))
      tmpd5 = funcvaluesd(costfuncdragviscous)
      funcvaluesd(costfuncdragviscous) = 0.0_8
      funcvaluesd(costfuncforcexviscous) = funcvaluesd(&
&       costfuncforcexviscous) + dragdirection(1)*tmpd5
      dragdirectiond(1) = dragdirectiond(1) + funcvalues(&
&       costfuncforcexviscous)*tmpd5
      funcvaluesd(costfuncforceyviscous) = funcvaluesd(&
&       costfuncforceyviscous) + dragdirection(2)*tmpd5
      dragdirectiond(2) = dragdirectiond(2) + funcvalues(&
&       costfuncforceyviscous)*tmpd5
      funcvaluesd(costfuncforcezviscous) = funcvaluesd(&
&       costfuncforcezviscous) + dragdirection(3)*tmpd5
      dragdirectiond(3) = dragdirectiond(3) + funcvalues(&
&       costfuncforcezviscous)*tmpd5
      call popreal8(funcvalues(costfuncdragpressure))
      tmpd4 = funcvaluesd(costfuncdragpressure)
      funcvaluesd(costfuncdragpressure) = 0.0_8
      funcvaluesd(costfuncforcexpressure) = funcvaluesd(&
&       costfuncforcexpressure) + dragdirection(1)*tmpd4
      dragdirectiond(1) = dragdirectiond(1) + funcvalues(&
&       costfuncforcexpressure)*tmpd4
      funcvaluesd(costfuncforceypressure) = funcvaluesd(&
&       costfuncforceypressure) + dragdirection(2)*tmpd4
      dragdirectiond(2) = dragdirectiond(2) + funcvalues(&
&       costfuncforceypressure)*tmpd4
      funcvaluesd(costfuncforcezpressure) = funcvaluesd(&
&       costfuncforcezpressure) + dragdirection(3)*tmpd4
      dragdirectiond(3) = dragdirectiond(3) + funcvalues(&
&       costfuncforcezpressure)*tmpd4
      call popreal8(funcvalues(costfuncdrag))
      tmpd3 = funcvaluesd(costfuncdrag)
      funcvaluesd(costfuncdrag) = 0.0_8
      funcvaluesd(costfuncforcex) = funcvaluesd(costfuncforcex) + &
&       dragdirection(1)*tmpd3
      dragdirectiond(1) = dragdirectiond(1) + funcvalues(costfuncforcex)&
&       *tmpd3
      funcvaluesd(costfuncforcey) = funcvaluesd(costfuncforcey) + &
&       dragdirection(2)*tmpd3
      dragdirectiond(2) = dragdirectiond(2) + funcvalues(costfuncforcey)&
&       *tmpd3
      funcvaluesd(costfuncforcez) = funcvaluesd(costfuncforcez) + &
&       dragdirection(3)*tmpd3
      dragdirectiond(3) = dragdirectiond(3) + funcvalues(costfuncforcez)&
&       *tmpd3
      call popreal8(funcvalues(costfuncliftmomentum))
      tmpd2 = funcvaluesd(costfuncliftmomentum)
      funcvaluesd(costfuncliftmomentum) = 0.0_8
      funcvaluesd(costfuncforcexmomentum) = funcvaluesd(&
&       costfuncforcexmomentum) + liftdirection(1)*tmpd2
      liftdirectiond(1) = liftdirectiond(1) + funcvalues(&
&       costfuncforcexmomentum)*tmpd2
      funcvaluesd(costfuncforceymomentum) = funcvaluesd(&
&       costfuncforceymomentum) + liftdirection(2)*tmpd2
      liftdirectiond(2) = liftdirectiond(2) + funcvalues(&
&       costfuncforceymomentum)*tmpd2
      funcvaluesd(costfuncforcezmomentum) = funcvaluesd(&
&       costfuncforcezmomentum) + liftdirection(3)*tmpd2
      liftdirectiond(3) = liftdirectiond(3) + funcvalues(&
&       costfuncforcezmomentum)*tmpd2
      call popreal8(funcvalues(costfuncliftviscous))
      tmpd1 = funcvaluesd(costfuncliftviscous)
      funcvaluesd(costfuncliftviscous) = 0.0_8
      funcvaluesd(costfuncforcexviscous) = funcvaluesd(&
&       costfuncforcexviscous) + liftdirection(1)*tmpd1
      liftdirectiond(1) = liftdirectiond(1) + funcvalues(&
&       costfuncforcexviscous)*tmpd1
      funcvaluesd(costfuncforceyviscous) = funcvaluesd(&
&       costfuncforceyviscous) + liftdirection(2)*tmpd1
      liftdirectiond(2) = liftdirectiond(2) + funcvalues(&
&       costfuncforceyviscous)*tmpd1
      funcvaluesd(costfuncforcezviscous) = funcvaluesd(&
&       costfuncforcezviscous) + liftdirection(3)*tmpd1
      liftdirectiond(3) = liftdirectiond(3) + funcvalues(&
&       costfuncforcezviscous)*tmpd1
      call popreal8(funcvalues(costfuncliftpressure))
      tmpd0 = funcvaluesd(costfuncliftpressure)
      funcvaluesd(costfuncliftpressure) = 0.0_8
      funcvaluesd(costfuncforcexpressure) = funcvaluesd(&
&       costfuncforcexpressure) + liftdirection(1)*tmpd0
      liftdirectiond(1) = liftdirectiond(1) + funcvalues(&
&       costfuncforcexpressure)*tmpd0
      funcvaluesd(costfuncforceypressure) = funcvaluesd(&
&       costfuncforceypressure) + liftdirection(2)*tmpd0
      liftdirectiond(2) = liftdirectiond(2) + funcvalues(&
&       costfuncforceypressure)*tmpd0
      funcvaluesd(costfuncforcezpressure) = funcvaluesd(&
&       costfuncforcezpressure) + liftdirection(3)*tmpd0
      liftdirectiond(3) = liftdirectiond(3) + funcvalues(&
&       costfuncforcezpressure)*tmpd0
      call popreal8(funcvalues(costfunclift))
      tmpd = funcvaluesd(costfunclift)
      funcvaluesd(costfunclift) = 0.0_8
      funcvaluesd(costfuncforcex) = funcvaluesd(costfuncforcex) + &
&       liftdirection(1)*tmpd
      liftdirectiond(1) = liftdirectiond(1) + funcvalues(costfuncforcex)&
&       *tmpd
      funcvaluesd(costfuncforcey) = funcvaluesd(costfuncforcey) + &
&       liftdirection(2)*tmpd
      liftdirectiond(2) = liftdirectiond(2) + funcvalues(costfuncforcey)&
&       *tmpd
      funcvaluesd(costfuncforcez) = funcvaluesd(costfuncforcez) + &
&       liftdirection(3)*tmpd
      liftdirectiond(3) = liftdirectiond(3) + funcvalues(costfuncforcez)&
&       *tmpd
      globalvalsd = 0.0_8
      momentd = 0.0_8
      cforced = 0.0_8
      cofxd = 0.0_8
      cofyd = 0.0_8
      cofzd = 0.0_8
      forced = 0.0_8
      cforcemd = 0.0_8
      forcemd = 0.0_8
      cforcepd = 0.0_8
      forcepd = 0.0_8
      cforcevd = 0.0_8
      cmomentd = 0.0_8
      forcevd = 0.0_8
      call popreal8array(cofz, 3*ntimeintervalsspectral)
      call popreal8array(cofy, 3*ntimeintervalsspectral)
      call popreal8array(cofx, 3*ntimeintervalsspectral)
!$bwd-of ii-loop 
      do sps=1,ntimeintervalsspectral
! ------------
! ------------
! 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) .ne. zero) then
          call pushcontrol1b(0)
        else
          call pushcontrol1b(1)
        end if
        if (force(2, sps) .ne. zero) then
          call pushcontrol1b(0)
        else
          call pushcontrol1b(1)
        end if
        if (force(3, sps) .ne. zero) then
          call pushcontrol1b(0)
        else
          call pushcontrol1b(1)
        end if
! fx
! fy
! fz
! ------------
! 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)
          call pushcontrol1b(0)
        else
          call pushcontrol1b(1)
        end if
! final part of the ks computation
        if (computecavitation) then
          call pushcontrol1b(0)
        else
          call pushcontrol1b(1)
        end if
! mass flow like objective
        mflow = globalvals(imassflow, sps)
        if (mflow .ne. zero) then
          call pushcontrol1b(0)
        else
          call pushcontrol1b(1)
        end if
! area averaged objectives
        garea = globalvals(iarea, sps)
        if (garea .ne. zero) then
          call pushcontrol1b(0)
        else
          call pushcontrol1b(1)
        end if
        mavgvid = ovrnts*funcvaluesd(costfuncmavgvi)
        mavgvzd = ovrnts*funcvaluesd(costfuncmavgvz)
        mavgvyd = ovrnts*funcvaluesd(costfuncmavgvy)
        mavgvxd = ovrnts*funcvaluesd(costfuncmavgvx)
        mavgad = ovrnts*funcvaluesd(costfuncmavga)
        mavgmnd = ovrnts*funcvaluesd(costfuncmavgmn)
        mavgpsd = ovrnts*funcvaluesd(costfuncmavgps)
        mavgrhod = ovrnts*funcvaluesd(costfuncmavgrho)
        mavgttotd = ovrnts*funcvaluesd(costfuncmavgttot)
        mavgptotd = ovrnts*funcvaluesd(costfuncmavgptot)
        mflowd = ovrnts*funcvaluesd(costfuncmdot)
        call popcontrol1b(branch)
        if (branch .eq. 0) then
          tempd = ovrnts*funcvaluesd(costfuncaavgps)/garea
          globalvalsd(iareaps, sps) = globalvalsd(iareaps, sps) + tempd
          garead = -(globalvals(iareaps, sps)*tempd/garea)
          tempd = ovrnts*funcvaluesd(costfuncaavgptot)/garea
          globalvalsd(iareaptot, sps) = globalvalsd(iareaptot, sps) + &
&           tempd
          garead = garead - globalvals(iareaptot, sps)*tempd/garea
        else
          garead = 0.0_8
        end if
        globalvalsd(iarea, sps) = globalvalsd(iarea, sps) + garead
        call popcontrol1b(branch)
        if (branch .eq. 0) then
          globalvalsd(imassvi, sps) = globalvalsd(imassvi, sps) + &
&           mavgvid/mflow
          mflowd = mflowd - globalvals(imassvi, sps)*mavgvid/mflow**2 - &
&           globalvals(imassvz, sps)*mavgvzd/mflow**2 - globalvals(&
&           imassvy, sps)*mavgvyd/mflow**2 - globalvals(imassvx, sps)*&
&           mavgvxd/mflow**2 - globalvals(imassa, sps)*mavgad/mflow**2 -&
&           globalvals(imassmn, sps)*mavgmnd/mflow**2 - globalvals(&
&           imassps, sps)*mavgpsd/mflow**2 - globalvals(imassrho, sps)*&
&           mavgrhod/mflow**2 - globalvals(imassttot, sps)*mavgttotd/&
&           mflow**2 - globalvals(imassptot, sps)*mavgptotd/mflow**2
          globalvalsd(imassvz, sps) = globalvalsd(imassvz, sps) + &
&           mavgvzd/mflow
          globalvalsd(imassvy, sps) = globalvalsd(imassvy, sps) + &
&           mavgvyd/mflow
          globalvalsd(imassvx, sps) = globalvalsd(imassvx, sps) + &
&           mavgvxd/mflow
          globalvalsd(imassa, sps) = globalvalsd(imassa, sps) + mavgad/&
&           mflow
          globalvalsd(imassmn, sps) = globalvalsd(imassmn, sps) + &
&           mavgmnd/mflow
          globalvalsd(imassps, sps) = globalvalsd(imassps, sps) + &
&           mavgpsd/mflow
          globalvalsd(imassrho, sps) = globalvalsd(imassrho, sps) + &
&           mavgrhod/mflow
          globalvalsd(imassttot, sps) = globalvalsd(imassttot, sps) + &
&           mavgttotd/mflow
          globalvalsd(imassptot, sps) = globalvalsd(imassptot, sps) + &
&           mavgptotd/mflow
        end if
        globalvalsd(imassflow, sps) = globalvalsd(imassflow, sps) + &
&         mflowd
        globalvalsd(icperror2, sps) = globalvalsd(icperror2, sps) + &
&         ovrnts*funcvaluesd(costfunccperror2)
        globalvalsd(ipower, sps) = globalvalsd(ipower, sps) + ovrnts*&
&         funcvaluesd(costfuncflowpower)
        globalvalsd(iarea, sps) = globalvalsd(iarea, sps) + ovrnts*&
&         funcvaluesd(costfuncarea)
        globalvalsd(isepavg+2, sps) = globalvalsd(isepavg+2, sps) + &
&         ovrnts*funcvaluesd(costfuncsepsensoravgz)
        globalvalsd(isepavg+1, sps) = globalvalsd(isepavg+1, sps) + &
&         ovrnts*funcvaluesd(costfuncsepsensoravgy)
        globalvalsd(isepavg, sps) = globalvalsd(isepavg, sps) + ovrnts*&
&         funcvaluesd(costfuncsepsensoravgx)
        globalvalsd(iaxismoment, sps) = globalvalsd(iaxismoment, sps) + &
&         ovrnts*funcvaluesd(costfuncaxismoment)
        call popcontrol1b(branch)
        if (branch .eq. 0) globalvalsd(icpmin, sps) = globalvalsd(icpmin&
&           , sps) - ovrnts*funcvaluesd(costfunccpmin)/(globalvals(&
&           icpmin, sps)*cpmin_rho)
        globalvalsd(icavitation, sps) = globalvalsd(icavitation, sps) + &
&         ovrnts*funcvaluesd(costfunccavitation)
        globalvalsd(isepsensor, sps) = globalvalsd(isepsensor, sps) + &
&         ovrnts*funcvaluesd(costfuncsepsensor)
        call popcontrol1b(branch)
        if (branch .eq. 0) then
          temp1 = 2*sepsensorkssharpness*(sepsensorksoffset+ks_comp)
          temp0 = one + exp(temp1)
          tempd1 = ovrnts*one*funcvaluesd(costfuncsepsensorksarea)/temp0
          globalvalsd(isepsensorarea, sps) = globalvalsd(isepsensorarea&
&           , sps) + ovrnts*funcvaluesd(costfuncsepsensorksarea)
          globalvalsd(isepsensorksarea, sps) = globalvalsd(&
&           isepsensorksarea, sps) + ks_comp*tempd1
          ks_compd = (globalvals(isepsensorksarea, sps)-&
&           sepsensorkssharpness*2*exp(temp1)*globalvals(&
&           isepsensorksarea, sps)*ks_comp/temp0)*tempd1 + funcvaluesd(&
&           costfuncsepsensorks)
          globalvalsd(isepsensorks, sps) = globalvalsd(isepsensorks, sps&
&           ) + ovrnts*ks_compd/(globalvals(isepsensorks, sps)*&
&           sepsenmaxrho)
        end if
        cmomentd(3, sps) = cmomentd(3, sps) + ovrnts*funcvaluesd(&
&         costfuncmomzcoef)
        cmomentd(2, sps) = cmomentd(2, sps) + ovrnts*funcvaluesd(&
&         costfuncmomycoef)
        cmomentd(1, sps) = cmomentd(1, sps) + ovrnts*funcvaluesd(&
&         costfuncmomxcoef)
        momentd(3, sps) = momentd(3, sps) + ovrnts*funcvaluesd(&
&         costfuncmomz)
        momentd(2, sps) = momentd(2, sps) + ovrnts*funcvaluesd(&
&         costfuncmomy)
        momentd(1, sps) = momentd(1, sps) + ovrnts*funcvaluesd(&
&         costfuncmomx)
        cofzd(3, sps) = cofzd(3, sps) + ovrnts*funcvaluesd(&
&         costfunccoforcezz)
        cofzd(2, sps) = cofzd(2, sps) + ovrnts*funcvaluesd(&
&         costfunccoforcezy)
        cofzd(1, sps) = cofzd(1, sps) + ovrnts*funcvaluesd(&
&         costfunccoforcezx)
        cofyd(3, sps) = cofyd(3, sps) + ovrnts*funcvaluesd(&
&         costfunccoforceyz)
        cofyd(2, sps) = cofyd(2, sps) + ovrnts*funcvaluesd(&
&         costfunccoforceyy)
        cofyd(1, sps) = cofyd(1, sps) + ovrnts*funcvaluesd(&
&         costfunccoforceyx)
        cofxd(3, sps) = cofxd(3, sps) + ovrnts*funcvaluesd(&
&         costfunccoforcexz)
        cofxd(2, sps) = cofxd(2, sps) + ovrnts*funcvaluesd(&
&         costfunccoforcexy)
        cofxd(1, sps) = cofxd(1, sps) + ovrnts*funcvaluesd(&
&         costfunccoforcexx)
        call popcontrol1b(branch)
        if (branch .eq. 0) then
          forced(3, sps) = forced(3, sps) - sum(cofz(:, sps)*cofzd(:, &
&           sps))/force(3, sps)**2
          cofzd(:, sps) = cofzd(:, sps)/force(3, sps)
        else
          cofzd(:, sps) = 0.0_8
        end if
        call popcontrol1b(branch)
        if (branch .eq. 0) then
          forced(2, sps) = forced(2, sps) - sum(cofy(:, sps)*cofyd(:, &
&           sps))/force(2, sps)**2
          cofyd(:, sps) = cofyd(:, sps)/force(2, sps)
        else
          cofyd(:, sps) = 0.0_8
        end if
        call popcontrol1b(branch)
        if (branch .eq. 0) then
          forced(1, sps) = forced(1, sps) - sum(cofx(:, sps)*cofxd(:, &
&           sps))/force(1, sps)**2
          cofxd(:, sps) = cofxd(:, sps)/force(1, sps)
        else
          cofxd(:, sps) = 0.0_8
        end if
        cforcemd(3, sps) = cforcemd(3, sps) + ovrnts*funcvaluesd(&
&         costfuncforcezcoefmomentum)
        cforcemd(2, sps) = cforcemd(2, sps) + ovrnts*funcvaluesd(&
&         costfuncforceycoefmomentum)
        cforcemd(1, sps) = cforcemd(1, sps) + ovrnts*funcvaluesd(&
&         costfuncforcexcoefmomentum)
        cforcevd(3, sps) = cforcevd(3, sps) + ovrnts*funcvaluesd(&
&         costfuncforcezcoefviscous)
        cforcevd(2, sps) = cforcevd(2, sps) + ovrnts*funcvaluesd(&
&         costfuncforceycoefviscous)
        cforcevd(1, sps) = cforcevd(1, sps) + ovrnts*funcvaluesd(&
&         costfuncforcexcoefviscous)
        cforcepd(3, sps) = cforcepd(3, sps) + ovrnts*funcvaluesd(&
&         costfuncforcezcoefpressure)
        cforcepd(2, sps) = cforcepd(2, sps) + ovrnts*funcvaluesd(&
&         costfuncforceycoefpressure)
        cforcepd(1, sps) = cforcepd(1, sps) + ovrnts*funcvaluesd(&
&         costfuncforcexcoefpressure)
        cforced(3, sps) = cforced(3, sps) + ovrnts*funcvaluesd(&
&         costfuncforcezcoef)
        cforced(2, sps) = cforced(2, sps) + ovrnts*funcvaluesd(&
&         costfuncforceycoef)
        cforced(1, sps) = cforced(1, sps) + ovrnts*funcvaluesd(&
&         costfuncforcexcoef)
        forcemd(3, sps) = forcemd(3, sps) + ovrnts*funcvaluesd(&
&         costfuncforcezmomentum)
        forcemd(2, sps) = forcemd(2, sps) + ovrnts*funcvaluesd(&
&         costfuncforceymomentum)
        forcemd(1, sps) = forcemd(1, sps) + ovrnts*funcvaluesd(&
&         costfuncforcexmomentum)
        forcevd(3, sps) = forcevd(3, sps) + ovrnts*funcvaluesd(&
&         costfuncforcezviscous)
        forcevd(2, sps) = forcevd(2, sps) + ovrnts*funcvaluesd(&
&         costfuncforceyviscous)
        forcevd(1, sps) = forcevd(1, sps) + ovrnts*funcvaluesd(&
&         costfuncforcexviscous)
        forcepd(3, sps) = forcepd(3, sps) + ovrnts*funcvaluesd(&
&         costfuncforcezpressure)
        forcepd(2, sps) = forcepd(2, sps) + ovrnts*funcvaluesd(&
&         costfuncforceypressure)
        forcepd(1, sps) = forcepd(1, sps) + ovrnts*funcvaluesd(&
&         costfuncforcexpressure)
        forced(3, sps) = forced(3, sps) + ovrnts*funcvaluesd(&
&         costfuncforcez)
        forced(2, sps) = forced(2, sps) + ovrnts*funcvaluesd(&
&         costfuncforcey)
        forced(1, sps) = forced(1, sps) + ovrnts*funcvaluesd(&
&         costfuncforcex)
      end do
      factd = sum(moment*cmomentd)
      momentd = momentd + fact*cmomentd
      call popreal8(fact)
      factd = factd/(lengthref*lref) + sum(forcem*cforcemd) + sum(forcev&
&       *cforcevd) + sum(forcep*cforcepd) + sum(force*cforced)
      forcemd = forcemd + fact*cforcemd
      forcevd = forcevd + fact*cforcevd
      forcepd = forcepd + fact*cforcepd
      forced = forced + fact*cforced
      temp = gammainf*surfaceref*(lref*lref)
      temp0 = temp*(machcoef*machcoef)*pref
      tempd = -(temp*two*factd/temp0**2)
      machcoefd = machcoefd + 2*machcoef*pref*tempd
      prefd = prefd + machcoef**2*tempd
      globalvalsd(imp:imp+2, :) = globalvalsd(imp:imp+2, :) + momentd
      globalvalsd(imv:imv+2, :) = globalvalsd(imv:imv+2, :) + momentd
      globalvalsd(iflowmm:iflowmm+2, :) = globalvalsd(iflowmm:iflowmm+2&
&       , :) + momentd
      globalvalsd(icoforcez:icoforcez+2, :) = globalvalsd(icoforcez:&
&       icoforcez+2, :) + cofzd
      globalvalsd(icoforcey:icoforcey+2, :) = globalvalsd(icoforcey:&
&       icoforcey+2, :) + cofyd
      globalvalsd(icoforcex:icoforcex+2, :) = globalvalsd(icoforcex:&
&       icoforcex+2, :) + cofxd
      globalvalsd(iflowfm:iflowfm+2, :) = globalvalsd(iflowfm:iflowfm+2&
&       , :) + forcemd
      globalvalsd(ifv:ifv+2, :) = globalvalsd(ifv:ifv+2, :) + forcevd
      globalvalsd(ifp:ifp+2, :) = globalvalsd(ifp:ifp+2, :) + forcepd + &
&       forced
      globalvalsd(ifv:ifv+2, :) = globalvalsd(ifv:ifv+2, :) + forced
      globalvalsd(iflowfm:iflowfm+2, :) = globalvalsd(iflowfm:iflowfm+2&
&       , :) + forced
      funcvaluesd = 0.0_8
    end if
  end subroutine getcostfunctions_b
 
  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_b, only : computetsderivatives
    use flowutils_b, only : getdirvector
    implicit none
! input/output
    real(kind=realtype), dimension(:, :), intent(in) :: globalvals
    real(kind=realtype), dimension(:), intent(out) :: 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
    intrinsic log
    intrinsic exp
    intrinsic sqrt
! 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
!$ad ii-loop
! here we finally assign the final function values
    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) .ne. zero) then
        cofx(:, sps) = cofx(:, sps)/force(1, sps)
      else
        cofx(:, sps) = zero
      end if
      if (force(2, sps) .ne. zero) then
        cofy(:, sps) = cofy(:, sps)/force(2, sps)
      else
        cofy(:, sps) = zero
      end if
      if (force(3, sps) .ne. 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 .ne. 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 .ne. 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 .ne. 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
    end if
  end subroutine getcostfunctions
 
!  differentiation of wallintegrationface in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: pinf pref *xx *pp1 *pp2 *ssi
!                *ww2 veldirfreestream machcoef pointref *(*viscsubface.tau)
!                *(*bcdata.fv) *(*bcdata.fp) *(*bcdata.area) localvalues
!   with respect to varying inputs: pinf pref *xx *pp1 *pp2 *ssi
!                *ww2 veldirfreestream machcoef pointref *(*viscsubface.tau)
!                *(*bcdata.fv) *(*bcdata.fp) *(*bcdata.area) localvalues
!   rw status of diff variables: pinf:incr pref:incr *xx:incr *pp1:incr
!                *pp2:incr *ssi:incr *ww2:incr veldirfreestream:incr
!                machcoef:incr pointref:incr *(*viscsubface.tau):incr
!                *(*bcdata.fv):in-out *(*bcdata.fp):in-out *(*bcdata.area):in-out
!                localvalues:in-out
!   plus diff mem management of: xx:in pp1:in pp2:in ssi:in ww2:in
!                viscsubface:in *viscsubface.tau:in bcdata:in *bcdata.fv:in
!                *bcdata.fp:in *bcdata.area:in
  subroutine wallintegrationface_b(localvalues, localvaluesd, 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, machcoefd, pointref, pointrefd,&
&   veldirfreestream, veldirfreestreamd, equations, momentaxis, &
&   cpmin_family, cpmin_rho, cavitationnumber, sepsenmaxfamily, &
&   sepsenmaxrho
    use bcpointers_b
    implicit none
! input/output variables
    real(kind=realtype), dimension(nlocalvalues), intent(inout) :: &
&   localvalues
    real(kind=realtype), dimension(nlocalvalues), intent(inout) :: &
&   localvaluesd
    integer(kind=inttype) :: mm
! local variables.
    real(kind=realtype), dimension(3) :: fp, fv, mp, mv
    real(kind=realtype), dimension(3) :: fpd, fvd, mpd, mvd
    real(kind=realtype), dimension(3) :: cofsumfx, cofsumfy, cofsumfz
    real(kind=realtype), dimension(3) :: cofsumfxd, cofsumfyd, cofsumfzd
    real(kind=realtype) :: yplusmax, sepsensorks, sepsensor, &
&   sepsensoravg(3), sepsensorarea, cavitation, cpmin_ks_sum, &
&   sepsensorksarea
    real(kind=realtype) :: sepsensorksd, sepsensord, sepsensoravgd(3), &
&   sepsensoraread, cavitationd, cpmin_ks_sumd, sepsensorksaread
    integer(kind=inttype) :: i, j, ii, blk
    real(kind=realtype) :: pm1, fx, fy, fz, fn
    real(kind=realtype) :: pm1d, fxd, fyd, fzd
    real(kind=realtype) :: vecttangential(3)
    real(kind=realtype) :: vecttangentiald(3)
    real(kind=realtype) :: vectdotproductfsnormal
    real(kind=realtype) :: vectdotproductfsnormald
    real(kind=realtype) :: xc, xco, yc, yco, zc, zco, qf(3), r(3), n(3)&
&   , l
    real(kind=realtype) :: xcd, xcod, ycd, ycod, zcd, zcod, rd(3)
    real(kind=realtype) :: fact, rho, mul, yplus, dwall
    real(kind=realtype) :: v(3), sensor, sensor1, cp, tmp, plocal, &
&   ks_exponent
    real(kind=realtype) :: vd(3), sensord, sensor1d, cpd, tmpd, plocald&
&   , ks_exponentd
    real(kind=realtype) :: tauxx, tauyy, tauzz
    real(kind=realtype) :: tauxxd, tauyyd, tauzzd
    real(kind=realtype) :: tauxy, tauxz, tauyz
    real(kind=realtype) :: tauxyd, tauxzd, tauyzd
    real(kind=realtype), dimension(3) :: refpoint
    real(kind=realtype), dimension(3) :: refpointd
    real(kind=realtype), dimension(3, 2) :: axispoints
    real(kind=realtype) :: mx, my, mz, cellarea, m0x, m0y, m0z, mvaxis, &
&   mpaxis
    real(kind=realtype) :: mxd, myd, mzd, cellaread, m0xd, m0yd, m0zd, &
&   mvaxisd, mpaxisd
    real(kind=realtype) :: cperror, cperror2
    real(kind=realtype) :: cperrord, cperror2d
    intrinsic mod
    intrinsic max
    intrinsic sqrt
    intrinsic cos
    intrinsic exp
    real(kind=realtype) :: temp
    real(kind=realtype) :: tempd
    real(kind=realtype), dimension(3) :: tmp0
    real(kind=realtype), dimension(3) :: tmpd0
    real(kind=realtype) :: temp0
    real(kind=realtype), dimension(3) :: tmp1
    real(kind=realtype), dimension(3) :: tmpd1
    real(kind=realtype) :: tempd0
    real(kind=realtype) :: temp1
    real(kind=realtype) :: tempd1
    integer :: branch
    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.
    call pushreal8array(n, 3)
    call pushreal8array(r, 3)
    call pushreal8array(v, 3)
    call pushreal8array(vecttangential, 3)
!$fwd-of ii-loop 
!
!         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
    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)
      if (bcdata(mm)%iblank(i, j) .lt. 0) then
        blk = 0
      else
        blk = bcdata(mm)%iblank(i, j)
      end if
      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.
!
    if (bctype(mm) .eq. nswalladiabatic .or. bctype(mm) .eq. &
&       nswallisothermal) then
      call pushcontrol1b(0)
    else
      call pushcontrol1b(1)
    end if
    cperror2d = localvaluesd(icperror2)
    mpaxisd = localvaluesd(iaxismoment)
    mvaxisd = localvaluesd(iaxismoment)
    sepsensoravgd = 0.0_8
    sepsensoravgd = localvaluesd(isepavg:isepavg+2)
    cpmin_ks_sumd = localvaluesd(icpmin)
    cavitationd = localvaluesd(icavitation)
    sepsensoraread = localvaluesd(isepsensorarea)
    sepsensorksaread = localvaluesd(isepsensorksarea)
    sepsensorksd = localvaluesd(isepsensorks)
    sepsensord = localvaluesd(isepsensor)
    cofsumfzd = 0.0_8
    cofsumfzd = localvaluesd(icoforcez:icoforcez+2)
    cofsumfyd = 0.0_8
    cofsumfyd = localvaluesd(icoforcey:icoforcey+2)
    cofsumfxd = 0.0_8
    cofsumfxd = localvaluesd(icoforcex:icoforcex+2)
    mvd = 0.0_8
    mvd = localvaluesd(imv:imv+2)
    mpd = 0.0_8
    mpd = localvaluesd(imp:imp+2)
    fvd = 0.0_8
    fvd = localvaluesd(ifv:ifv+2)
    fpd = 0.0_8
    fpd = localvaluesd(ifp:ifp+2)
    call popcontrol1b(branch)
    if (branch .eq. 0) then
      rd = 0.0_8
      refpointd = 0.0_8
!$bwd-of 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 (bcdata(mm)%iblank(i, j) .lt. 0) then
          blk = 0
        else
          blk = bcdata(mm)%iblank(i, j)
        end if
        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.
! 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
! force-y
! force-z
! 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
! save the face based forces for the slice operations
! 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).
! compute the local value of y+. due to the usage
! of pointers there is on offset of -1 in dd2wall..
        mxd = blk*mvd(1)
        myd = blk*mvd(2)
        mzd = blk*mvd(3)
        tempd1 = blk*mvaxisd
        m0xd = n(1)*tempd1
        m0yd = n(2)*tempd1
        m0zd = n(3)*tempd1
        fzd = bcdatad(mm)%fv(i, j, 3) + r(2)*m0xd - r(1)*m0yd + zco*blk*&
&         cofsumfzd(3) + yco*blk*cofsumfzd(2) + xco*blk*cofsumfzd(1) + &
&         yc*mxd - xc*myd + blk*fvd(3)
        bcdatad(mm)%fv(i, j, 3) = 0.0_8
        fyd = bcdatad(mm)%fv(i, j, 2) + r(1)*m0zd + zco*blk*cofsumfyd(3)&
&         - r(3)*m0xd + yco*blk*cofsumfyd(2) + xco*blk*cofsumfyd(1) + xc&
&         *mzd + blk*fvd(2) - zc*mxd
        bcdatad(mm)%fv(i, j, 2) = 0.0_8
        fxd = bcdatad(mm)%fv(i, j, 1) + r(3)*m0yd - r(2)*m0zd + zco*blk*&
&         cofsumfxd(3) + yco*blk*cofsumfxd(2) + xco*blk*cofsumfxd(1) + &
&         zc*myd - yc*mzd + blk*fvd(1)
        bcdatad(mm)%fv(i, j, 1) = 0.0_8
        rd(1) = rd(1) + fy*m0zd - fz*m0yd
        rd(2) = rd(2) + fz*m0xd - fx*m0zd
        rd(3) = rd(3) + fx*m0yd - fy*m0xd
        tempd1 = fourth*rd(3)
        rd(3) = 0.0_8
        xxd(i, j, 3) = xxd(i, j, 3) + tempd1
        xxd(i+1, j, 3) = xxd(i+1, j, 3) + tempd1
        xxd(i, j+1, 3) = xxd(i, j+1, 3) + tempd1
        xxd(i+1, j+1, 3) = xxd(i+1, j+1, 3) + tempd1
        tempd1 = fourth*rd(2)
        rd(2) = 0.0_8
        xxd(i, j, 2) = xxd(i, j, 2) + tempd1
        xxd(i+1, j, 2) = xxd(i+1, j, 2) + tempd1
        xxd(i, j+1, 2) = xxd(i, j+1, 2) + tempd1
        xxd(i+1, j+1, 2) = xxd(i+1, j+1, 2) + tempd1
        tempd1 = fourth*rd(1)
        rd(1) = 0.0_8
        xxd(i, j, 1) = xxd(i, j, 1) + tempd1
        xxd(i+1, j, 1) = xxd(i+1, j, 1) + tempd1
        xxd(i, j+1, 1) = xxd(i, j+1, 1) + tempd1
        xxd(i+1, j+1, 1) = xxd(i+1, j+1, 1) + tempd1
        zcod = fz*blk*cofsumfzd(3) + fy*blk*cofsumfyd(3) + fx*blk*&
&         cofsumfxd(3)
        ycod = fz*blk*cofsumfzd(2) + fy*blk*cofsumfyd(2) + fx*blk*&
&         cofsumfxd(2)
        xcod = fz*blk*cofsumfzd(1) + fy*blk*cofsumfyd(1) + fx*blk*&
&         cofsumfxd(1)
        tempd1 = fourth*zcod
        xxd(i, j, 3) = xxd(i, j, 3) + tempd1
        xxd(i+1, j, 3) = xxd(i+1, j, 3) + tempd1
        xxd(i, j+1, 3) = xxd(i, j+1, 3) + tempd1
        xxd(i+1, j+1, 3) = xxd(i+1, j+1, 3) + tempd1
        tempd1 = fourth*ycod
        xxd(i, j, 2) = xxd(i, j, 2) + tempd1
        xxd(i+1, j, 2) = xxd(i+1, j, 2) + tempd1
        xxd(i, j+1, 2) = xxd(i, j+1, 2) + tempd1
        xxd(i+1, j+1, 2) = xxd(i+1, j+1, 2) + tempd1
        tempd1 = fourth*xcod
        xxd(i, j, 1) = xxd(i, j, 1) + tempd1
        xxd(i+1, j, 1) = xxd(i+1, j, 1) + tempd1
        xxd(i, j+1, 1) = xxd(i, j+1, 1) + tempd1
        xxd(i+1, j+1, 1) = xxd(i+1, j+1, 1) + tempd1
        xcd = fy*mzd - fz*myd
        ycd = fz*mxd - fx*mzd
        zcd = fx*myd - fy*mxd
        tempd1 = fourth*zcd
        refpointd(3) = refpointd(3) - zcd
        xxd(i, j, 3) = xxd(i, j, 3) + tempd1
        xxd(i+1, j, 3) = xxd(i+1, j, 3) + tempd1
        xxd(i, j+1, 3) = xxd(i, j+1, 3) + tempd1
        xxd(i+1, j+1, 3) = xxd(i+1, j+1, 3) + tempd1
        tempd1 = fourth*ycd
        refpointd(2) = refpointd(2) - ycd
        xxd(i, j, 2) = xxd(i, j, 2) + tempd1
        xxd(i+1, j, 2) = xxd(i+1, j, 2) + tempd1
        xxd(i, j+1, 2) = xxd(i, j+1, 2) + tempd1
        xxd(i+1, j+1, 2) = xxd(i+1, j+1, 2) + tempd1
        tempd1 = fourth*xcd
        refpointd(1) = refpointd(1) - xcd
        xxd(i, j, 1) = xxd(i, j, 1) + tempd1
        xxd(i+1, j, 1) = xxd(i+1, j, 1) + tempd1
        xxd(i, j+1, 1) = xxd(i, j+1, 1) + tempd1
        xxd(i+1, j+1, 1) = xxd(i+1, j+1, 1) + tempd1
        tempd1 = -(pref*fact*fzd)
        prefd = prefd - (tauxz*ssi(i, j, 1)+tauyz*ssi(i, j, 2)+tauzz*ssi&
&         (i, j, 3))*fact*fzd - (tauxy*ssi(i, j, 1)+tauyy*ssi(i, j, 2)+&
&         tauyz*ssi(i, j, 3))*fact*fyd - (tauxx*ssi(i, j, 1)+tauxy*ssi(i&
&         , j, 2)+tauxz*ssi(i, j, 3))*fact*fxd
        tauxzd = ssi(i, j, 1)*tempd1
        ssid(i, j, 1) = ssid(i, j, 1) + tauxz*tempd1
        tauyzd = ssi(i, j, 2)*tempd1
        ssid(i, j, 2) = ssid(i, j, 2) + tauyz*tempd1
        tauzzd = ssi(i, j, 3)*tempd1
        ssid(i, j, 3) = ssid(i, j, 3) + tauzz*tempd1
        tempd1 = -(pref*fact*fyd)
        tauxyd = ssi(i, j, 1)*tempd1
        ssid(i, j, 1) = ssid(i, j, 1) + tauxy*tempd1
        tauyyd = ssi(i, j, 2)*tempd1
        ssid(i, j, 2) = ssid(i, j, 2) + tauyy*tempd1
        tauyzd = tauyzd + ssi(i, j, 3)*tempd1
        ssid(i, j, 3) = ssid(i, j, 3) + tauyz*tempd1
        tempd1 = -(pref*fact*fxd)
        tauxxd = ssi(i, j, 1)*tempd1
        ssid(i, j, 1) = ssid(i, j, 1) + tauxx*tempd1
        tauxyd = tauxyd + ssi(i, j, 2)*tempd1
        ssid(i, j, 2) = ssid(i, j, 2) + tauxy*tempd1
        tauxzd = tauxzd + ssi(i, j, 3)*tempd1
        ssid(i, j, 3) = ssid(i, j, 3) + tauxz*tempd1
        viscsubfaced(mm)%tau(i, j, 6) = viscsubfaced(mm)%tau(i, j, 6) + &
&         tauyzd
        viscsubfaced(mm)%tau(i, j, 5) = viscsubfaced(mm)%tau(i, j, 5) + &
&         tauxzd
        viscsubfaced(mm)%tau(i, j, 4) = viscsubfaced(mm)%tau(i, j, 4) + &
&         tauxyd
        viscsubfaced(mm)%tau(i, j, 3) = viscsubfaced(mm)%tau(i, j, 3) + &
&         tauzzd
        viscsubfaced(mm)%tau(i, j, 2) = viscsubfaced(mm)%tau(i, j, 2) + &
&         tauyyd
        viscsubfaced(mm)%tau(i, j, 1) = viscsubfaced(mm)%tau(i, j, 1) + &
&         tauxxd
      end do
    else
      bcdatad(mm)%fv = 0.0_8
      rd = 0.0_8
      refpointd = 0.0_8
    end if
    vd = 0.0_8
    vecttangentiald = 0.0_8
    call popreal8array(vecttangential, 3)
    call popreal8array(v, 3)
    call popreal8array(r, 3)
    call popreal8array(n, 3)
!$bwd-of 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)
      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)
      if (bcdata(mm)%iblank(i, j) .lt. 0) then
        blk = 0
      else
        blk = bcdata(mm)%iblank(i, j)
      end if
      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
! 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
! force-y
! force-z
! 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
! save the face-based forces and area
      cellarea = sqrt(ssi(i, j, 1)**2 + ssi(i, j, 2)**2 + ssi(i, j, 3)**&
&       2)
! get normalized surface velocity:
      v(1) = ww2(i, j, ivx)
      v(2) = ww2(i, j, ivy)
      v(3) = ww2(i, j, ivz)
      tmp0 = v/(sqrt(v(1)**2+v(2)**2+v(3)**2)+1e-16)
      call pushreal8array(v, 3)
      v = tmp0
      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)
        tmp1 = vecttangential/(sqrt(vecttangential(1)**2+vecttangential(&
&         2)**2+vecttangential(3)**2)+1e-16)
        call pushreal8array(vecttangential, 3)
        vecttangential = tmp1
! 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 pushcontrol1b(0)
      else
        call pushcontrol1b(1)
      end if
! dot product with free stream
      call pushreal8(sensor)
      sensor = -(v(1)*veldirfreestream(1)+v(2)*veldirfreestream(2)+v(3)*&
&       veldirfreestream(3))
!now run through a smooth heaviside function:
      call pushreal8(sensor)
      sensor = one/(one+exp(-(2*sepsensorsharpness*(sensor-&
&       sepsensoroffset))))
! and integrate over the area of this cell and save, blanking as we go.
      call pushreal8(sensor)
      sensor = sensor*cellarea*blk
! also accumulate into the sepsensoravg
! x-y-zco are computed above for center of force computations
      if (computecavitation) then
        plocal = pp2(i, j)
        tmp = two/(gammainf*machcoef*machcoef)
        cp = tmp*(plocal-pinf)
        sensor1 = -cp - cavitationnumber
        call pushreal8(sensor1)
        sensor1 = sensor1**cavexponent/(one+exp(2*cavsensorsharpness*(-&
&         sensor1+cavsensoroffset)))
! also do the ks-based cpmin computation
        sensor1d = cavitationd
        cellaread = sensor1*blk*sensor1d
        sensor1d = cellarea*blk*sensor1d
        call popreal8(sensor1)
        temp1 = 2*cavsensorsharpness*(cavsensoroffset-sensor1)
        temp0 = one + exp(temp1)
        if (sensor1 .le. 0.0_8 .and. (cavexponent .eq. 0.0_8 .or. &
&           cavexponent .ne. int(cavexponent))) then
          sensor1d = cavsensorsharpness*2*exp(temp1)*sensor1**&
&           cavexponent*sensor1d/temp0**2
        else
          sensor1d = (cavexponent*sensor1**(cavexponent-1)/temp0+&
&           cavsensorsharpness*2*exp(temp1)*sensor1**cavexponent/temp0**&
&           2)*sensor1d
        end if
        ks_exponentd = blk*cpmin_ks_sumd
        cpd = -(cpmin_rho*exp(cpmin_rho*(cpmin_family(spectralsol)-cp))*&
&         ks_exponentd) - sensor1d
        tmpd = (plocal-pinf)*cpd
        plocald = tmp*cpd
        pinfd = pinfd - tmp*cpd
        temp1 = gammainf*(machcoef*machcoef)
        machcoefd = machcoefd - 2*machcoef*gammainf*two*tmpd/temp1**2
        tmp = two/(gammainf*pinf*machcoef*machcoef)
        pp2d(i, j) = pp2d(i, j) + plocald
      else
        cellaread = 0.0_8
      end if
      sensord = zco*sepsensoravgd(3) + yco*sepsensoravgd(2) + xco*&
&       sepsensoravgd(1) + sepsensord
      zcod = sensor*sepsensoravgd(3)
      ycod = sensor*sepsensoravgd(2)
      xcod = sensor*sepsensoravgd(1)
      call popreal8(sensor)
      cellaread = cellaread + sensor*blk*sensord
      sensord = cellarea*blk*sensord
      call popreal8(sensor)
      temp1 = -(2*sepsensorsharpness*(sensor-sepsensoroffset))
      temp0 = one + exp(temp1)
      sensord = sepsensorsharpness*2*exp(temp1)*one*sensord/temp0**2
      call popreal8(sensor)
      vd(1) = vd(1) - veldirfreestream(1)*sensord
      veldirfreestreamd(1) = veldirfreestreamd(1) - v(1)*sensord
      vd(2) = vd(2) - veldirfreestream(2)*sensord
      veldirfreestreamd(2) = veldirfreestreamd(2) - v(2)*sensord
      vd(3) = vd(3) - veldirfreestream(3)*sensord
      veldirfreestreamd(3) = veldirfreestreamd(3) - v(3)*sensord
      call popcontrol1b(branch)
      if (branch .eq. 0) then
        ks_exponentd = blk*sepsensorksd
        temp0 = -(2*sepsensorkssharpness*(sepsensorksoffset+sensor))
        temp = one + exp(temp0)
        tempd1 = blk*one*sepsensoraread/temp
        cellaread = cellaread + tempd1 + blk*sepsensorksaread
        sensord = sepsensorkssharpness*2*exp(temp0)*cellarea*tempd1/temp
        call ksaggregationfunction_b(sensor, sensord, sepsenmaxfamily(&
&                              spectralsol), sepsenmaxrho, ks_exponent, &
&                              ks_exponentd)
        sensord = -(sensord/(cos(zero)-cos(degtorad*sepsensorksphi)+&
&         1e-16))
        vd(1) = vd(1) + vecttangential(1)*sensord
        vecttangentiald(1) = vecttangentiald(1) + v(1)*sensord
        vd(2) = vd(2) + vecttangential(2)*sensord
        vecttangentiald(2) = vecttangentiald(2) + v(2)*sensord
        vd(3) = vd(3) + vecttangential(3)*sensord
        vecttangentiald(3) = vecttangentiald(3) + v(3)*sensord
        call popreal8array(vecttangential, 3)
        tmpd1 = vecttangentiald
        temp0 = vecttangential(1)*vecttangential(1) + vecttangential(2)*&
&         vecttangential(2) + vecttangential(3)*vecttangential(3)
        temp = sqrt(temp0)
        vecttangentiald = tmpd1/(temp+1e-16)
        if (temp0 .eq. 0.0_8) then
          tempd0 = 0.0_8
        else
          tempd0 = -(sum(vecttangential*tmpd1)/(2.0*temp*(temp+1e-16)**2&
&           ))
        end if
        vecttangentiald(1) = vecttangentiald(1) + 2*vecttangential(1)*&
&         tempd0
        vecttangentiald(2) = vecttangentiald(2) + 2*vecttangential(2)*&
&         tempd0
        vecttangentiald(3) = vecttangentiald(3) + 2*vecttangential(3)*&
&         tempd0
        vectdotproductfsnormald = -(bcdata(mm)%norm(i, j, 3)*&
&         vecttangentiald(3)) - bcdata(mm)%norm(i, j, 2)*vecttangentiald&
&         (2) - bcdata(mm)%norm(i, j, 1)*vecttangentiald(1)
        veldirfreestreamd(3) = veldirfreestreamd(3) + vecttangentiald(3)&
&         + bcdata(mm)%norm(i, j, 3)*vectdotproductfsnormald
        vecttangentiald(3) = 0.0_8
        veldirfreestreamd(2) = veldirfreestreamd(2) + vecttangentiald(2)&
&         + bcdata(mm)%norm(i, j, 2)*vectdotproductfsnormald
        vecttangentiald(2) = 0.0_8
        veldirfreestreamd(1) = veldirfreestreamd(1) + vecttangentiald(1)&
&         + bcdata(mm)%norm(i, j, 1)*vectdotproductfsnormald
        vecttangentiald(1) = 0.0_8
      end if
      mxd = blk*mpd(1)
      myd = blk*mpd(2)
      mzd = blk*mpd(3)
      call popreal8array(v, 3)
      tmpd0 = vd
      temp = v(1)*v(1) + v(2)*v(2) + v(3)*v(3)
      temp0 = sqrt(temp)
      vd = tmpd0/(temp0+1e-16)
      if (temp .eq. 0.0_8) then
        tempd = 0.0_8
      else
        tempd = -(sum(v*tmpd0)/(2.0*temp0*(temp0+1e-16)**2))
      end if
      vd(1) = vd(1) + 2*v(1)*tempd
      vd(2) = vd(2) + 2*v(2)*tempd
      vd(3) = vd(3) + 2*v(3)*tempd
      ww2d(i, j, ivz) = ww2d(i, j, ivz) + vd(3)
      vd(3) = 0.0_8
      ww2d(i, j, ivy) = ww2d(i, j, ivy) + vd(2)
      vd(2) = 0.0_8
      ww2d(i, j, ivx) = ww2d(i, j, ivx) + vd(1)
      vd(1) = 0.0_8
      cellaread = cellaread + bcdatad(mm)%area(i, j)
      bcdatad(mm)%area(i, j) = 0.0_8
      if (ssi(i, j, 1)**2 + ssi(i, j, 2)**2 + ssi(i, j, 3)**2 .eq. 0.0_8&
&     ) then
        tempd = 0.0_8
      else
        tempd = cellaread/(2.0*sqrt(ssi(i, j, 1)**2+ssi(i, j, 2)**2+ssi(&
&         i, j, 3)**2))
      end if
      ssid(i, j, 1) = ssid(i, j, 1) + 2*ssi(i, j, 1)*tempd
      ssid(i, j, 2) = ssid(i, j, 2) + 2*ssi(i, j, 2)*tempd
      ssid(i, j, 3) = ssid(i, j, 3) + 2*ssi(i, j, 3)*tempd
      tempd = blk*mpaxisd
      m0xd = n(1)*tempd
      m0yd = n(2)*tempd
      fzd = bcdatad(mm)%fp(i, j, 3) + r(2)*m0xd - r(1)*m0yd + zco*blk*&
&       cofsumfzd(3) + yco*blk*cofsumfzd(2) + xco*blk*cofsumfzd(1) + yc*&
&       mxd - xc*myd + blk*fpd(3)
      bcdatad(mm)%fp(i, j, 3) = 0.0_8
      m0zd = n(3)*tempd
      fyd = bcdatad(mm)%fp(i, j, 2) + r(1)*m0zd + zco*blk*cofsumfyd(3) -&
&       r(3)*m0xd + yco*blk*cofsumfyd(2) + xco*blk*cofsumfyd(1) + xc*mzd&
&       + blk*fpd(2) - zc*mxd
      bcdatad(mm)%fp(i, j, 2) = 0.0_8
      fxd = bcdatad(mm)%fp(i, j, 1) + r(3)*m0yd - r(2)*m0zd + zco*blk*&
&       cofsumfxd(3) + yco*blk*cofsumfxd(2) + xco*blk*cofsumfxd(1) + zc*&
&       myd - yc*mzd + blk*fpd(1)
      bcdatad(mm)%fp(i, j, 1) = 0.0_8
      rd(1) = rd(1) + fy*m0zd - fz*m0yd
      rd(2) = rd(2) + fz*m0xd - fx*m0zd
      rd(3) = rd(3) + fx*m0yd - fy*m0xd
      tempd = fourth*rd(3)
      rd(3) = 0.0_8
      xxd(i, j, 3) = xxd(i, j, 3) + tempd
      xxd(i+1, j, 3) = xxd(i+1, j, 3) + tempd
      xxd(i, j+1, 3) = xxd(i, j+1, 3) + tempd
      xxd(i+1, j+1, 3) = xxd(i+1, j+1, 3) + tempd
      tempd = fourth*rd(2)
      rd(2) = 0.0_8
      xxd(i, j, 2) = xxd(i, j, 2) + tempd
      xxd(i+1, j, 2) = xxd(i+1, j, 2) + tempd
      xxd(i, j+1, 2) = xxd(i, j+1, 2) + tempd
      xxd(i+1, j+1, 2) = xxd(i+1, j+1, 2) + tempd
      tempd = fourth*rd(1)
      rd(1) = 0.0_8
      xxd(i, j, 1) = xxd(i, j, 1) + tempd
      xxd(i+1, j, 1) = xxd(i+1, j, 1) + tempd
      xxd(i, j+1, 1) = xxd(i, j+1, 1) + tempd
      xxd(i+1, j+1, 1) = xxd(i+1, j+1, 1) + tempd
      zcod = zcod + fz*blk*cofsumfzd(3) + fy*blk*cofsumfyd(3) + fx*blk*&
&       cofsumfxd(3)
      ycod = ycod + fz*blk*cofsumfzd(2) + fy*blk*cofsumfyd(2) + fx*blk*&
&       cofsumfxd(2)
      xcod = xcod + fz*blk*cofsumfzd(1) + fy*blk*cofsumfyd(1) + fx*blk*&
&       cofsumfxd(1)
      tempd = fourth*zcod
      xxd(i, j, 3) = xxd(i, j, 3) + tempd
      xxd(i+1, j, 3) = xxd(i+1, j, 3) + tempd
      xxd(i, j+1, 3) = xxd(i, j+1, 3) + tempd
      xxd(i+1, j+1, 3) = xxd(i+1, j+1, 3) + tempd
      tempd = fourth*ycod
      xxd(i, j, 2) = xxd(i, j, 2) + tempd
      xxd(i+1, j, 2) = xxd(i+1, j, 2) + tempd
      xxd(i, j+1, 2) = xxd(i, j+1, 2) + tempd
      xxd(i+1, j+1, 2) = xxd(i+1, j+1, 2) + tempd
      tempd = fourth*xcod
      xxd(i, j, 1) = xxd(i, j, 1) + tempd
      xxd(i+1, j, 1) = xxd(i+1, j, 1) + tempd
      xxd(i, j+1, 1) = xxd(i, j+1, 1) + tempd
      xxd(i+1, j+1, 1) = xxd(i+1, j+1, 1) + tempd
      xcd = fy*mzd - fz*myd
      ycd = fz*mxd - fx*mzd
      zcd = fx*myd - fy*mxd
      pm1d = ssi(i, j, 3)*fzd + ssi(i, j, 2)*fyd + ssi(i, j, 1)*fxd
      ssid(i, j, 3) = ssid(i, j, 3) + pm1*fzd
      ssid(i, j, 2) = ssid(i, j, 2) + pm1*fyd
      ssid(i, j, 1) = ssid(i, j, 1) + pm1*fxd
      tempd = fourth*zcd
      refpointd(3) = refpointd(3) - zcd
      xxd(i, j, 3) = xxd(i, j, 3) + tempd
      xxd(i+1, j, 3) = xxd(i+1, j, 3) + tempd
      xxd(i, j+1, 3) = xxd(i, j+1, 3) + tempd
      xxd(i+1, j+1, 3) = xxd(i+1, j+1, 3) + tempd
      tempd = fourth*ycd
      refpointd(2) = refpointd(2) - ycd
      xxd(i, j, 2) = xxd(i, j, 2) + tempd
      xxd(i+1, j, 2) = xxd(i+1, j, 2) + tempd
      xxd(i, j+1, 2) = xxd(i, j+1, 2) + tempd
      xxd(i+1, j+1, 2) = xxd(i+1, j+1, 2) + tempd
      tempd = fourth*xcd
      refpointd(1) = refpointd(1) - xcd
      xxd(i, j, 1) = xxd(i, j, 1) + tempd
      xxd(i+1, j, 1) = xxd(i+1, j, 1) + tempd
      xxd(i, j+1, 1) = xxd(i, j+1, 1) + tempd
      xxd(i+1, j+1, 1) = xxd(i+1, j+1, 1) + tempd
      cperrord = 2*cperror*cperror2d
      cpd = cperrord
      tempd = tmp*cpd
      tmpd = (half*(pp2(i, j)+pp1(i, j))-pinf)*cpd
      pp2d(i, j) = pp2d(i, j) + half*tempd
      pp1d(i, j) = pp1d(i, j) + half*tempd
      pinfd = pinfd - tempd
      temp = gammainf*pinf*(machcoef*machcoef)
      tempd = -(two*tmpd/temp**2)
      pinfd = pinfd + gammainf*machcoef**2*tempd
      machcoefd = machcoefd + 2*machcoef*gammainf*pinf*tempd
      tempd = pref*fact*pm1d
      prefd = prefd + (half*(pp2(i, j)+pp1(i, j))-pinf)*fact*pm1d
      pp2d(i, j) = pp2d(i, j) + half*tempd
      pp1d(i, j) = pp1d(i, j) + half*tempd
      pinfd = pinfd - tempd
    end do
    pointrefd(3) = pointrefd(3) + lref*refpointd(3)
    refpointd(3) = 0.0_8
    pointrefd(2) = pointrefd(2) + lref*refpointd(2)
    refpointd(2) = 0.0_8
    pointrefd(1) = pointrefd(1) + lref*refpointd(1)
  end subroutine wallintegrationface_b
 
  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_b
    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
    intrinsic mod
    intrinsic max
    intrinsic sqrt
    intrinsic cos
    intrinsic exp
    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
!$ad ii-loop
!
!         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
    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)
      if (bcdata(mm)%iblank(i, j) .lt. 0) then
        blk = 0
      else
        blk = bcdata(mm)%iblank(i, j)
      end if
      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.
!
    if (bctype(mm) .eq. nswalladiabatic .or. bctype(mm) .eq. &
&       nswallisothermal) then
! initialize dwall for the laminar case and set the pointer
! for the unit normals.
      dwall = zero
!$ad ii-loop
! 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+.
      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 (bcdata(mm)%iblank(i, j) .lt. 0) then
          blk = 0
        else
          blk = bcdata(mm)%iblank(i, j)
        end if
        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..
      end do
    else
! if we had no viscous force, set the viscous component to zero
      bcdata(mm)%fv = zero
    end if
! 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
  end subroutine wallintegrationface
 
!  differentiation of ksaggregationfunction in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: ks_g g
!   with respect to varying inputs: g
  subroutine ksaggregationfunction_b(g, gd, max_g, g_rho, ks_g, ks_gd)
    use precision
    implicit none
    real(kind=realtype), intent(inout) :: ks_g
    real(kind=realtype), intent(inout) :: ks_gd
    real(kind=realtype), intent(in) :: g, max_g, g_rho
    real(kind=realtype) :: gd
    intrinsic exp
    gd = gd + g_rho*exp(g_rho*(g-max_g))*ks_gd
  end subroutine ksaggregationfunction_b
 
  subroutine ksaggregationfunction(g, max_g, g_rho, ks_g)
    use precision
    implicit none
    real(kind=realtype), intent(inout) :: ks_g
    real(kind=realtype), intent(in) :: g, max_g, g_rho
    intrinsic exp
    ks_g = exp(g_rho*(g-max_g))
  end subroutine ksaggregationfunction
 
!  differentiation of flowintegrationface in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: timeref tref pref rgas rhoref
!                *xx *pp1 *pp2 *ssi *ww1 *ww2 pointref localvalues
!   with respect to varying inputs: timeref tref pref rgas rhoref
!                *xx *pp1 *pp2 *ssi *ww1 *ww2 pointref localvalues
!   rw status of diff variables: timeref:incr tref:incr pref:incr
!                rgas:incr rhoref:incr *xx:incr *pp1:incr *pp2:incr
!                *ssi:incr *ww1:incr *ww2:incr pointref:incr localvalues:in-out
!   plus diff mem management of: xx:in pp1:in pp2:in ssi:in ww1:in
!                ww2:in
  subroutine flowintegrationface_b(isinflow, localvalues, localvaluesd, &
&   mm)
    use constants
    use blockpointers, only : bctype, bcfaceid, bcdata, bcdatad, &
&   addgridvelocities
    use flowvarrefstate, only : pref, prefd, pinf, pinfd, rhoref, &
&   rhorefd, timeref, timerefd, lref, tref, trefd, rgas, rgasd, uref, &
&   urefd, uinf, uinfd, rhoinf, rhoinfd, gammainf
    use inputphysics, only : pointref, pointrefd, flowtype, rgasdim
    use flowutils_b, only : computeptot, computeptot_b, computettot, &
&   computettot_b
    use bcpointers_b, only : ssi, ssid, sface, ww1, ww1d, ww2, ww2d, pp1&
&   , pp1d, pp2, pp2d, xx, xxd, gamma1, gamma2
    use utils_b, only : mynorm2
    implicit none
! input/output variables
    logical, intent(in) :: isinflow
    real(kind=realtype), dimension(nlocalvalues), intent(inout) :: &
&   localvalues
    real(kind=realtype), dimension(nlocalvalues), intent(inout) :: &
&   localvaluesd
    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) :: massflowrated, mass_ptotd, mass_ttotd, &
&   mass_psd, mass_mnd, mass_ad, mass_rhod, mass_vxd, mass_vyd, mass_vzd&
&   , mass_nxd, mass_nyd, mass_nzd, mass_vid
    real(kind=realtype) :: area_ptot, area_ps
    real(kind=realtype) :: area_ptotd, area_psd
    real(kind=realtype) :: govgm1, gm1ovg, viconst, vilocal, pratio
    real(kind=realtype) :: vilocald, pratiod
    real(kind=realtype) :: mredim
    real(kind=realtype) :: mredimd
    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) :: xcd, xcod, ycd, ycod, zcd, zcod, mxd, myd, &
&   mzd
    real(kind=realtype) :: sf, vmag, vnm, vnmfreestreamref, vxm, vym, &
&   vzm, fx, fy, fz, u, v, w
    real(kind=realtype) :: vmagd, vnmd, vxmd, vymd, vzmd, fxd, fyd, fzd
    real(kind=realtype) :: pm, ptot, ttot, rhom, gammam, am
    real(kind=realtype) :: pmd, ptotd, ttotd, rhomd, amd
    real(kind=realtype) :: area, cellarea, overcellarea
    real(kind=realtype) :: aread, cellaread, overcellaread
    real(kind=realtype), dimension(3) :: fp, mp, fmom, mmom, refpoint, &
&   sfacecoordref
    real(kind=realtype), dimension(3) :: fpd, mpd, fmomd, mmomd, &
&   refpointd, sfacecoordrefd
    real(kind=realtype), dimension(3) :: cofsumfx, cofsumfy, cofsumfz
    real(kind=realtype), dimension(3) :: cofsumfxd, cofsumfyd, cofsumfzd
    real(kind=realtype) :: mnm, massflowratelocal
    real(kind=realtype) :: mnmd, massflowratelocald
    intrinsic sqrt
    intrinsic mod
    intrinsic max
    intrinsic min
    real(kind=realtype) :: temp
    real(kind=realtype) :: tempd
    real(kind=realtype) :: tempd0
    integer :: branch
    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 .eq. internalflow) internalflowfact = -one
    inflowfact = one
    if (isinflow) inflowfact = -one
! 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)
    mass_vid = localvaluesd(imassvi)
    mass_nzd = localvaluesd(imassnz)
    mass_nyd = localvaluesd(imassny)
    mass_nxd = localvaluesd(imassnx)
    mass_vzd = localvaluesd(imassvz)
    mass_vyd = localvaluesd(imassvy)
    mass_vxd = localvaluesd(imassvx)
    area_psd = localvaluesd(iareaps)
    area_ptotd = localvaluesd(iareaptot)
    cofsumfzd = 0.0_8
    cofsumfzd = localvaluesd(icoforcez:icoforcez+2)
    cofsumfyd = 0.0_8
    cofsumfyd = localvaluesd(icoforcey:icoforcey+2)
    cofsumfxd = 0.0_8
    cofsumfxd = localvaluesd(icoforcex:icoforcex+2)
    mmomd = 0.0_8
    mmomd = localvaluesd(iflowmm:iflowmm+2)
    mpd = 0.0_8
    mpd = localvaluesd(iflowmp:iflowmp+2)
    fmomd = 0.0_8
    fmomd = localvaluesd(iflowfm:iflowfm+2)
    fpd = 0.0_8
    fpd = localvaluesd(ifp:ifp+2)
    mass_mnd = localvaluesd(imassmn)
    mass_psd = localvaluesd(imassps)
    mass_ttotd = localvaluesd(imassttot)
    mass_ptotd = localvaluesd(imassptot)
    mass_ad = localvaluesd(imassa)
    mass_rhod = localvaluesd(imassrho)
    aread = localvaluesd(iarea)
    massflowrated = localvaluesd(imassflow)
    mredimd = 0.0_8
    ptotd = 0.0_8
    refpointd = 0.0_8
    sfacecoordrefd = 0.0_8
    ttotd = 0.0_8
!$bwd-of 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
      if (bcdata(mm)%iblank(i, j) .lt. 0) then
        blk = 0
      else
        blk = bcdata(mm)%iblank(i, j)
      end if
      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)
      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
! re-dimentionalize quantities
      call pushreal8(pm)
      pm = pm*pref
      sfacecoordref(1) = sf*ssi(i, j, 1)*overcellarea
      sfacecoordref(2) = sf*ssi(i, j, 2)*overcellarea
      sfacecoordref(3) = sf*ssi(i, j, 3)*overcellarea
      govgm1 = gammainf/(gammainf-one)
      gm1ovg = one/govgm1
      viconst = two*govgm1*rgasdim
      if (one .gt. one/ptot) then
        pratio = one/ptot
        call pushcontrol1b(0)
      else
        call pushcontrol1b(1)
        pratio = one
      end if
      vilocal = sqrt(viconst*(one-pratio**gm1ovg)*ttot*tref)
      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.
      call pushreal8(pm)
      pm = -((pm-pinf*pref)*fact*blk)
! update the pressure force and moment coefficients.
! 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
! force-y
! force-z
! 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
      call pushreal8(massflowratelocal)
      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
! 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
! force-y
! force-z
      mxd = mmomd(1)
      myd = mmomd(2)
      mzd = mmomd(3)
      zcod = fz*cofsumfzd(3) + fy*cofsumfyd(3) + fx*cofsumfxd(3)
      fzd = zco*cofsumfzd(3) + yco*cofsumfzd(2) + xco*cofsumfzd(1) + yc*&
&       mxd - xc*myd + fmomd(3)
      ycod = fz*cofsumfzd(2) + fy*cofsumfyd(2) + fx*cofsumfxd(2)
      xcod = fz*cofsumfzd(1) + fy*cofsumfyd(1) + fx*cofsumfxd(1)
      fyd = zco*cofsumfyd(3) + yco*cofsumfyd(2) + xco*cofsumfyd(1) + xc*&
&       mzd + fmomd(2) - zc*mxd
      fxd = zco*cofsumfxd(3) + yco*cofsumfxd(2) + xco*cofsumfxd(1) + zc*&
&       myd - yc*mzd + fmomd(1)
      xcd = fy*mzd - fz*myd
      ycd = fz*mxd - fx*mzd
      zcd = fx*myd - fy*mxd
      ssid(i, j, 3) = ssid(i, j, 3) + massflowratelocal*vzm*fzd
      tempd0 = ssi(i, j, 3)*fzd
      massflowratelocald = vzm*tempd0
      vzmd = massflowratelocal*tempd0
      ssid(i, j, 2) = ssid(i, j, 2) + massflowratelocal*vym*fyd
      tempd0 = ssi(i, j, 2)*fyd
      massflowratelocald = massflowratelocald + vym*tempd0
      vymd = massflowratelocal*tempd0
      ssid(i, j, 1) = ssid(i, j, 1) + massflowratelocal*vxm*fxd
      tempd0 = ssi(i, j, 1)*fxd
      massflowratelocald = massflowratelocald + vxm*tempd0
      vxmd = massflowratelocal*tempd0
      call popreal8(massflowratelocal)
      tempd0 = fact*blk*internalflowfact*inflowfact*massflowratelocald/(&
&       timeref*cellarea)
      massflowratelocald = tempd0
      tempd = -(massflowratelocal*tempd0/(timeref*cellarea))
      timerefd = timerefd + cellarea*tempd
      cellaread = timeref*tempd
      fz = pm*ssi(i, j, 3)
      fy = pm*ssi(i, j, 2)
      fx = pm*ssi(i, j, 1)
      zcod = zcod + fz*cofsumfzd(3) + fy*cofsumfyd(3) + fx*cofsumfxd(3)
      ycod = ycod + fz*cofsumfzd(2) + fy*cofsumfyd(2) + fx*cofsumfxd(2)
      xcod = xcod + fz*cofsumfzd(1) + fy*cofsumfyd(1) + fx*cofsumfxd(1)
      tempd0 = fourth*zcod
      xxd(i, j, 3) = xxd(i, j, 3) + tempd0
      xxd(i+1, j, 3) = xxd(i+1, j, 3) + tempd0
      xxd(i, j+1, 3) = xxd(i, j+1, 3) + tempd0
      xxd(i+1, j+1, 3) = xxd(i+1, j+1, 3) + tempd0
      tempd0 = fourth*ycod
      xxd(i, j, 2) = xxd(i, j, 2) + tempd0
      xxd(i+1, j, 2) = xxd(i+1, j, 2) + tempd0
      xxd(i, j+1, 2) = xxd(i, j+1, 2) + tempd0
      xxd(i+1, j+1, 2) = xxd(i+1, j+1, 2) + tempd0
      tempd0 = fourth*xcod
      xxd(i, j, 1) = xxd(i, j, 1) + tempd0
      xxd(i+1, j, 1) = xxd(i+1, j, 1) + tempd0
      xxd(i, j+1, 1) = xxd(i, j+1, 1) + tempd0
      xxd(i+1, j+1, 1) = xxd(i+1, j+1, 1) + tempd0
      mzd = mpd(3)
      myd = mpd(2)
      fxd = zco*cofsumfxd(3) + yco*cofsumfxd(2) + xco*cofsumfxd(1) + zc*&
&       myd - yc*mzd + fpd(1)
      mxd = mpd(1)
      fzd = zco*cofsumfzd(3) + yco*cofsumfzd(2) + xco*cofsumfzd(1) + yc*&
&       mxd - xc*myd + fpd(3)
      fyd = zco*cofsumfyd(3) + yco*cofsumfyd(2) + xco*cofsumfyd(1) + xc*&
&       mzd + fpd(2) - zc*mxd
      xcd = xcd + fy*mzd - fz*myd
      ycd = ycd + fz*mxd - fx*mzd
      zcd = zcd + fx*myd - fy*mxd
      pmd = ssi(i, j, 3)*fzd + ssi(i, j, 2)*fyd + ssi(i, j, 1)*fxd
      ssid(i, j, 3) = ssid(i, j, 3) + pm*fzd
      ssid(i, j, 2) = ssid(i, j, 2) + pm*fyd
      ssid(i, j, 1) = ssid(i, j, 1) + pm*fxd
      call popreal8(pm)
      tempd0 = -(fact*blk*pmd)
      pmd = tempd0
      prefd = prefd - pinf*tempd0
      tempd0 = fourth*zcd
      refpointd(3) = refpointd(3) - zcd
      xxd(i, j, 3) = xxd(i, j, 3) + tempd0
      xxd(i+1, j, 3) = xxd(i+1, j, 3) + tempd0
      xxd(i, j+1, 3) = xxd(i, j+1, 3) + tempd0
      xxd(i+1, j+1, 3) = xxd(i+1, j+1, 3) + tempd0
      tempd0 = fourth*ycd
      refpointd(2) = refpointd(2) - ycd
      xxd(i, j, 2) = xxd(i, j, 2) + tempd0
      xxd(i+1, j, 2) = xxd(i+1, j, 2) + tempd0
      xxd(i, j+1, 2) = xxd(i, j+1, 2) + tempd0
      xxd(i+1, j+1, 2) = xxd(i+1, j+1, 2) + tempd0
      tempd0 = fourth*xcd
      refpointd(1) = refpointd(1) - xcd
      xxd(i, j, 1) = xxd(i, j, 1) + tempd0
      xxd(i+1, j, 1) = xxd(i+1, j, 1) + tempd0
      xxd(i, j+1, 1) = xxd(i, j+1, 1) + tempd0
      xxd(i+1, j+1, 1) = xxd(i+1, j+1, 1) + tempd0
      ssid(i, j, 3) = ssid(i, j, 3) + overcellarea*massflowratelocal*&
&       mass_nzd
      tempd0 = ssi(i, j, 3)*mass_nzd
      overcellaread = massflowratelocal*tempd0
      massflowratelocald = massflowratelocald + overcellarea*tempd0
      ssid(i, j, 2) = ssid(i, j, 2) + overcellarea*massflowratelocal*&
&       mass_nyd
      tempd0 = ssi(i, j, 2)*mass_nyd
      overcellaread = overcellaread + massflowratelocal*tempd0
      massflowratelocald = massflowratelocald + overcellarea*tempd0
      ssid(i, j, 1) = ssid(i, j, 1) + overcellarea*massflowratelocal*&
&       mass_nxd
      tempd0 = ssi(i, j, 1)*mass_nxd
      overcellaread = overcellaread + massflowratelocal*tempd0
      massflowratelocald = massflowratelocald + overcellarea*tempd0 + &
&       vilocal*mass_vid
      vilocald = massflowratelocal*mass_vid
      temp = one - pratio**gm1ovg
      if (viconst*(temp*(ttot*tref)) .eq. 0.0_8) then
        tempd0 = 0.0_8
      else
        tempd0 = viconst*vilocald/(2.0*sqrt(viconst*(temp*(ttot*tref))))
      end if
      if (pratio .le. 0.0_8 .and. (gm1ovg .eq. 0.0_8 .or. gm1ovg .ne. &
&         int(gm1ovg))) then
        pratiod = 0.0_8
      else
        pratiod = -(gm1ovg*pratio**(gm1ovg-1)*ttot*tref*tempd0)
      end if
      ttotd = ttotd + tref*temp*tempd0
      trefd = trefd + ttot*temp*tempd0
      call popcontrol1b(branch)
      if (branch .eq. 0) ptotd = ptotd - one*pratiod/ptot**2
      tempd = blk*area_ptotd
      vzmd = vzmd + uref*massflowratelocal*mass_vzd
      sfacecoordrefd(3) = sfacecoordrefd(3) - massflowratelocal*mass_vzd
      massflowratelocald = massflowratelocald + (uref*vzm-sfacecoordref(&
&       3))*mass_vzd + (uref*vym-sfacecoordref(2))*mass_vyd + (uref*vxm-&
&       sfacecoordref(1))*mass_vxd + mnm*mass_mnd + pm*mass_psd + am*&
&       uref*mass_ad + rhom*rhoref*mass_rhod + ttot*tref*mass_ttotd + &
&       ptot*pref*mass_ptotd + massflowrated
      vymd = vymd + uref*massflowratelocal*mass_vyd
      sfacecoordrefd(2) = sfacecoordrefd(2) - massflowratelocal*mass_vyd
      vxmd = vxmd + uref*massflowratelocal*mass_vxd
      sfacecoordrefd(1) = sfacecoordrefd(1) - massflowratelocal*mass_vxd
      ssid(i, j, 3) = ssid(i, j, 3) + overcellarea*sf*sfacecoordrefd(3)
      overcellaread = overcellaread + ssi(i, j, 3)*sf*sfacecoordrefd(3) &
&       + ssi(i, j, 2)*sf*sfacecoordrefd(2) + ssi(i, j, 1)*sf*&
&       sfacecoordrefd(1)
      sfacecoordrefd(3) = 0.0_8
      ssid(i, j, 2) = ssid(i, j, 2) + overcellarea*sf*sfacecoordrefd(2)
      sfacecoordrefd(2) = 0.0_8
      ssid(i, j, 1) = ssid(i, j, 1) + overcellarea*sf*sfacecoordrefd(1)
      sfacecoordrefd(1) = 0.0_8
      pmd = pmd + cellarea*blk*area_psd + massflowratelocal*mass_psd
      cellaread = cellaread + pm*blk*area_psd + ptot*pref*tempd + blk*&
&       aread - overcellaread/cellarea**2
      ptotd = ptotd + pref*cellarea*tempd + massflowratelocal*pref*&
&       mass_ptotd
      mnmd = massflowratelocal*mass_mnd
      amd = massflowratelocal*uref*mass_ad - vmag*mnmd/am**2
      rhorefd = rhorefd + rhom*massflowratelocal*mass_rhod
      ttotd = ttotd + massflowratelocal*tref*mass_ttotd
      trefd = trefd + ttot*massflowratelocal*mass_ttotd
      call popreal8(pm)
      prefd = prefd + ptot*cellarea*tempd + ptot*massflowratelocal*&
&       mass_ptotd + pm*pmd
      pmd = pref*pmd
      tempd = blk*fact*massflowratelocald
      rhomd = massflowratelocal*rhoref*mass_rhod + vnm*mredim*tempd
      vnmd = rhom*mredim*tempd
      mredimd = mredimd + rhom*vnm*tempd
      call computettot_b(rhom, rhomd, vxm, vxmd, vym, vymd, vzm, vzmd, &
&                  pm, pmd, ttot, ttotd)
      call computeptot_b(rhom, rhomd, vxm, vxmd, vym, vymd, vzm, vzmd, &
&                  pm, pmd, ptot, ptotd)
      if (ssi(i, j, 1)**2 + ssi(i, j, 2)**2 + ssi(i, j, 3)**2 .eq. 0.0_8&
&     ) then
        tempd = 0.0_8
      else
        tempd = cellaread/(2.0*sqrt(ssi(i, j, 1)**2+ssi(i, j, 2)**2+ssi(&
&         i, j, 3)**2))
      end if
      ssid(i, j, 3) = ssid(i, j, 3) + 2*ssi(i, j, 3)*tempd + vzm*vnmd
      ssid(i, j, 2) = ssid(i, j, 2) + 2*ssi(i, j, 2)*tempd + vym*vnmd
      ssid(i, j, 1) = ssid(i, j, 1) + 2*ssi(i, j, 1)*tempd + vxm*vnmd
      vmagd = mnmd/am
      if (gammam*(pm/rhom) .eq. 0.0_8) then
        tempd = 0.0_8
      else
        tempd = gammam*amd/(rhom*2.0*sqrt(gammam*(pm/rhom)))
      end if
      pmd = pmd + tempd
      rhomd = rhomd - pm*tempd/rhom
      if (vxm**2 + vym**2 + vzm**2 .eq. 0.0_8) then
        tempd = 0.0_8
      else
        tempd = vmagd/(2.0*sqrt(vxm**2+vym**2+vzm**2))
      end if
      vxmd = vxmd + 2*vxm*tempd + ssi(i, j, 1)*vnmd
      vymd = vymd + 2*vym*tempd + ssi(i, j, 2)*vnmd
      vzmd = vzmd + 2*vzm*tempd + ssi(i, j, 3)*vnmd
      pp1d(i, j) = pp1d(i, j) + half*pmd
      pp2d(i, j) = pp2d(i, j) + half*pmd
      ww1d(i, j, irho) = ww1d(i, j, irho) + half*rhomd
      ww2d(i, j, irho) = ww2d(i, j, irho) + half*rhomd
      ww1d(i, j, ivz) = ww1d(i, j, ivz) + half*vzmd
      ww2d(i, j, ivz) = ww2d(i, j, ivz) + half*vzmd
      ww1d(i, j, ivy) = ww1d(i, j, ivy) + half*vymd
      ww2d(i, j, ivy) = ww2d(i, j, ivy) + half*vymd
      ww1d(i, j, ivx) = ww1d(i, j, ivx) + half*vxmd
      ww2d(i, j, ivx) = ww2d(i, j, ivx) + half*vxmd
    end do
    if (pref*rhoref .eq. 0.0_8) then
      tempd = 0.0_8
    else
      tempd = mredimd/(2.0*sqrt(pref*rhoref))
    end if
    prefd = prefd + rhoref*tempd
    rhorefd = rhorefd + pref*tempd
    pointrefd(3) = pointrefd(3) + lref*refpointd(3)
    refpointd(3) = 0.0_8
    pointrefd(2) = pointrefd(2) + lref*refpointd(2)
    refpointd(2) = 0.0_8
    pointrefd(1) = pointrefd(1) + lref*refpointd(1)
  end subroutine flowintegrationface_b
 
  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_b, only : computeptot, computettot
    use bcpointers_b, only : ssi, sface, ww1, ww2, pp1, pp2, xx, gamma1,&
&   gamma2
    use utils_b, 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
    intrinsic sqrt
    intrinsic mod
    intrinsic max
    intrinsic min
    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 .eq. internalflow) internalflowfact = -one
    inflowfact = one
    if (isinflow) inflowfact = -one
! 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
      if (bcdata(mm)%iblank(i, j) .lt. 0) then
        blk = 0
      else
        blk = bcdata(mm)%iblank(i, j)
      end if
      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
      if (one .gt. one/ptot) then
        pratio = one/ptot
      else
        pratio = one
      end if
      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               |
!                                                                      |
! ----------------------------------------------------------------------
 
end module surfaceintegrations_b