surfaceIntegrations_d.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_d
  implicit none
 
contains
!  differentiation of getcostfunctions in forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: funcvalues
!   with respect to varying inputs: pref machcoef dragdirection
!                liftdirection globalvals
!   rw status of diff variables: pref:in machcoef:in dragdirection:in
!                liftdirection:in globalvals:in funcvalues:out
  subroutine getcostfunctions_d(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_d, only : computetsderivatives
    use flowutils_d, only : getdirvector
    implicit none
! input/output
    real(kind=realtype), dimension(:, :), intent(in) :: globalvals
    real(kind=realtype), dimension(:, :), intent(in) :: globalvalsd
    real(kind=realtype), dimension(:), intent(out) :: funcvalues
    real(kind=realtype), dimension(:), intent(out) :: 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) :: arg1
    real(kind=realtype) :: arg1d
    real(kind=realtype) :: temp
    real(kind=realtype) :: temp0
    real(kind=realtype), dimension(3) :: temp1
! factor used for time-averaged quantities.
    ovrnts = one/ntimeintervalsspectral
! sum pressure and viscous contributions
    forced = globalvalsd(ifp:ifp+2, :) + globalvalsd(ifv:ifv+2, :) + &
&     globalvalsd(iflowfm:iflowfm+2, :)
    force = globalvals(ifp:ifp+2, :) + globalvals(ifv:ifv+2, :) + &
&     globalvals(iflowfm:iflowfm+2, :)
    forcepd = globalvalsd(ifp:ifp+2, :)
    forcep = globalvals(ifp:ifp+2, :)
    forcevd = globalvalsd(ifv:ifv+2, :)
    forcev = globalvals(ifv:ifv+2, :)
    forcemd = globalvalsd(iflowfm:iflowfm+2, :)
    forcem = globalvals(iflowfm:iflowfm+2, :)
    cofxd = globalvalsd(icoforcex:icoforcex+2, :)
    cofx = globalvals(icoforcex:icoforcex+2, :)
    cofyd = globalvalsd(icoforcey:icoforcey+2, :)
    cofy = globalvals(icoforcey:icoforcey+2, :)
    cofzd = globalvalsd(icoforcez:icoforcez+2, :)
    cofz = globalvals(icoforcez:icoforcez+2, :)
    momentd = globalvalsd(imp:imp+2, :) + globalvalsd(imv:imv+2, :) + &
&     globalvalsd(iflowmm:iflowmm+2, :)
    moment = globalvals(imp:imp+2, :) + globalvals(imv:imv+2, :) + &
&     globalvals(iflowmm:iflowmm+2, :)
    temp = gammainf*surfaceref*(lref*lref)
    temp0 = two/(temp*(machcoef*machcoef)*pref)
    factd = -(temp0*(pref*2*machcoef*machcoefd+machcoef**2*prefd)/(&
&     machcoef**2*pref))
    fact = temp0
    cforced = force*factd + fact*forced
    cforce = fact*force
    cforcepd = forcep*factd + fact*forcepd
    cforcep = fact*forcep
    cforcevd = forcev*factd + fact*forcevd
    cforcev = fact*forcev
    cforcemd = forcem*factd + fact*forcemd
    cforcem = fact*forcem
! moment factor has an extra lengthref
    factd = factd/(lengthref*lref)
    fact = fact/(lengthref*lref)
    cmomentd = moment*factd + fact*momentd
    cmoment = fact*moment
! zero values since we are summing.
    funcvalues = zero
    funcvaluesd = 0.0_8
! here we finally assign the final function values
    do sps=1,ntimeintervalsspectral
      funcvaluesd(costfuncforcex) = funcvaluesd(costfuncforcex) + ovrnts&
&       *forced(1, sps)
      funcvalues(costfuncforcex) = funcvalues(costfuncforcex) + ovrnts*&
&       force(1, sps)
      funcvaluesd(costfuncforcey) = funcvaluesd(costfuncforcey) + ovrnts&
&       *forced(2, sps)
      funcvalues(costfuncforcey) = funcvalues(costfuncforcey) + ovrnts*&
&       force(2, sps)
      funcvaluesd(costfuncforcez) = funcvaluesd(costfuncforcez) + ovrnts&
&       *forced(3, sps)
      funcvalues(costfuncforcez) = funcvalues(costfuncforcez) + ovrnts*&
&       force(3, sps)
      funcvaluesd(costfuncforcexpressure) = funcvaluesd(&
&       costfuncforcexpressure) + ovrnts*forcepd(1, sps)
      funcvalues(costfuncforcexpressure) = funcvalues(&
&       costfuncforcexpressure) + ovrnts*forcep(1, sps)
      funcvaluesd(costfuncforceypressure) = funcvaluesd(&
&       costfuncforceypressure) + ovrnts*forcepd(2, sps)
      funcvalues(costfuncforceypressure) = funcvalues(&
&       costfuncforceypressure) + ovrnts*forcep(2, sps)
      funcvaluesd(costfuncforcezpressure) = funcvaluesd(&
&       costfuncforcezpressure) + ovrnts*forcepd(3, sps)
      funcvalues(costfuncforcezpressure) = funcvalues(&
&       costfuncforcezpressure) + ovrnts*forcep(3, sps)
      funcvaluesd(costfuncforcexviscous) = funcvaluesd(&
&       costfuncforcexviscous) + ovrnts*forcevd(1, sps)
      funcvalues(costfuncforcexviscous) = funcvalues(&
&       costfuncforcexviscous) + ovrnts*forcev(1, sps)
      funcvaluesd(costfuncforceyviscous) = funcvaluesd(&
&       costfuncforceyviscous) + ovrnts*forcevd(2, sps)
      funcvalues(costfuncforceyviscous) = funcvalues(&
&       costfuncforceyviscous) + ovrnts*forcev(2, sps)
      funcvaluesd(costfuncforcezviscous) = funcvaluesd(&
&       costfuncforcezviscous) + ovrnts*forcevd(3, sps)
      funcvalues(costfuncforcezviscous) = funcvalues(&
&       costfuncforcezviscous) + ovrnts*forcev(3, sps)
      funcvaluesd(costfuncforcexmomentum) = funcvaluesd(&
&       costfuncforcexmomentum) + ovrnts*forcemd(1, sps)
      funcvalues(costfuncforcexmomentum) = funcvalues(&
&       costfuncforcexmomentum) + ovrnts*forcem(1, sps)
      funcvaluesd(costfuncforceymomentum) = funcvaluesd(&
&       costfuncforceymomentum) + ovrnts*forcemd(2, sps)
      funcvalues(costfuncforceymomentum) = funcvalues(&
&       costfuncforceymomentum) + ovrnts*forcem(2, sps)
      funcvaluesd(costfuncforcezmomentum) = funcvaluesd(&
&       costfuncforcezmomentum) + ovrnts*forcemd(3, sps)
      funcvalues(costfuncforcezmomentum) = funcvalues(&
&       costfuncforcezmomentum) + ovrnts*forcem(3, sps)
! ------------
      funcvaluesd(costfuncforcexcoef) = funcvaluesd(costfuncforcexcoef) &
&       + ovrnts*cforced(1, sps)
      funcvalues(costfuncforcexcoef) = funcvalues(costfuncforcexcoef) + &
&       ovrnts*cforce(1, sps)
      funcvaluesd(costfuncforceycoef) = funcvaluesd(costfuncforceycoef) &
&       + ovrnts*cforced(2, sps)
      funcvalues(costfuncforceycoef) = funcvalues(costfuncforceycoef) + &
&       ovrnts*cforce(2, sps)
      funcvaluesd(costfuncforcezcoef) = funcvaluesd(costfuncforcezcoef) &
&       + ovrnts*cforced(3, sps)
      funcvalues(costfuncforcezcoef) = funcvalues(costfuncforcezcoef) + &
&       ovrnts*cforce(3, sps)
      funcvaluesd(costfuncforcexcoefpressure) = funcvaluesd(&
&       costfuncforcexcoefpressure) + ovrnts*cforcepd(1, sps)
      funcvalues(costfuncforcexcoefpressure) = funcvalues(&
&       costfuncforcexcoefpressure) + ovrnts*cforcep(1, sps)
      funcvaluesd(costfuncforceycoefpressure) = funcvaluesd(&
&       costfuncforceycoefpressure) + ovrnts*cforcepd(2, sps)
      funcvalues(costfuncforceycoefpressure) = funcvalues(&
&       costfuncforceycoefpressure) + ovrnts*cforcep(2, sps)
      funcvaluesd(costfuncforcezcoefpressure) = funcvaluesd(&
&       costfuncforcezcoefpressure) + ovrnts*cforcepd(3, sps)
      funcvalues(costfuncforcezcoefpressure) = funcvalues(&
&       costfuncforcezcoefpressure) + ovrnts*cforcep(3, sps)
      funcvaluesd(costfuncforcexcoefviscous) = funcvaluesd(&
&       costfuncforcexcoefviscous) + ovrnts*cforcevd(1, sps)
      funcvalues(costfuncforcexcoefviscous) = funcvalues(&
&       costfuncforcexcoefviscous) + ovrnts*cforcev(1, sps)
      funcvaluesd(costfuncforceycoefviscous) = funcvaluesd(&
&       costfuncforceycoefviscous) + ovrnts*cforcevd(2, sps)
      funcvalues(costfuncforceycoefviscous) = funcvalues(&
&       costfuncforceycoefviscous) + ovrnts*cforcev(2, sps)
      funcvaluesd(costfuncforcezcoefviscous) = funcvaluesd(&
&       costfuncforcezcoefviscous) + ovrnts*cforcevd(3, sps)
      funcvalues(costfuncforcezcoefviscous) = funcvalues(&
&       costfuncforcezcoefviscous) + ovrnts*cforcev(3, sps)
      funcvaluesd(costfuncforcexcoefmomentum) = funcvaluesd(&
&       costfuncforcexcoefmomentum) + ovrnts*cforcemd(1, sps)
      funcvalues(costfuncforcexcoefmomentum) = funcvalues(&
&       costfuncforcexcoefmomentum) + ovrnts*cforcem(1, sps)
      funcvaluesd(costfuncforceycoefmomentum) = funcvaluesd(&
&       costfuncforceycoefmomentum) + ovrnts*cforcemd(2, sps)
      funcvalues(costfuncforceycoefmomentum) = funcvalues(&
&       costfuncforceycoefmomentum) + ovrnts*cforcem(2, sps)
      funcvaluesd(costfuncforcezcoefmomentum) = funcvaluesd(&
&       costfuncforcezcoefmomentum) + ovrnts*cforcemd(3, 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
        temp1 = cofx(:, sps)/force(1, sps)
        cofxd(:, sps) = (cofxd(:, sps)-temp1*forced(1, sps))/force(1, &
&         sps)
        cofx(:, sps) = temp1
      else
        cofxd(:, sps) = 0.0_8
        cofx(:, sps) = zero
      end if
      if (force(2, sps) .ne. zero) then
        temp1 = cofy(:, sps)/force(2, sps)
        cofyd(:, sps) = (cofyd(:, sps)-temp1*forced(2, sps))/force(2, &
&         sps)
        cofy(:, sps) = temp1
      else
        cofyd(:, sps) = 0.0_8
        cofy(:, sps) = zero
      end if
      if (force(3, sps) .ne. zero) then
        temp1 = cofz(:, sps)/force(3, sps)
        cofzd(:, sps) = (cofzd(:, sps)-temp1*forced(3, sps))/force(3, &
&         sps)
        cofz(:, sps) = temp1
      else
        cofzd(:, sps) = 0.0_8
        cofz(:, sps) = zero
      end if
! fx
      funcvaluesd(costfunccoforcexx) = funcvaluesd(costfunccoforcexx) + &
&       ovrnts*cofxd(1, sps)
      funcvalues(costfunccoforcexx) = funcvalues(costfunccoforcexx) + &
&       ovrnts*cofx(1, sps)
      funcvaluesd(costfunccoforcexy) = funcvaluesd(costfunccoforcexy) + &
&       ovrnts*cofxd(2, sps)
      funcvalues(costfunccoforcexy) = funcvalues(costfunccoforcexy) + &
&       ovrnts*cofx(2, sps)
      funcvaluesd(costfunccoforcexz) = funcvaluesd(costfunccoforcexz) + &
&       ovrnts*cofxd(3, sps)
      funcvalues(costfunccoforcexz) = funcvalues(costfunccoforcexz) + &
&       ovrnts*cofx(3, sps)
! fy
      funcvaluesd(costfunccoforceyx) = funcvaluesd(costfunccoforceyx) + &
&       ovrnts*cofyd(1, sps)
      funcvalues(costfunccoforceyx) = funcvalues(costfunccoforceyx) + &
&       ovrnts*cofy(1, sps)
      funcvaluesd(costfunccoforceyy) = funcvaluesd(costfunccoforceyy) + &
&       ovrnts*cofyd(2, sps)
      funcvalues(costfunccoforceyy) = funcvalues(costfunccoforceyy) + &
&       ovrnts*cofy(2, sps)
      funcvaluesd(costfunccoforceyz) = funcvaluesd(costfunccoforceyz) + &
&       ovrnts*cofyd(3, sps)
      funcvalues(costfunccoforceyz) = funcvalues(costfunccoforceyz) + &
&       ovrnts*cofy(3, sps)
! fz
      funcvaluesd(costfunccoforcezx) = funcvaluesd(costfunccoforcezx) + &
&       ovrnts*cofzd(1, sps)
      funcvalues(costfunccoforcezx) = funcvalues(costfunccoforcezx) + &
&       ovrnts*cofz(1, sps)
      funcvaluesd(costfunccoforcezy) = funcvaluesd(costfunccoforcezy) + &
&       ovrnts*cofzd(2, sps)
      funcvalues(costfunccoforcezy) = funcvalues(costfunccoforcezy) + &
&       ovrnts*cofz(2, sps)
      funcvaluesd(costfunccoforcezz) = funcvaluesd(costfunccoforcezz) + &
&       ovrnts*cofzd(3, sps)
      funcvalues(costfunccoforcezz) = funcvalues(costfunccoforcezz) + &
&       ovrnts*cofz(3, sps)
! ------------
      funcvaluesd(costfuncmomx) = funcvaluesd(costfuncmomx) + ovrnts*&
&       momentd(1, sps)
      funcvalues(costfuncmomx) = funcvalues(costfuncmomx) + ovrnts*&
&       moment(1, sps)
      funcvaluesd(costfuncmomy) = funcvaluesd(costfuncmomy) + ovrnts*&
&       momentd(2, sps)
      funcvalues(costfuncmomy) = funcvalues(costfuncmomy) + ovrnts*&
&       moment(2, sps)
      funcvaluesd(costfuncmomz) = funcvaluesd(costfuncmomz) + ovrnts*&
&       momentd(3, sps)
      funcvalues(costfuncmomz) = funcvalues(costfuncmomz) + ovrnts*&
&       moment(3, sps)
      funcvaluesd(costfuncmomxcoef) = funcvaluesd(costfuncmomxcoef) + &
&       ovrnts*cmomentd(1, sps)
      funcvalues(costfuncmomxcoef) = funcvalues(costfuncmomxcoef) + &
&       ovrnts*cmoment(1, sps)
      funcvaluesd(costfuncmomycoef) = funcvaluesd(costfuncmomycoef) + &
&       ovrnts*cmomentd(2, sps)
      funcvalues(costfuncmomycoef) = funcvalues(costfuncmomycoef) + &
&       ovrnts*cmoment(2, sps)
      funcvaluesd(costfuncmomzcoef) = funcvaluesd(costfuncmomzcoef) + &
&       ovrnts*cmomentd(3, 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_compd = ovrnts*globalvalsd(isepsensorks, sps)/(sepsenmaxrho*&
&         globalvals(isepsensorks, sps))
        ks_comp = ovrnts*(sepsenmaxfamily(sps)+log(globalvals(&
&         isepsensorks, sps))/sepsenmaxrho)
        funcvaluesd(costfuncsepsensorks) = funcvaluesd(&
&         costfuncsepsensorks) + ks_compd
        funcvalues(costfuncsepsensorks) = funcvalues(costfuncsepsensorks&
&         ) + ks_comp
        arg1d = sepsensorkssharpness*2*ks_compd
        arg1 = 2*sepsensorkssharpness*(ks_comp+sepsensorksoffset)
        temp0 = one + exp(arg1)
        temp = globalvals(isepsensorksarea, sps)*ks_comp/temp0
        funcvaluesd(costfuncsepsensorksarea) = funcvaluesd(&
&         costfuncsepsensorksarea) + ovrnts*one*(ks_comp*globalvalsd(&
&         isepsensorksarea, sps)+globalvals(isepsensorksarea, sps)*&
&         ks_compd-temp*exp(arg1)*arg1d)/temp0 + ovrnts*globalvalsd(&
&         isepsensorarea, sps)
        funcvalues(costfuncsepsensorksarea) = funcvalues(&
&         costfuncsepsensorksarea) + ovrnts*one*temp + ovrnts*globalvals&
&         (isepsensorarea, sps)
      end if
      funcvaluesd(costfuncsepsensor) = funcvaluesd(costfuncsepsensor) + &
&       ovrnts*globalvalsd(isepsensor, sps)
      funcvalues(costfuncsepsensor) = funcvalues(costfuncsepsensor) + &
&       ovrnts*globalvals(isepsensor, sps)
      funcvaluesd(costfunccavitation) = funcvaluesd(costfunccavitation) &
&       + ovrnts*globalvalsd(icavitation, 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.
        funcvaluesd(costfunccpmin) = funcvaluesd(costfunccpmin) - ovrnts&
&         *globalvalsd(icpmin, sps)/(cpmin_rho*globalvals(icpmin, sps))
        funcvalues(costfunccpmin) = funcvalues(costfunccpmin) + ovrnts*(&
&         cpmin_family(sps)-log(globalvals(icpmin, sps))/cpmin_rho)
      end if
      funcvaluesd(costfuncaxismoment) = funcvaluesd(costfuncaxismoment) &
&       + ovrnts*globalvalsd(iaxismoment, sps)
      funcvalues(costfuncaxismoment) = funcvalues(costfuncaxismoment) + &
&       ovrnts*globalvals(iaxismoment, sps)
      funcvaluesd(costfuncsepsensoravgx) = funcvaluesd(&
&       costfuncsepsensoravgx) + ovrnts*globalvalsd(isepavg, sps)
      funcvalues(costfuncsepsensoravgx) = funcvalues(&
&       costfuncsepsensoravgx) + ovrnts*globalvals(isepavg, sps)
      funcvaluesd(costfuncsepsensoravgy) = funcvaluesd(&
&       costfuncsepsensoravgy) + ovrnts*globalvalsd(isepavg+1, sps)
      funcvalues(costfuncsepsensoravgy) = funcvalues(&
&       costfuncsepsensoravgy) + ovrnts*globalvals(isepavg+1, sps)
      funcvaluesd(costfuncsepsensoravgz) = funcvaluesd(&
&       costfuncsepsensoravgz) + ovrnts*globalvalsd(isepavg+2, sps)
      funcvalues(costfuncsepsensoravgz) = funcvalues(&
&       costfuncsepsensoravgz) + ovrnts*globalvals(isepavg+2, sps)
      funcvaluesd(costfuncarea) = funcvaluesd(costfuncarea) + ovrnts*&
&       globalvalsd(iarea, sps)
      funcvalues(costfuncarea) = funcvalues(costfuncarea) + ovrnts*&
&       globalvals(iarea, sps)
      funcvaluesd(costfuncflowpower) = funcvaluesd(costfuncflowpower) + &
&       ovrnts*globalvalsd(ipower, sps)
      funcvalues(costfuncflowpower) = funcvalues(costfuncflowpower) + &
&       ovrnts*globalvals(ipower, sps)
      funcvaluesd(costfunccperror2) = funcvaluesd(costfunccperror2) + &
&       ovrnts*globalvalsd(icperror2, sps)
      funcvalues(costfunccperror2) = funcvalues(costfunccperror2) + &
&       ovrnts*globalvals(icperror2, sps)
! mass flow like objective
      mflowd = globalvalsd(imassflow, sps)
      mflow = globalvals(imassflow, sps)
      if (mflow .ne. zero) then
        temp0 = globalvals(imassptot, sps)/mflow
        mavgptotd = (globalvalsd(imassptot, sps)-temp0*mflowd)/mflow
        mavgptot = temp0
        temp0 = globalvals(imassttot, sps)/mflow
        mavgttotd = (globalvalsd(imassttot, sps)-temp0*mflowd)/mflow
        mavgttot = temp0
        temp0 = globalvals(imassrho, sps)/mflow
        mavgrhod = (globalvalsd(imassrho, sps)-temp0*mflowd)/mflow
        mavgrho = temp0
        temp0 = globalvals(imassps, sps)/mflow
        mavgpsd = (globalvalsd(imassps, sps)-temp0*mflowd)/mflow
        mavgps = temp0
        temp0 = globalvals(imassmn, sps)/mflow
        mavgmnd = (globalvalsd(imassmn, sps)-temp0*mflowd)/mflow
        mavgmn = temp0
        temp0 = globalvals(imassa, sps)/mflow
        mavgad = (globalvalsd(imassa, sps)-temp0*mflowd)/mflow
        mavga = temp0
        temp0 = globalvals(imassvx, sps)/mflow
        mavgvxd = (globalvalsd(imassvx, sps)-temp0*mflowd)/mflow
        mavgvx = temp0
        temp0 = globalvals(imassvy, sps)/mflow
        mavgvyd = (globalvalsd(imassvy, sps)-temp0*mflowd)/mflow
        mavgvy = temp0
        temp0 = globalvals(imassvz, sps)/mflow
        mavgvzd = (globalvalsd(imassvz, sps)-temp0*mflowd)/mflow
        mavgvz = temp0
        temp0 = globalvals(imassvi, sps)/mflow
        mavgvid = (globalvalsd(imassvi, sps)-temp0*mflowd)/mflow
        mavgvi = temp0
        arg1 = globalvals(imassnx, sps)**2 + globalvals(imassny, sps)**2&
&         + globalvals(imassnz, sps)**2
        mag = sqrt(arg1)
      else
        mavgptot = zero
        mavgttot = zero
        mavgrho = zero
        mavgps = zero
        mavgmn = zero
        mavga = zero
        mavgvx = zero
        mavgvy = zero
        mavgvz = zero
        mavgvi = zero
        mavgvyd = 0.0_8
        mavgvzd = 0.0_8
        mavgttotd = 0.0_8
        mavgad = 0.0_8
        mavgpsd = 0.0_8
        mavgrhod = 0.0_8
        mavgmnd = 0.0_8
        mavgptotd = 0.0_8
        mavgvid = 0.0_8
        mavgvxd = 0.0_8
      end if
! area averaged objectives
      garead = globalvalsd(iarea, sps)
      garea = globalvals(iarea, sps)
      if (garea .ne. zero) then
! area averaged pressure
        temp0 = globalvals(iareaptot, sps)/garea
        funcvaluesd(costfuncaavgptot) = funcvaluesd(costfuncaavgptot) + &
&         ovrnts*(globalvalsd(iareaptot, sps)-temp0*garead)/garea
        funcvalues(costfuncaavgptot) = funcvalues(costfuncaavgptot) + &
&         ovrnts*temp0
        temp0 = globalvals(iareaps, sps)/garea
        funcvaluesd(costfuncaavgps) = funcvaluesd(costfuncaavgps) + &
&         ovrnts*(globalvalsd(iareaps, sps)-temp0*garead)/garea
        funcvalues(costfuncaavgps) = funcvalues(costfuncaavgps) + ovrnts&
&         *temp0
      end if
      funcvaluesd(costfuncmdot) = funcvaluesd(costfuncmdot) + ovrnts*&
&       mflowd
      funcvalues(costfuncmdot) = funcvalues(costfuncmdot) + ovrnts*mflow
      funcvaluesd(costfuncmavgptot) = funcvaluesd(costfuncmavgptot) + &
&       ovrnts*mavgptotd
      funcvalues(costfuncmavgptot) = funcvalues(costfuncmavgptot) + &
&       ovrnts*mavgptot
      funcvaluesd(costfuncmavgttot) = funcvaluesd(costfuncmavgttot) + &
&       ovrnts*mavgttotd
      funcvalues(costfuncmavgttot) = funcvalues(costfuncmavgttot) + &
&       ovrnts*mavgttot
      funcvaluesd(costfuncmavgrho) = funcvaluesd(costfuncmavgrho) + &
&       ovrnts*mavgrhod
      funcvalues(costfuncmavgrho) = funcvalues(costfuncmavgrho) + ovrnts&
&       *mavgrho
      funcvaluesd(costfuncmavgps) = funcvaluesd(costfuncmavgps) + ovrnts&
&       *mavgpsd
      funcvalues(costfuncmavgps) = funcvalues(costfuncmavgps) + ovrnts*&
&       mavgps
      funcvaluesd(costfuncmavgmn) = funcvaluesd(costfuncmavgmn) + ovrnts&
&       *mavgmnd
      funcvalues(costfuncmavgmn) = funcvalues(costfuncmavgmn) + ovrnts*&
&       mavgmn
      funcvaluesd(costfuncmavga) = funcvaluesd(costfuncmavga) + ovrnts*&
&       mavgad
      funcvalues(costfuncmavga) = funcvalues(costfuncmavga) + ovrnts*&
&       mavga
      funcvaluesd(costfuncmavgvx) = funcvaluesd(costfuncmavgvx) + ovrnts&
&       *mavgvxd
      funcvalues(costfuncmavgvx) = funcvalues(costfuncmavgvx) + ovrnts*&
&       mavgvx
      funcvaluesd(costfuncmavgvy) = funcvaluesd(costfuncmavgvy) + ovrnts&
&       *mavgvyd
      funcvalues(costfuncmavgvy) = funcvalues(costfuncmavgvy) + ovrnts*&
&       mavgvy
      funcvaluesd(costfuncmavgvz) = funcvaluesd(costfuncmavgvz) + ovrnts&
&       *mavgvzd
      funcvalues(costfuncmavgvz) = funcvalues(costfuncmavgvz) + ovrnts*&
&       mavgvz
      funcvaluesd(costfuncmavgvi) = funcvaluesd(costfuncmavgvi) + ovrnts&
&       *mavgvid
      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.
    funcvaluesd(costfunclift) = liftdirection(1)*funcvaluesd(&
&     costfuncforcex) + funcvalues(costfuncforcex)*liftdirectiond(1) + &
&     liftdirection(2)*funcvaluesd(costfuncforcey) + funcvalues(&
&     costfuncforcey)*liftdirectiond(2) + liftdirection(3)*funcvaluesd(&
&     costfuncforcez) + funcvalues(costfuncforcez)*liftdirectiond(3)
    funcvalues(costfunclift) = funcvalues(costfuncforcex)*liftdirection(&
&     1) + funcvalues(costfuncforcey)*liftdirection(2) + funcvalues(&
&     costfuncforcez)*liftdirection(3)
    funcvaluesd(costfuncliftpressure) = liftdirection(1)*funcvaluesd(&
&     costfuncforcexpressure) + funcvalues(costfuncforcexpressure)*&
&     liftdirectiond(1) + liftdirection(2)*funcvaluesd(&
&     costfuncforceypressure) + funcvalues(costfuncforceypressure)*&
&     liftdirectiond(2) + liftdirection(3)*funcvaluesd(&
&     costfuncforcezpressure) + funcvalues(costfuncforcezpressure)*&
&     liftdirectiond(3)
    funcvalues(costfuncliftpressure) = funcvalues(costfuncforcexpressure&
&     )*liftdirection(1) + funcvalues(costfuncforceypressure)*&
&     liftdirection(2) + funcvalues(costfuncforcezpressure)*&
&     liftdirection(3)
    funcvaluesd(costfuncliftviscous) = liftdirection(1)*funcvaluesd(&
&     costfuncforcexviscous) + funcvalues(costfuncforcexviscous)*&
&     liftdirectiond(1) + liftdirection(2)*funcvaluesd(&
&     costfuncforceyviscous) + funcvalues(costfuncforceyviscous)*&
&     liftdirectiond(2) + liftdirection(3)*funcvaluesd(&
&     costfuncforcezviscous) + funcvalues(costfuncforcezviscous)*&
&     liftdirectiond(3)
    funcvalues(costfuncliftviscous) = funcvalues(costfuncforcexviscous)*&
&     liftdirection(1) + funcvalues(costfuncforceyviscous)*liftdirection&
&     (2) + funcvalues(costfuncforcezviscous)*liftdirection(3)
    funcvaluesd(costfuncliftmomentum) = liftdirection(1)*funcvaluesd(&
&     costfuncforcexmomentum) + funcvalues(costfuncforcexmomentum)*&
&     liftdirectiond(1) + liftdirection(2)*funcvaluesd(&
&     costfuncforceymomentum) + funcvalues(costfuncforceymomentum)*&
&     liftdirectiond(2) + liftdirection(3)*funcvaluesd(&
&     costfuncforcezmomentum) + funcvalues(costfuncforcezmomentum)*&
&     liftdirectiond(3)
    funcvalues(costfuncliftmomentum) = funcvalues(costfuncforcexmomentum&
&     )*liftdirection(1) + funcvalues(costfuncforceymomentum)*&
&     liftdirection(2) + funcvalues(costfuncforcezmomentum)*&
&     liftdirection(3)
!-----
    funcvaluesd(costfuncdrag) = dragdirection(1)*funcvaluesd(&
&     costfuncforcex) + funcvalues(costfuncforcex)*dragdirectiond(1) + &
&     dragdirection(2)*funcvaluesd(costfuncforcey) + funcvalues(&
&     costfuncforcey)*dragdirectiond(2) + dragdirection(3)*funcvaluesd(&
&     costfuncforcez) + funcvalues(costfuncforcez)*dragdirectiond(3)
    funcvalues(costfuncdrag) = funcvalues(costfuncforcex)*dragdirection(&
&     1) + funcvalues(costfuncforcey)*dragdirection(2) + funcvalues(&
&     costfuncforcez)*dragdirection(3)
    funcvaluesd(costfuncdragpressure) = dragdirection(1)*funcvaluesd(&
&     costfuncforcexpressure) + funcvalues(costfuncforcexpressure)*&
&     dragdirectiond(1) + dragdirection(2)*funcvaluesd(&
&     costfuncforceypressure) + funcvalues(costfuncforceypressure)*&
&     dragdirectiond(2) + dragdirection(3)*funcvaluesd(&
&     costfuncforcezpressure) + funcvalues(costfuncforcezpressure)*&
&     dragdirectiond(3)
    funcvalues(costfuncdragpressure) = funcvalues(costfuncforcexpressure&
&     )*dragdirection(1) + funcvalues(costfuncforceypressure)*&
&     dragdirection(2) + funcvalues(costfuncforcezpressure)*&
&     dragdirection(3)
    funcvaluesd(costfuncdragviscous) = dragdirection(1)*funcvaluesd(&
&     costfuncforcexviscous) + funcvalues(costfuncforcexviscous)*&
&     dragdirectiond(1) + dragdirection(2)*funcvaluesd(&
&     costfuncforceyviscous) + funcvalues(costfuncforceyviscous)*&
&     dragdirectiond(2) + dragdirection(3)*funcvaluesd(&
&     costfuncforcezviscous) + funcvalues(costfuncforcezviscous)*&
&     dragdirectiond(3)
    funcvalues(costfuncdragviscous) = funcvalues(costfuncforcexviscous)*&
&     dragdirection(1) + funcvalues(costfuncforceyviscous)*dragdirection&
&     (2) + funcvalues(costfuncforcezviscous)*dragdirection(3)
    funcvaluesd(costfuncdragmomentum) = dragdirection(1)*funcvaluesd(&
&     costfuncforcexmomentum) + funcvalues(costfuncforcexmomentum)*&
&     dragdirectiond(1) + dragdirection(2)*funcvaluesd(&
&     costfuncforceymomentum) + funcvalues(costfuncforceymomentum)*&
&     dragdirectiond(2) + dragdirection(3)*funcvaluesd(&
&     costfuncforcezmomentum) + funcvalues(costfuncforcezmomentum)*&
&     dragdirectiond(3)
    funcvalues(costfuncdragmomentum) = funcvalues(costfuncforcexmomentum&
&     )*dragdirection(1) + funcvalues(costfuncforceymomentum)*&
&     dragdirection(2) + funcvalues(costfuncforcezmomentum)*&
&     dragdirection(3)
!-----
    funcvaluesd(costfuncliftcoef) = liftdirection(1)*funcvaluesd(&
&     costfuncforcexcoef) + funcvalues(costfuncforcexcoef)*&
&     liftdirectiond(1) + liftdirection(2)*funcvaluesd(&
&     costfuncforceycoef) + funcvalues(costfuncforceycoef)*&
&     liftdirectiond(2) + liftdirection(3)*funcvaluesd(&
&     costfuncforcezcoef) + funcvalues(costfuncforcezcoef)*&
&     liftdirectiond(3)
    funcvalues(costfuncliftcoef) = funcvalues(costfuncforcexcoef)*&
&     liftdirection(1) + funcvalues(costfuncforceycoef)*liftdirection(2)&
&     + funcvalues(costfuncforcezcoef)*liftdirection(3)
    funcvaluesd(costfuncliftcoefpressure) = liftdirection(1)*funcvaluesd&
&     (costfuncforcexcoefpressure) + funcvalues(&
&     costfuncforcexcoefpressure)*liftdirectiond(1) + liftdirection(2)*&
&     funcvaluesd(costfuncforceycoefpressure) + funcvalues(&
&     costfuncforceycoefpressure)*liftdirectiond(2) + liftdirection(3)*&
&     funcvaluesd(costfuncforcezcoefpressure) + funcvalues(&
&     costfuncforcezcoefpressure)*liftdirectiond(3)
    funcvalues(costfuncliftcoefpressure) = funcvalues(&
&     costfuncforcexcoefpressure)*liftdirection(1) + funcvalues(&
&     costfuncforceycoefpressure)*liftdirection(2) + funcvalues(&
&     costfuncforcezcoefpressure)*liftdirection(3)
    funcvaluesd(costfuncliftcoefviscous) = liftdirection(1)*funcvaluesd(&
&     costfuncforcexcoefviscous) + funcvalues(costfuncforcexcoefviscous)&
&     *liftdirectiond(1) + liftdirection(2)*funcvaluesd(&
&     costfuncforceycoefviscous) + funcvalues(costfuncforceycoefviscous)&
&     *liftdirectiond(2) + liftdirection(3)*funcvaluesd(&
&     costfuncforcezcoefviscous) + funcvalues(costfuncforcezcoefviscous)&
&     *liftdirectiond(3)
    funcvalues(costfuncliftcoefviscous) = funcvalues(&
&     costfuncforcexcoefviscous)*liftdirection(1) + funcvalues(&
&     costfuncforceycoefviscous)*liftdirection(2) + funcvalues(&
&     costfuncforcezcoefviscous)*liftdirection(3)
    funcvaluesd(costfuncliftcoefmomentum) = liftdirection(1)*funcvaluesd&
&     (costfuncforcexcoefmomentum) + funcvalues(&
&     costfuncforcexcoefmomentum)*liftdirectiond(1) + liftdirection(2)*&
&     funcvaluesd(costfuncforceycoefmomentum) + funcvalues(&
&     costfuncforceycoefmomentum)*liftdirectiond(2) + liftdirection(3)*&
&     funcvaluesd(costfuncforcezcoefmomentum) + funcvalues(&
&     costfuncforcezcoefmomentum)*liftdirectiond(3)
    funcvalues(costfuncliftcoefmomentum) = funcvalues(&
&     costfuncforcexcoefmomentum)*liftdirection(1) + funcvalues(&
&     costfuncforceycoefmomentum)*liftdirection(2) + funcvalues(&
&     costfuncforcezcoefmomentum)*liftdirection(3)
!-----
    funcvaluesd(costfuncdragcoef) = dragdirection(1)*funcvaluesd(&
&     costfuncforcexcoef) + funcvalues(costfuncforcexcoef)*&
&     dragdirectiond(1) + dragdirection(2)*funcvaluesd(&
&     costfuncforceycoef) + funcvalues(costfuncforceycoef)*&
&     dragdirectiond(2) + dragdirection(3)*funcvaluesd(&
&     costfuncforcezcoef) + funcvalues(costfuncforcezcoef)*&
&     dragdirectiond(3)
    funcvalues(costfuncdragcoef) = funcvalues(costfuncforcexcoef)*&
&     dragdirection(1) + funcvalues(costfuncforceycoef)*dragdirection(2)&
&     + funcvalues(costfuncforcezcoef)*dragdirection(3)
    funcvaluesd(costfuncdragcoefpressure) = dragdirection(1)*funcvaluesd&
&     (costfuncforcexcoefpressure) + funcvalues(&
&     costfuncforcexcoefpressure)*dragdirectiond(1) + dragdirection(2)*&
&     funcvaluesd(costfuncforceycoefpressure) + funcvalues(&
&     costfuncforceycoefpressure)*dragdirectiond(2) + dragdirection(3)*&
&     funcvaluesd(costfuncforcezcoefpressure) + funcvalues(&
&     costfuncforcezcoefpressure)*dragdirectiond(3)
    funcvalues(costfuncdragcoefpressure) = funcvalues(&
&     costfuncforcexcoefpressure)*dragdirection(1) + funcvalues(&
&     costfuncforceycoefpressure)*dragdirection(2) + funcvalues(&
&     costfuncforcezcoefpressure)*dragdirection(3)
    funcvaluesd(costfuncdragcoefviscous) = dragdirection(1)*funcvaluesd(&
&     costfuncforcexcoefviscous) + funcvalues(costfuncforcexcoefviscous)&
&     *dragdirectiond(1) + dragdirection(2)*funcvaluesd(&
&     costfuncforceycoefviscous) + funcvalues(costfuncforceycoefviscous)&
&     *dragdirectiond(2) + dragdirection(3)*funcvaluesd(&
&     costfuncforcezcoefviscous) + funcvalues(costfuncforcezcoefviscous)&
&     *dragdirectiond(3)
    funcvalues(costfuncdragcoefviscous) = funcvalues(&
&     costfuncforcexcoefviscous)*dragdirection(1) + funcvalues(&
&     costfuncforceycoefviscous)*dragdirection(2) + funcvalues(&
&     costfuncforcezcoefviscous)*dragdirection(3)
    funcvaluesd(costfuncdragcoefmomentum) = dragdirection(1)*funcvaluesd&
&     (costfuncforcexcoefmomentum) + funcvalues(&
&     costfuncforcexcoefmomentum)*dragdirectiond(1) + dragdirection(2)*&
&     funcvaluesd(costfuncforceycoefmomentum) + funcvalues(&
&     costfuncforceycoefmomentum)*dragdirectiond(2) + dragdirection(3)*&
&     funcvaluesd(costfuncforcezcoefmomentum) + funcvalues(&
&     costfuncforcezcoefmomentum)*dragdirectiond(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
    fxliftd = liftdirection(1)*funcvaluesd(costfuncforcex) + funcvalues(&
&     costfuncforcex)*liftdirectiond(1)
    fxlift = funcvalues(costfuncforcex)*liftdirection(1)
    fyliftd = liftdirection(2)*funcvaluesd(costfuncforcey) + funcvalues(&
&     costfuncforcey)*liftdirectiond(2)
    fylift = funcvalues(costfuncforcey)*liftdirection(2)
    fzliftd = liftdirection(3)*funcvaluesd(costfuncforcez) + funcvalues(&
&     costfuncforcez)*liftdirectiond(3)
    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
      temp0 = (fxlift*funcvalues(costfunccoforcexx)+fylift*funcvalues(&
&       costfunccoforceyx)+fzlift*funcvalues(costfunccoforcezx))/(fxlift&
&       +fylift+fzlift)
      funcvaluesd(costfunccofliftx) = (funcvalues(costfunccoforcexx)*&
&       fxliftd+fxlift*funcvaluesd(costfunccoforcexx)+funcvalues(&
&       costfunccoforceyx)*fyliftd+fylift*funcvaluesd(costfunccoforceyx)&
&       +funcvalues(costfunccoforcezx)*fzliftd+fzlift*funcvaluesd(&
&       costfunccoforcezx)-temp0*(fxliftd+fyliftd+fzliftd))/(fxlift+&
&       fylift+fzlift)
      funcvalues(costfunccofliftx) = temp0
      temp0 = (fxlift*funcvalues(costfunccoforcexy)+fylift*funcvalues(&
&       costfunccoforceyy)+fzlift*funcvalues(costfunccoforcezy))/(fxlift&
&       +fylift+fzlift)
      funcvaluesd(costfunccoflifty) = (funcvalues(costfunccoforcexy)*&
&       fxliftd+fxlift*funcvaluesd(costfunccoforcexy)+funcvalues(&
&       costfunccoforceyy)*fyliftd+fylift*funcvaluesd(costfunccoforceyy)&
&       +funcvalues(costfunccoforcezy)*fzliftd+fzlift*funcvaluesd(&
&       costfunccoforcezy)-temp0*(fxliftd+fyliftd+fzliftd))/(fxlift+&
&       fylift+fzlift)
      funcvalues(costfunccoflifty) = temp0
      temp0 = (fxlift*funcvalues(costfunccoforcexz)+fylift*funcvalues(&
&       costfunccoforceyz)+fzlift*funcvalues(costfunccoforcezz))/(fxlift&
&       +fylift+fzlift)
      funcvaluesd(costfunccofliftz) = (funcvalues(costfunccoforcexz)*&
&       fxliftd+fxlift*funcvaluesd(costfunccoforcexz)+funcvalues(&
&       costfunccoforceyz)*fyliftd+fylift*funcvaluesd(costfunccoforceyz)&
&       +funcvalues(costfunccoforcezz)*fzliftd+fzlift*funcvaluesd(&
&       costfunccoforcezz)-temp0*(fxliftd+fyliftd+fzliftd))/(fxlift+&
&       fylift+fzlift)
      funcvalues(costfunccofliftz) = temp0
    else
      funcvaluesd(costfunccofliftx) = 0.0_8
      funcvalues(costfunccofliftx) = zero
      funcvaluesd(costfunccoflifty) = 0.0_8
      funcvalues(costfunccoflifty) = zero
      funcvaluesd(costfunccofliftz) = 0.0_8
      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_d
 
  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_d, only : computetsderivatives
    use flowutils_d, 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
    real(kind=realtype) :: arg1
! 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
        arg1 = 2*sepsensorkssharpness*(ks_comp+sepsensorksoffset)
        funcvalues(costfuncsepsensorksarea) = funcvalues(&
&         costfuncsepsensorksarea) + ovrnts*globalvals(isepsensorksarea&
&         , sps)*ks_comp*one/(one+exp(arg1)) + 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
        arg1 = globalvals(imassnx, sps)**2 + globalvals(imassny, sps)**2&
&         + globalvals(imassnz, sps)**2
        mag = sqrt(arg1)
      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 forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: *(*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:in pref:in *xx:in *pp1:in
!                *pp2:in *ssi:in *ww2:in veldirfreestream:in machcoef:in
!                pointref:in *(*viscsubface.tau):in *(*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_d(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_d
    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) :: arg1
    real(kind=realtype) :: arg1d
    real(kind=realtype) :: result1
    real(kind=realtype) :: result1d
    real(kind=realtype) :: temp
    real(kind=realtype) :: temp0
    real(kind=realtype) :: temp1
    real(kind=realtype) :: tempd
    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.
    refpointd = 0.0_8
    refpointd(1) = lref*pointrefd(1)
    refpoint(1) = lref*pointref(1)
    refpointd(2) = lref*pointrefd(2)
    refpoint(2) = lref*pointref(2)
    refpointd(3) = lref*pointrefd(3)
    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
    sepsensoravgd = 0.0_8
    cofsumfxd = 0.0_8
    cofsumfyd = 0.0_8
    cofsumfzd = 0.0_8
    rd = 0.0_8
    cpmin_ks_sumd = 0.0_8
    vd = 0.0_8
    mpaxisd = 0.0_8
    vecttangentiald = 0.0_8
    sepsensorksd = 0.0_8
    sepsensorksaread = 0.0_8
    cperror2d = 0.0_8
    fpd = 0.0_8
    sepsensoraread = 0.0_8
    mpd = 0.0_8
    cavitationd = 0.0_8
    sepsensord = 0.0_8
!
!         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.
      temp = half*(pp2(i, j)+pp1(i, j)) - pinf
      pm1d = fact*(pref*(half*(pp2d(i, j)+pp1d(i, j))-pinfd)+temp*prefd)
      pm1 = fact*(temp*pref)
      temp = two/(gammainf*pinf*(machcoef*machcoef))
      tmpd = -(temp*(machcoef**2*gammainf*pinfd+gammainf*pinf*2*machcoef&
&       *machcoefd)/(gammainf*pinf*machcoef**2))
      tmp = temp
      temp = half*(pp2(i, j)+pp1(i, j)) - pinf
      cpd = tmp*(half*(pp2d(i, j)+pp1d(i, j))-pinfd) + temp*tmpd
      cp = temp*tmp
      cperrord = cpd
      cperror = cp - bcdata(mm)%cptarget(i, j)
      cperror2d = cperror2d + 2*cperror*cperrord
      cperror2 = cperror2 + cperror*cperror
      xcd = fourth*(xxd(i, j, 1)+xxd(i+1, j, 1)+xxd(i, j+1, 1)+xxd(i+1, &
&       j+1, 1)) - refpointd(1)
      xc = fourth*(xx(i, j, 1)+xx(i+1, j, 1)+xx(i, j+1, 1)+xx(i+1, j+1, &
&       1)) - refpoint(1)
      ycd = fourth*(xxd(i, j, 2)+xxd(i+1, j, 2)+xxd(i, j+1, 2)+xxd(i+1, &
&       j+1, 2)) - refpointd(2)
      yc = fourth*(xx(i, j, 2)+xx(i+1, j, 2)+xx(i, j+1, 2)+xx(i+1, j+1, &
&       2)) - refpoint(2)
      zcd = fourth*(xxd(i, j, 3)+xxd(i+1, j, 3)+xxd(i, j+1, 3)+xxd(i+1, &
&       j+1, 3)) - refpointd(3)
      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
      fxd = ssi(i, j, 1)*pm1d + pm1*ssid(i, j, 1)
      fx = pm1*ssi(i, j, 1)
      fyd = ssi(i, j, 2)*pm1d + pm1*ssid(i, j, 2)
      fy = pm1*ssi(i, j, 2)
      fzd = ssi(i, j, 3)*pm1d + pm1*ssid(i, j, 3)
      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
      fpd(1) = fpd(1) + blk*fxd
      fp(1) = fp(1) + fx*blk
      fpd(2) = fpd(2) + blk*fyd
      fp(2) = fp(2) + fy*blk
      fpd(3) = fpd(3) + blk*fzd
      fp(3) = fp(3) + fz*blk
      mxd = fz*ycd + yc*fzd - fy*zcd - zc*fyd
      mx = yc*fz - zc*fy
      myd = fx*zcd + zc*fxd - fz*xcd - xc*fzd
      my = zc*fx - xc*fz
      mzd = fy*xcd + xc*fyd - fx*ycd - yc*fxd
      mz = xc*fy - yc*fx
      mpd(1) = mpd(1) + blk*mxd
      mp(1) = mp(1) + mx*blk
      mpd(2) = mpd(2) + blk*myd
      mp(2) = mp(2) + my*blk
      mpd(3) = mpd(3) + blk*mzd
      mp(3) = mp(3) + mz*blk
! the force integral for the center of pressure computation.
! we need the cell centers wrt origin
      xcod = fourth*(xxd(i, j, 1)+xxd(i+1, j, 1)+xxd(i, j+1, 1)+xxd(i+1&
&       , j+1, 1))
      xco = fourth*(xx(i, j, 1)+xx(i+1, j, 1)+xx(i, j+1, 1)+xx(i+1, j+1&
&       , 1))
      ycod = fourth*(xxd(i, j, 2)+xxd(i+1, j, 2)+xxd(i, j+1, 2)+xxd(i+1&
&       , j+1, 2))
      yco = fourth*(xx(i, j, 2)+xx(i+1, j, 2)+xx(i, j+1, 2)+xx(i+1, j+1&
&       , 2))
      zcod = fourth*(xxd(i, j, 3)+xxd(i+1, j, 3)+xxd(i, j+1, 3)+xxd(i+1&
&       , j+1, 3))
      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
      cofsumfxd(1) = cofsumfxd(1) + blk*(fx*xcod+xco*fxd)
      cofsumfx(1) = cofsumfx(1) + xco*fx*blk
      cofsumfxd(2) = cofsumfxd(2) + blk*(fx*ycod+yco*fxd)
      cofsumfx(2) = cofsumfx(2) + yco*fx*blk
      cofsumfxd(3) = cofsumfxd(3) + blk*(fx*zcod+zco*fxd)
      cofsumfx(3) = cofsumfx(3) + zco*fx*blk
! force-y
      cofsumfyd(1) = cofsumfyd(1) + blk*(fy*xcod+xco*fyd)
      cofsumfy(1) = cofsumfy(1) + xco*fy*blk
      cofsumfyd(2) = cofsumfyd(2) + blk*(fy*ycod+yco*fyd)
      cofsumfy(2) = cofsumfy(2) + yco*fy*blk
      cofsumfyd(3) = cofsumfyd(3) + blk*(fy*zcod+zco*fyd)
      cofsumfy(3) = cofsumfy(3) + zco*fy*blk
! force-z
      cofsumfzd(1) = cofsumfzd(1) + blk*(fz*xcod+xco*fzd)
      cofsumfz(1) = cofsumfz(1) + xco*fz*blk
      cofsumfzd(2) = cofsumfzd(2) + blk*(fz*ycod+yco*fzd)
      cofsumfz(2) = cofsumfz(2) + yco*fz*blk
      cofsumfzd(3) = cofsumfzd(3) + blk*(fz*zcod+zco*fzd)
      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
      rd(1) = fourth*(xxd(i, j, 1)+xxd(i+1, j, 1)+xxd(i, j+1, 1)+xxd(i+1&
&       , j+1, 1))
      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)
      rd(2) = fourth*(xxd(i, j, 2)+xxd(i+1, j, 2)+xxd(i, j+1, 2)+xxd(i+1&
&       , j+1, 2))
      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)
      rd(3) = fourth*(xxd(i, j, 3)+xxd(i+1, j, 3)+xxd(i, j+1, 3)+xxd(i+1&
&       , j+1, 3))
      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)
      arg1 = (axispoints(1, 2)-axispoints(1, 1))**2 + (axispoints(2, 2)-&
&       axispoints(2, 1))**2 + (axispoints(3, 2)-axispoints(3, 1))**2
      l = sqrt(arg1)
      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
      m0xd = fz*rd(2) + r(2)*fzd - fy*rd(3) - r(3)*fyd
      m0x = r(2)*fz - r(3)*fy
      m0yd = fx*rd(3) + r(3)*fxd - fz*rd(1) - r(1)*fzd
      m0y = r(3)*fx - r(1)*fz
      m0zd = fy*rd(1) + r(1)*fyd - fx*rd(2) - r(2)*fxd
      m0z = r(1)*fy - r(2)*fx
      mpaxisd = mpaxisd + blk*(n(1)*m0xd+n(2)*m0yd+n(3)*m0zd)
      mpaxis = mpaxis + (m0x*n(1)+m0y*n(2)+m0z*n(3))*blk
! save the face-based forces and area
      bcdatad(mm)%fp(i, j, 1) = fxd
      bcdata(mm)%fp(i, j, 1) = fx
      bcdatad(mm)%fp(i, j, 2) = fyd
      bcdata(mm)%fp(i, j, 2) = fy
      bcdatad(mm)%fp(i, j, 3) = fzd
      bcdata(mm)%fp(i, j, 3) = fz
      arg1d = 2*ssi(i, j, 1)*ssid(i, j, 1) + 2*ssi(i, j, 2)*ssid(i, j, 2&
&       ) + 2*ssi(i, j, 3)*ssid(i, j, 3)
      arg1 = ssi(i, j, 1)**2 + ssi(i, j, 2)**2 + ssi(i, j, 3)**2
      temp = sqrt(arg1)
      if (arg1 .eq. 0.0_8) then
        cellaread = 0.0_8
      else
        cellaread = arg1d/(2.0*temp)
      end if
      cellarea = temp
      bcdatad(mm)%area(i, j) = cellaread
      bcdata(mm)%area(i, j) = cellarea
! get normalized surface velocity:
      vd(1) = ww2d(i, j, ivx)
      v(1) = ww2(i, j, ivx)
      vd(2) = ww2d(i, j, ivy)
      v(2) = ww2(i, j, ivy)
      vd(3) = ww2d(i, j, ivz)
      v(3) = ww2(i, j, ivz)
      arg1d = 2*v(1)*vd(1) + 2*v(2)*vd(2) + 2*v(3)*vd(3)
      arg1 = v(1)**2 + v(2)**2 + v(3)**2
      temp = sqrt(arg1)
      if (arg1 .eq. 0.0_8) then
        result1d = 0.0_8
      else
        result1d = arg1d/(2.0*temp)
      end if
      result1 = temp
      vd = (vd-v*result1d/(result1+1e-16))/(result1+1e-16)
      v = v/(result1+1e-16)
      if (computesepsensorks) then
! freestream projection over the surface.
        temp = bcdata(mm)%norm(i, j, 1)
        temp0 = bcdata(mm)%norm(i, j, 2)
        temp1 = bcdata(mm)%norm(i, j, 3)
        vectdotproductfsnormald = temp*veldirfreestreamd(1) + temp0*&
&         veldirfreestreamd(2) + temp1*veldirfreestreamd(3)
        vectdotproductfsnormal = temp*veldirfreestream(1) + temp0*&
&         veldirfreestream(2) + temp1*veldirfreestream(3)
! tangential vector on the surface, which is the freestream projected vector
        temp1 = bcdata(mm)%norm(i, j, 1)
        vecttangentiald(1) = veldirfreestreamd(1) - temp1*&
&         vectdotproductfsnormald
        vecttangential(1) = veldirfreestream(1) - temp1*&
&         vectdotproductfsnormal
        temp1 = bcdata(mm)%norm(i, j, 2)
        vecttangentiald(2) = veldirfreestreamd(2) - temp1*&
&         vectdotproductfsnormald
        vecttangential(2) = veldirfreestream(2) - temp1*&
&         vectdotproductfsnormal
        temp1 = bcdata(mm)%norm(i, j, 3)
        vecttangentiald(3) = veldirfreestreamd(3) - temp1*&
&         vectdotproductfsnormald
        vecttangential(3) = veldirfreestream(3) - temp1*&
&         vectdotproductfsnormal
        arg1d = 2*vecttangential(1)*vecttangentiald(1) + 2*&
&         vecttangential(2)*vecttangentiald(2) + 2*vecttangential(3)*&
&         vecttangentiald(3)
        arg1 = vecttangential(1)**2 + vecttangential(2)**2 + &
&         vecttangential(3)**2
        temp1 = sqrt(arg1)
        if (arg1 .eq. 0.0_8) then
          result1d = 0.0_8
        else
          result1d = arg1d/(2.0*temp1)
        end if
        result1 = temp1
        vecttangentiald = (vecttangentiald-vecttangential*result1d/(&
&         result1+1e-16))/(result1+1e-16)
        vecttangential = vecttangential/(result1+1e-16)
! computing separation sensor
! velocity dot products
        sensord = vecttangential(1)*vd(1) + v(1)*vecttangentiald(1) + &
&         vecttangential(2)*vd(2) + v(2)*vecttangentiald(2) + &
&         vecttangential(3)*vd(3) + v(3)*vecttangentiald(3)
        sensor = v(1)*vecttangential(1) + v(2)*vecttangential(2) + v(3)*&
&         vecttangential(3)
! sepsensor value
        temp1 = cos(zero) - cos(degtorad*sepsensorksphi) + 1e-16
        sensord = -(sensord/temp1)
        sensor = (cos(degtorad*sepsensorksphi)-sensor)/temp1
! also do the ks-based spensenor max computation
        call ksaggregationfunction_d(sensor, sensord, sepsenmaxfamily(&
&                              spectralsol), sepsenmaxrho, ks_exponent, &
&                              ks_exponentd)
        sepsensorksaread = sepsensorksaread + blk*cellaread
        sepsensorksarea = sepsensorksarea + blk*cellarea
        arg1d = -(sepsensorkssharpness*2*sensord)
        arg1 = -(2*sepsensorkssharpness*(sensor+sepsensorksoffset))
        temp1 = one + exp(arg1)
        sepsensoraread = blk*one*(cellaread-cellarea*exp(arg1)*arg1d/&
&         temp1)/temp1 + sepsensoraread
        sepsensorarea = blk*one*(cellarea/temp1) + sepsensorarea
        sepsensorksd = sepsensorksd + blk*ks_exponentd
        sepsensorks = sepsensorks + ks_exponent*blk
      end if
! dot product with free stream
      sensord = -(veldirfreestream(1)*vd(1)) - v(1)*veldirfreestreamd(1)&
&       - veldirfreestream(2)*vd(2) - v(2)*veldirfreestreamd(2) - &
&       veldirfreestream(3)*vd(3) - v(3)*veldirfreestreamd(3)
      sensor = -(v(1)*veldirfreestream(1)+v(2)*veldirfreestream(2)+v(3)*&
&       veldirfreestream(3))
!now run through a smooth heaviside function:
      arg1d = -(sepsensorsharpness*2*sensord)
      arg1 = -(2*sepsensorsharpness*(sensor-sepsensoroffset))
      temp1 = one/(one+exp(arg1))
      sensord = -(temp1*exp(arg1)*arg1d/(one+exp(arg1)))
      sensor = temp1
! and integrate over the area of this cell and save, blanking as we go.
      sensord = blk*(cellarea*sensord+sensor*cellaread)
      sensor = sensor*cellarea*blk
      sepsensord = sepsensord + sensord
      sepsensor = sepsensor + sensor
! also accumulate into the sepsensoravg
! x-y-zco are computed above for center of force computations
      sepsensoravgd(1) = sepsensoravgd(1) + xco*sensord + sensor*xcod
      sepsensoravg(1) = sepsensoravg(1) + sensor*xco
      sepsensoravgd(2) = sepsensoravgd(2) + yco*sensord + sensor*ycod
      sepsensoravg(2) = sepsensoravg(2) + sensor*yco
      sepsensoravgd(3) = sepsensoravgd(3) + zco*sensord + sensor*zcod
      sepsensoravg(3) = sepsensoravg(3) + sensor*zco
      if (computecavitation) then
        plocald = pp2d(i, j)
        plocal = pp2(i, j)
        temp1 = two/(gammainf*(machcoef*machcoef))
        tmpd = -(temp1*2*machcoefd/machcoef)
        tmp = temp1
        cpd = (plocal-pinf)*tmpd + tmp*(plocald-pinfd)
        cp = tmp*(plocal-pinf)
        sensor1d = -cpd
        sensor1 = -cp - cavitationnumber
        arg1d = -(cavsensorsharpness*2*sensor1d)
        arg1 = 2*cavsensorsharpness*(-sensor1+cavsensoroffset)
        temp1 = one + exp(arg1)
        temp = sensor1**cavexponent/temp1
        if (sensor1 .le. 0.0_8 .and. (cavexponent .eq. 0.0_8 .or. &
&           cavexponent .ne. int(cavexponent))) then
          tempd = 0.0_8
        else
          tempd = cavexponent*sensor1**(cavexponent-1)*sensor1d
        end if
        sensor1d = (tempd-temp*exp(arg1)*arg1d)/temp1
        sensor1 = temp
        sensor1d = blk*(cellarea*sensor1d+sensor1*cellaread)
        sensor1 = sensor1*cellarea*blk
        cavitationd = cavitationd + sensor1d
        cavitation = cavitation + sensor1
! also do the ks-based cpmin computation
        temp1 = cpmin_rho*(cpmin_family(spectralsol)-cp)
        ks_exponentd = -(exp(temp1)*cpmin_rho*cpd)
        ks_exponent = exp(temp1)
        cpmin_ks_sumd = cpmin_ks_sumd + blk*ks_exponentd
        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
      mvaxisd = 0.0_8
      fvd = 0.0_8
      mvd = 0.0_8
! 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
        tauxxd = viscsubfaced(mm)%tau(i, j, 1)
        tauxx = viscsubface(mm)%tau(i, j, 1)
        tauyyd = viscsubfaced(mm)%tau(i, j, 2)
        tauyy = viscsubface(mm)%tau(i, j, 2)
        tauzzd = viscsubfaced(mm)%tau(i, j, 3)
        tauzz = viscsubface(mm)%tau(i, j, 3)
        tauxyd = viscsubfaced(mm)%tau(i, j, 4)
        tauxy = viscsubface(mm)%tau(i, j, 4)
        tauxzd = viscsubfaced(mm)%tau(i, j, 5)
        tauxz = viscsubface(mm)%tau(i, j, 5)
        tauyzd = viscsubfaced(mm)%tau(i, j, 6)
        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.
        temp1 = tauxx*ssi(i, j, 1) + tauxy*ssi(i, j, 2) + tauxz*ssi(i, j&
&         , 3)
        fxd = -(fact*(pref*(ssi(i, j, 1)*tauxxd+tauxx*ssid(i, j, 1)+ssi(&
&         i, j, 2)*tauxyd+tauxy*ssid(i, j, 2)+ssi(i, j, 3)*tauxzd+tauxz*&
&         ssid(i, j, 3))+temp1*prefd))
        fx = -(fact*(temp1*pref))
        temp1 = tauxy*ssi(i, j, 1) + tauyy*ssi(i, j, 2) + tauyz*ssi(i, j&
&         , 3)
        fyd = -(fact*(pref*(ssi(i, j, 1)*tauxyd+tauxy*ssid(i, j, 1)+ssi(&
&         i, j, 2)*tauyyd+tauyy*ssid(i, j, 2)+ssi(i, j, 3)*tauyzd+tauyz*&
&         ssid(i, j, 3))+temp1*prefd))
        fy = -(fact*(temp1*pref))
        temp1 = tauxz*ssi(i, j, 1) + tauyz*ssi(i, j, 2) + tauzz*ssi(i, j&
&         , 3)
        fzd = -(fact*(pref*(ssi(i, j, 1)*tauxzd+tauxz*ssid(i, j, 1)+ssi(&
&         i, j, 2)*tauyzd+tauyz*ssid(i, j, 2)+ssi(i, j, 3)*tauzzd+tauzz*&
&         ssid(i, j, 3))+temp1*prefd))
        fz = -(fact*(temp1*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.
        xcd = fourth*(xxd(i, j, 1)+xxd(i+1, j, 1)+xxd(i, j+1, 1)+xxd(i+1&
&         , j+1, 1)) - refpointd(1)
        xc = fourth*(xx(i, j, 1)+xx(i+1, j, 1)+xx(i, j+1, 1)+xx(i+1, j+1&
&         , 1)) - refpoint(1)
        ycd = fourth*(xxd(i, j, 2)+xxd(i+1, j, 2)+xxd(i, j+1, 2)+xxd(i+1&
&         , j+1, 2)) - refpointd(2)
        yc = fourth*(xx(i, j, 2)+xx(i+1, j, 2)+xx(i, j+1, 2)+xx(i+1, j+1&
&         , 2)) - refpoint(2)
        zcd = fourth*(xxd(i, j, 3)+xxd(i+1, j, 3)+xxd(i, j+1, 3)+xxd(i+1&
&         , j+1, 3)) - refpointd(3)
        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.
        fvd(1) = fvd(1) + blk*fxd
        fv(1) = fv(1) + fx*blk
        fvd(2) = fvd(2) + blk*fyd
        fv(2) = fv(2) + fy*blk
        fvd(3) = fvd(3) + blk*fzd
        fv(3) = fv(3) + fz*blk
        mxd = fz*ycd + yc*fzd - fy*zcd - zc*fyd
        mx = yc*fz - zc*fy
        myd = fx*zcd + zc*fxd - fz*xcd - xc*fzd
        my = zc*fx - xc*fz
        mzd = fy*xcd + xc*fyd - fx*ycd - yc*fxd
        mz = xc*fy - yc*fx
        mvd(1) = mvd(1) + blk*mxd
        mv(1) = mv(1) + mx*blk
        mvd(2) = mvd(2) + blk*myd
        mv(2) = mv(2) + my*blk
        mvd(3) = mvd(3) + blk*mzd
        mv(3) = mv(3) + mz*blk
! the force integral for the center of pressure computation.
! we need the cell centers wrt origin
        xcod = fourth*(xxd(i, j, 1)+xxd(i+1, j, 1)+xxd(i, j+1, 1)+xxd(i+&
&         1, j+1, 1))
        xco = fourth*(xx(i, j, 1)+xx(i+1, j, 1)+xx(i, j+1, 1)+xx(i+1, j+&
&         1, 1))
        ycod = fourth*(xxd(i, j, 2)+xxd(i+1, j, 2)+xxd(i, j+1, 2)+xxd(i+&
&         1, j+1, 2))
        yco = fourth*(xx(i, j, 2)+xx(i+1, j, 2)+xx(i, j+1, 2)+xx(i+1, j+&
&         1, 2))
        zcod = fourth*(xxd(i, j, 3)+xxd(i+1, j, 3)+xxd(i, j+1, 3)+xxd(i+&
&         1, j+1, 3))
        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
        cofsumfxd(1) = cofsumfxd(1) + blk*(fx*xcod+xco*fxd)
        cofsumfx(1) = cofsumfx(1) + xco*fx*blk
        cofsumfxd(2) = cofsumfxd(2) + blk*(fx*ycod+yco*fxd)
        cofsumfx(2) = cofsumfx(2) + yco*fx*blk
        cofsumfxd(3) = cofsumfxd(3) + blk*(fx*zcod+zco*fxd)
        cofsumfx(3) = cofsumfx(3) + zco*fx*blk
! force-y
        cofsumfyd(1) = cofsumfyd(1) + blk*(fy*xcod+xco*fyd)
        cofsumfy(1) = cofsumfy(1) + xco*fy*blk
        cofsumfyd(2) = cofsumfyd(2) + blk*(fy*ycod+yco*fyd)
        cofsumfy(2) = cofsumfy(2) + yco*fy*blk
        cofsumfyd(3) = cofsumfyd(3) + blk*(fy*zcod+zco*fyd)
        cofsumfy(3) = cofsumfy(3) + zco*fy*blk
! force-z
        cofsumfzd(1) = cofsumfzd(1) + blk*(fz*xcod+xco*fzd)
        cofsumfz(1) = cofsumfz(1) + xco*fz*blk
        cofsumfzd(2) = cofsumfzd(2) + blk*(fz*ycod+yco*fzd)
        cofsumfz(2) = cofsumfz(2) + yco*fz*blk
        cofsumfzd(3) = cofsumfzd(3) + blk*(fz*zcod+zco*fzd)
        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
        rd(1) = fourth*(xxd(i, j, 1)+xxd(i+1, j, 1)+xxd(i, j+1, 1)+xxd(i&
&         +1, j+1, 1))
        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)
        rd(2) = fourth*(xxd(i, j, 2)+xxd(i+1, j, 2)+xxd(i, j+1, 2)+xxd(i&
&         +1, j+1, 2))
        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)
        rd(3) = fourth*(xxd(i, j, 3)+xxd(i+1, j, 3)+xxd(i, j+1, 3)+xxd(i&
&         +1, j+1, 3))
        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)
        arg1 = (axispoints(1, 2)-axispoints(1, 1))**2 + (axispoints(2, 2&
&         )-axispoints(2, 1))**2 + (axispoints(3, 2)-axispoints(3, 1))**&
&         2
        l = sqrt(arg1)
        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
        m0xd = fz*rd(2) + r(2)*fzd - fy*rd(3) - r(3)*fyd
        m0x = r(2)*fz - r(3)*fy
        m0yd = fx*rd(3) + r(3)*fxd - fz*rd(1) - r(1)*fzd
        m0y = r(3)*fx - r(1)*fz
        m0zd = fy*rd(1) + r(1)*fyd - fx*rd(2) - r(2)*fxd
        m0z = r(1)*fy - r(2)*fx
        mvaxisd = mvaxisd + blk*(n(1)*m0xd+n(2)*m0yd+n(3)*m0zd)
        mvaxis = mvaxis + (m0x*n(1)+m0y*n(2)+m0z*n(3))*blk
! save the face based forces for the slice operations
        bcdatad(mm)%fv(i, j, 1) = fxd
        bcdata(mm)%fv(i, j, 1) = fx
        bcdatad(mm)%fv(i, j, 2) = fyd
        bcdata(mm)%fv(i, j, 2) = fy
        bcdatad(mm)%fv(i, j, 3) = fzd
        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
      bcdatad(mm)%fv = 0.0_8
      bcdata(mm)%fv = zero
      mvaxisd = 0.0_8
      fvd = 0.0_8
      mvd = 0.0_8
    end if
! increment the local values array with the values we computed here.
    localvaluesd(ifp:ifp+2) = localvaluesd(ifp:ifp+2) + fpd
    localvalues(ifp:ifp+2) = localvalues(ifp:ifp+2) + fp
    localvaluesd(ifv:ifv+2) = localvaluesd(ifv:ifv+2) + fvd
    localvalues(ifv:ifv+2) = localvalues(ifv:ifv+2) + fv
    localvaluesd(imp:imp+2) = localvaluesd(imp:imp+2) + mpd
    localvalues(imp:imp+2) = localvalues(imp:imp+2) + mp
    localvaluesd(imv:imv+2) = localvaluesd(imv:imv+2) + mvd
    localvalues(imv:imv+2) = localvalues(imv:imv+2) + mv
    localvaluesd(icoforcex:icoforcex+2) = localvaluesd(icoforcex:&
&     icoforcex+2) + cofsumfxd
    localvalues(icoforcex:icoforcex+2) = localvalues(icoforcex:icoforcex&
&     +2) + cofsumfx
    localvaluesd(icoforcey:icoforcey+2) = localvaluesd(icoforcey:&
&     icoforcey+2) + cofsumfyd
    localvalues(icoforcey:icoforcey+2) = localvalues(icoforcey:icoforcey&
&     +2) + cofsumfy
    localvaluesd(icoforcez:icoforcez+2) = localvaluesd(icoforcez:&
&     icoforcez+2) + cofsumfzd
    localvalues(icoforcez:icoforcez+2) = localvalues(icoforcez:icoforcez&
&     +2) + cofsumfz
    localvaluesd(isepsensor) = localvaluesd(isepsensor) + sepsensord
    localvalues(isepsensor) = localvalues(isepsensor) + sepsensor
    localvaluesd(isepsensorks) = localvaluesd(isepsensorks) + &
&     sepsensorksd
    localvalues(isepsensorks) = localvalues(isepsensorks) + sepsensorks
    localvaluesd(isepsensorksarea) = localvaluesd(isepsensorksarea) + &
&     sepsensorksaread
    localvalues(isepsensorksarea) = localvalues(isepsensorksarea) + &
&     sepsensorksarea
    localvaluesd(isepsensorarea) = localvaluesd(isepsensorarea) + &
&     sepsensoraread
    localvalues(isepsensorarea) = localvalues(isepsensorarea) + &
&     sepsensorarea
    localvaluesd(icavitation) = localvaluesd(icavitation) + cavitationd
    localvalues(icavitation) = localvalues(icavitation) + cavitation
    localvaluesd(icpmin) = localvaluesd(icpmin) + cpmin_ks_sumd
    localvalues(icpmin) = localvalues(icpmin) + cpmin_ks_sum
    localvaluesd(isepavg:isepavg+2) = localvaluesd(isepavg:isepavg+2) + &
&     sepsensoravgd
    localvalues(isepavg:isepavg+2) = localvalues(isepavg:isepavg+2) + &
&     sepsensoravg
    localvaluesd(iaxismoment) = localvaluesd(iaxismoment) + mpaxisd + &
&     mvaxisd
    localvalues(iaxismoment) = localvalues(iaxismoment) + mpaxis + &
&     mvaxis
    localvaluesd(icperror2) = localvaluesd(icperror2) + cperror2d
    localvalues(icperror2) = localvalues(icperror2) + cperror2
  end subroutine wallintegrationface_d
 
  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_d
    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
    real(kind=realtype) :: arg1
    real(kind=realtype) :: result1
    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)
      arg1 = (axispoints(1, 2)-axispoints(1, 1))**2 + (axispoints(2, 2)-&
&       axispoints(2, 1))**2 + (axispoints(3, 2)-axispoints(3, 1))**2
      l = sqrt(arg1)
      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
      arg1 = ssi(i, j, 1)**2 + ssi(i, j, 2)**2 + ssi(i, j, 3)**2
      cellarea = sqrt(arg1)
      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)
      arg1 = v(1)**2 + v(2)**2 + v(3)**2
      result1 = sqrt(arg1)
      v = v/(result1+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)
        arg1 = vecttangential(1)**2 + vecttangential(2)**2 + &
&         vecttangential(3)**2
        result1 = sqrt(arg1)
        vecttangential = vecttangential/(result1+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
        arg1 = -(2*sepsensorkssharpness*(sensor+sepsensorksoffset))
        sepsensorarea = cellarea*blk*one/(one+exp(arg1)) + 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:
      arg1 = -(2*sepsensorsharpness*(sensor-sepsensoroffset))
      sensor = one/(one+exp(arg1))
! 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
        arg1 = 2*cavsensorsharpness*(-sensor1+cavsensoroffset)
        sensor1 = sensor1**cavexponent/(one+exp(arg1))
        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)
        arg1 = (axispoints(1, 2)-axispoints(1, 1))**2 + (axispoints(2, 2&
&         )-axispoints(2, 1))**2 + (axispoints(3, 2)-axispoints(3, 1))**&
&         2
        l = sqrt(arg1)
        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 forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: ks_g
!   with respect to varying inputs: g
  subroutine ksaggregationfunction_d(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), intent(in) :: gd
    intrinsic exp
    ks_gd = exp(g_rho*(g-max_g))*g_rho*gd
    ks_g = exp(g_rho*(g-max_g))
  end subroutine ksaggregationfunction_d
 
  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 forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: 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:in tref:in pref:in rgas:in
!                rhoref:in *xx:in *pp1:in *pp2:in *ssi:in *ww1:in
!                *ww2:in pointref:in localvalues:in-out
!   plus diff mem management of: xx:in pp1:in pp2:in ssi:in ww1:in
!                ww2:in
  subroutine flowintegrationface_d(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_d, only : computeptot, computeptot_d, computettot, &
&   computettot_d
    use bcpointers_d, only : ssi, ssid, sface, ww1, ww1d, ww2, ww2d, pp1&
&   , pp1d, pp2, pp2d, xx, xxd, gamma1, gamma2
    use utils_d, 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) :: arg1
    real(kind=realtype) :: arg1d
    real(kind=realtype) :: result1
    real(kind=realtype) :: result1d
    real(kind=realtype) :: temp
    real(kind=realtype) :: tempd
    real(kind=realtype) :: temp0
    refpointd = 0.0_8
    refpointd(1) = lref*pointrefd(1)
    refpoint(1) = lref*pointref(1)
    refpointd(2) = lref*pointrefd(2)
    refpoint(2) = lref*pointref(2)
    refpointd(3) = lref*pointrefd(3)
    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
    temp = sqrt(pref*rhoref)
    if (pref*rhoref .eq. 0.0_8) then
      mredimd = 0.0_8
    else
      mredimd = (rhoref*prefd+pref*rhorefd)/(2.0*temp)
    end if
    mredim = temp
    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
    mass_ptotd = 0.0_8
    mass_vid = 0.0_8
    cofsumfxd = 0.0_8
    cofsumfyd = 0.0_8
    cofsumfzd = 0.0_8
    aread = 0.0_8
    mmomd = 0.0_8
    mass_vxd = 0.0_8
    mass_vyd = 0.0_8
    mass_ad = 0.0_8
    mass_vzd = 0.0_8
    ptotd = 0.0_8
    mass_psd = 0.0_8
    mass_mnd = 0.0_8
    sfacecoordrefd = 0.0_8
    area_ptotd = 0.0_8
    mass_rhod = 0.0_8
    mass_ttotd = 0.0_8
    mass_nxd = 0.0_8
    fpd = 0.0_8
    mass_nyd = 0.0_8
    mass_nzd = 0.0_8
    fmomd = 0.0_8
    area_psd = 0.0_8
    ttotd = 0.0_8
    massflowrated = 0.0_8
    mpd = 0.0_8
    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
      vxmd = half*(ww1d(i, j, ivx)+ww2d(i, j, ivx))
      vxm = half*(ww1(i, j, ivx)+ww2(i, j, ivx))
      vymd = half*(ww1d(i, j, ivy)+ww2d(i, j, ivy))
      vym = half*(ww1(i, j, ivy)+ww2(i, j, ivy))
      vzmd = half*(ww1d(i, j, ivz)+ww2d(i, j, ivz))
      vzm = half*(ww1(i, j, ivz)+ww2(i, j, ivz))
      rhomd = half*(ww1d(i, j, irho)+ww2d(i, j, irho))
      rhom = half*(ww1(i, j, irho)+ww2(i, j, irho))
      pmd = half*(pp1d(i, j)+pp2d(i, j))
      pm = half*(pp1(i, j)+pp2(i, j))
      gammam = half*(gamma1(i, j)+gamma2(i, j))
      vnmd = ssi(i, j, 1)*vxmd + vxm*ssid(i, j, 1) + ssi(i, j, 2)*vymd +&
&       vym*ssid(i, j, 2) + ssi(i, j, 3)*vzmd + vzm*ssid(i, j, 3)
      vnm = vxm*ssi(i, j, 1) + vym*ssi(i, j, 2) + vzm*ssi(i, j, 3) - sf
      arg1d = 2*vxm*vxmd + 2*vym*vymd + 2*vzm*vzmd
      arg1 = vxm**2 + vym**2 + vzm**2
      temp = sqrt(arg1)
      if (arg1 .eq. 0.0_8) then
        result1d = 0.0_8
      else
        result1d = arg1d/(2.0*temp)
      end if
      result1 = temp
      vmagd = result1d
      vmag = result1 - sf
      arg1d = gammam*(pmd-pm*rhomd/rhom)/rhom
      arg1 = gammam*pm/rhom
      temp = sqrt(arg1)
      if (arg1 .eq. 0.0_8) then
        amd = 0.0_8
      else
        amd = arg1d/(2.0*temp)
      end if
      am = temp
      mnmd = (vmagd-vmag*amd/am)/am
      mnm = vmag/am
      arg1d = 2*ssi(i, j, 1)*ssid(i, j, 1) + 2*ssi(i, j, 2)*ssid(i, j, 2&
&       ) + 2*ssi(i, j, 3)*ssid(i, j, 3)
      arg1 = ssi(i, j, 1)**2 + ssi(i, j, 2)**2 + ssi(i, j, 3)**2
      temp = sqrt(arg1)
      if (arg1 .eq. 0.0_8) then
        cellaread = 0.0_8
      else
        cellaread = arg1d/(2.0*temp)
      end if
      cellarea = temp
      aread = aread + blk*cellaread
      area = area + cellarea*blk
      overcellaread = -(cellaread/cellarea**2)
      overcellarea = 1/cellarea
      call computeptot_d(rhom, rhomd, vxm, vxmd, vym, vymd, vzm, vzmd, &
&                  pm, pmd, ptot, ptotd)
      call computettot_d(rhom, rhomd, vxm, vxmd, vym, vymd, vzm, vzmd, &
&                  pm, pmd, ttot, ttotd)
      massflowratelocald = blk*fact*(mredim*(vnm*rhomd+rhom*vnmd)+rhom*&
&       vnm*mredimd)
      massflowratelocal = rhom*vnm*blk*fact*mredim
      massflowrated = massflowrated + massflowratelocald
      massflowrate = massflowrate + massflowratelocal
! re-dimentionalize quantities
      pmd = pref*pmd + pm*prefd
      pm = pm*pref
      mass_ptotd = mass_ptotd + pref*(massflowratelocal*ptotd+ptot*&
&       massflowratelocald) + ptot*massflowratelocal*prefd
      mass_ptot = mass_ptot + ptot*massflowratelocal*pref
      mass_ttotd = mass_ttotd + tref*(massflowratelocal*ttotd+ttot*&
&       massflowratelocald) + ttot*massflowratelocal*trefd
      mass_ttot = mass_ttot + ttot*massflowratelocal*tref
      mass_rhod = mass_rhod + rhoref*(massflowratelocal*rhomd+rhom*&
&       massflowratelocald) + rhom*massflowratelocal*rhorefd
      mass_rho = mass_rho + rhom*massflowratelocal*rhoref
      mass_ad = mass_ad + uref*(massflowratelocal*amd+am*&
&       massflowratelocald)
      mass_a = mass_a + am*massflowratelocal*uref
      mass_psd = mass_psd + massflowratelocal*pmd + pm*&
&       massflowratelocald
      mass_ps = mass_ps + pm*massflowratelocal
      mass_mnd = mass_mnd + massflowratelocal*mnmd + mnm*&
&       massflowratelocald
      mass_mn = mass_mn + mnm*massflowratelocal
      area_ptotd = area_ptotd + blk*(cellarea*(pref*ptotd+ptot*prefd)+&
&       ptot*pref*cellaread)
      area_ptot = area_ptot + ptot*pref*cellarea*blk
      area_psd = area_psd + blk*(cellarea*pmd+pm*cellaread)
      area_ps = area_ps + pm*cellarea*blk
      sfacecoordrefd(1) = sf*(overcellarea*ssid(i, j, 1)+ssi(i, j, 1)*&
&       overcellaread)
      sfacecoordref(1) = sf*ssi(i, j, 1)*overcellarea
      sfacecoordrefd(2) = sf*(overcellarea*ssid(i, j, 2)+ssi(i, j, 2)*&
&       overcellaread)
      sfacecoordref(2) = sf*ssi(i, j, 2)*overcellarea
      sfacecoordrefd(3) = sf*(overcellarea*ssid(i, j, 3)+ssi(i, j, 3)*&
&       overcellaread)
      sfacecoordref(3) = sf*ssi(i, j, 3)*overcellarea
      temp = uref*vxm - sfacecoordref(1)
      mass_vxd = mass_vxd + massflowratelocal*(uref*vxmd-sfacecoordrefd(&
&       1)) + temp*massflowratelocald
      mass_vx = mass_vx + temp*massflowratelocal
      temp = uref*vym - sfacecoordref(2)
      mass_vyd = mass_vyd + massflowratelocal*(uref*vymd-sfacecoordrefd(&
&       2)) + temp*massflowratelocald
      mass_vy = mass_vy + temp*massflowratelocal
      temp = uref*vzm - sfacecoordref(3)
      mass_vzd = mass_vzd + massflowratelocal*(uref*vzmd-sfacecoordrefd(&
&       3)) + temp*massflowratelocald
      mass_vz = mass_vz + temp*massflowratelocal
      govgm1 = gammainf/(gammainf-one)
      gm1ovg = one/govgm1
      viconst = two*govgm1*rgasdim
      if (one .gt. one/ptot) then
        pratiod = -(one*ptotd/ptot**2)
        pratio = one/ptot
      else
        pratio = one
        pratiod = 0.0_8
      end if
      temp0 = one - pratio**gm1ovg
      if (pratio .le. 0.0_8 .and. (gm1ovg .eq. 0.0_8 .or. gm1ovg .ne. &
&         int(gm1ovg))) then
        tempd = 0.0_8
      else
        tempd = gm1ovg*pratio**(gm1ovg-1)*pratiod
      end if
      arg1d = viconst*(temp0*(tref*ttotd+ttot*trefd)-ttot*tref*tempd)
      arg1 = viconst*(temp0*(ttot*tref))
      temp0 = sqrt(arg1)
      if (arg1 .eq. 0.0_8) then
        vilocald = 0.0_8
      else
        vilocald = arg1d/(2.0*temp0)
      end if
      vilocal = temp0
      mass_vid = mass_vid + massflowratelocal*vilocald + vilocal*&
&       massflowratelocald
      mass_vi = mass_vi + vilocal*massflowratelocal
      mass_nxd = mass_nxd + overcellarea*massflowratelocal*ssid(i, j, 1)&
&       + ssi(i, j, 1)*(massflowratelocal*overcellaread+overcellarea*&
&       massflowratelocald)
      mass_nx = mass_nx + ssi(i, j, 1)*overcellarea*massflowratelocal
      mass_nyd = mass_nyd + overcellarea*massflowratelocal*ssid(i, j, 2)&
&       + ssi(i, j, 2)*(massflowratelocal*overcellaread+overcellarea*&
&       massflowratelocald)
      mass_ny = mass_ny + ssi(i, j, 2)*overcellarea*massflowratelocal
      mass_nzd = mass_nzd + overcellarea*massflowratelocal*ssid(i, j, 3)&
&       + ssi(i, j, 3)*(massflowratelocal*overcellaread+overcellarea*&
&       massflowratelocald)
      mass_nz = mass_nz + ssi(i, j, 3)*overcellarea*massflowratelocal
      xcd = fourth*(xxd(i, j, 1)+xxd(i+1, j, 1)+xxd(i, j+1, 1)+xxd(i+1, &
&       j+1, 1)) - refpointd(1)
      xc = fourth*(xx(i, j, 1)+xx(i+1, j, 1)+xx(i, j+1, 1)+xx(i+1, j+1, &
&       1)) - refpoint(1)
      ycd = fourth*(xxd(i, j, 2)+xxd(i+1, j, 2)+xxd(i, j+1, 2)+xxd(i+1, &
&       j+1, 2)) - refpointd(2)
      yc = fourth*(xx(i, j, 2)+xx(i+1, j, 2)+xx(i, j+1, 2)+xx(i+1, j+1, &
&       2)) - refpoint(2)
      zcd = fourth*(xxd(i, j, 3)+xxd(i+1, j, 3)+xxd(i, j+1, 3)+xxd(i+1, &
&       j+1, 3)) - refpointd(3)
      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.
      pmd = -(fact*blk*(pmd-pinf*prefd))
      pm = -((pm-pinf*pref)*fact*blk)
      fxd = ssi(i, j, 1)*pmd + pm*ssid(i, j, 1)
      fx = pm*ssi(i, j, 1)
      fyd = ssi(i, j, 2)*pmd + pm*ssid(i, j, 2)
      fy = pm*ssi(i, j, 2)
      fzd = ssi(i, j, 3)*pmd + pm*ssid(i, j, 3)
      fz = pm*ssi(i, j, 3)
! update the pressure force and moment coefficients.
      fpd(1) = fpd(1) + fxd
      fp(1) = fp(1) + fx
      fpd(2) = fpd(2) + fyd
      fp(2) = fp(2) + fy
      fpd(3) = fpd(3) + fzd
      fp(3) = fp(3) + fz
      mxd = fz*ycd + yc*fzd - fy*zcd - zc*fyd
      mx = yc*fz - zc*fy
      myd = fx*zcd + zc*fxd - fz*xcd - xc*fzd
      my = zc*fx - xc*fz
      mzd = fy*xcd + xc*fyd - fx*ycd - yc*fxd
      mz = xc*fy - yc*fx
      mpd(1) = mpd(1) + mxd
      mp(1) = mp(1) + mx
      mpd(2) = mpd(2) + myd
      mp(2) = mp(2) + my
      mpd(3) = mpd(3) + mzd
      mp(3) = mp(3) + mz
! the force integral for the center of pressure computation.
! we need the cell centers wrt origin
      xcod = fourth*(xxd(i, j, 1)+xxd(i+1, j, 1)+xxd(i, j+1, 1)+xxd(i+1&
&       , j+1, 1))
      xco = fourth*(xx(i, j, 1)+xx(i+1, j, 1)+xx(i, j+1, 1)+xx(i+1, j+1&
&       , 1))
      ycod = fourth*(xxd(i, j, 2)+xxd(i+1, j, 2)+xxd(i, j+1, 2)+xxd(i+1&
&       , j+1, 2))
      yco = fourth*(xx(i, j, 2)+xx(i+1, j, 2)+xx(i, j+1, 2)+xx(i+1, j+1&
&       , 2))
      zcod = fourth*(xxd(i, j, 3)+xxd(i+1, j, 3)+xxd(i, j+1, 3)+xxd(i+1&
&       , j+1, 3))
      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
      cofsumfxd(1) = cofsumfxd(1) + fx*xcod + xco*fxd
      cofsumfx(1) = cofsumfx(1) + xco*fx
      cofsumfxd(2) = cofsumfxd(2) + fx*ycod + yco*fxd
      cofsumfx(2) = cofsumfx(2) + yco*fx
      cofsumfxd(3) = cofsumfxd(3) + fx*zcod + zco*fxd
      cofsumfx(3) = cofsumfx(3) + zco*fx
! force-y
      cofsumfyd(1) = cofsumfyd(1) + fy*xcod + xco*fyd
      cofsumfy(1) = cofsumfy(1) + xco*fy
      cofsumfyd(2) = cofsumfyd(2) + fy*ycod + yco*fyd
      cofsumfy(2) = cofsumfy(2) + yco*fy
      cofsumfyd(3) = cofsumfyd(3) + fy*zcod + zco*fyd
      cofsumfy(3) = cofsumfy(3) + zco*fy
! force-z
      cofsumfzd(1) = cofsumfzd(1) + fz*xcod + xco*fzd
      cofsumfz(1) = cofsumfz(1) + xco*fz
      cofsumfzd(2) = cofsumfzd(2) + fz*ycod + yco*fzd
      cofsumfz(2) = cofsumfz(2) + yco*fz
      cofsumfzd(3) = cofsumfzd(3) + fz*zcod + zco*fzd
      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
      temp0 = fact*blk*internalflowfact*inflowfact
      temp = massflowratelocal/(timeref*cellarea)
      massflowratelocald = temp0*(massflowratelocald-temp*(cellarea*&
&       timerefd+timeref*cellaread))/(timeref*cellarea)
      massflowratelocal = temp0*temp
      fxd = massflowratelocal*vxm*ssid(i, j, 1) + ssi(i, j, 1)*(vxm*&
&       massflowratelocald+massflowratelocal*vxmd)
      fx = massflowratelocal*ssi(i, j, 1)*vxm
      fyd = massflowratelocal*vym*ssid(i, j, 2) + ssi(i, j, 2)*(vym*&
&       massflowratelocald+massflowratelocal*vymd)
      fy = massflowratelocal*ssi(i, j, 2)*vym
      fzd = massflowratelocal*vzm*ssid(i, j, 3) + ssi(i, j, 3)*(vzm*&
&       massflowratelocald+massflowratelocal*vzmd)
      fz = massflowratelocal*ssi(i, j, 3)*vzm
      fmomd(1) = fmomd(1) + fxd
      fmom(1) = fmom(1) + fx
      fmomd(2) = fmomd(2) + fyd
      fmom(2) = fmom(2) + fy
      fmomd(3) = fmomd(3) + fzd
      fmom(3) = fmom(3) + fz
      mxd = fz*ycd + yc*fzd - fy*zcd - zc*fyd
      mx = yc*fz - zc*fy
      myd = fx*zcd + zc*fxd - fz*xcd - xc*fzd
      my = zc*fx - xc*fz
      mzd = fy*xcd + xc*fyd - fx*ycd - yc*fxd
      mz = xc*fy - yc*fx
      mmomd(1) = mmomd(1) + mxd
      mmom(1) = mmom(1) + mx
      mmomd(2) = mmomd(2) + myd
      mmom(2) = mmom(2) + my
      mmomd(3) = mmomd(3) + mzd
      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
      cofsumfxd(1) = cofsumfxd(1) + fx*xcod + xco*fxd
      cofsumfx(1) = cofsumfx(1) + xco*fx
      cofsumfxd(2) = cofsumfxd(2) + fx*ycod + yco*fxd
      cofsumfx(2) = cofsumfx(2) + yco*fx
      cofsumfxd(3) = cofsumfxd(3) + fx*zcod + zco*fxd
      cofsumfx(3) = cofsumfx(3) + zco*fx
! force-y
      cofsumfyd(1) = cofsumfyd(1) + fy*xcod + xco*fyd
      cofsumfy(1) = cofsumfy(1) + xco*fy
      cofsumfyd(2) = cofsumfyd(2) + fy*ycod + yco*fyd
      cofsumfy(2) = cofsumfy(2) + yco*fy
      cofsumfyd(3) = cofsumfyd(3) + fy*zcod + zco*fyd
      cofsumfy(3) = cofsumfy(3) + zco*fy
! force-z
      cofsumfzd(1) = cofsumfzd(1) + fz*xcod + xco*fzd
      cofsumfz(1) = cofsumfz(1) + xco*fz
      cofsumfzd(2) = cofsumfzd(2) + fz*ycod + yco*fzd
      cofsumfz(2) = cofsumfz(2) + yco*fz
      cofsumfzd(3) = cofsumfzd(3) + fz*zcod + zco*fzd
      cofsumfz(3) = cofsumfz(3) + zco*fz
    end do
! increment the local values array with what we computed here
    localvaluesd(imassflow) = localvaluesd(imassflow) + massflowrated
    localvalues(imassflow) = localvalues(imassflow) + massflowrate
    localvaluesd(iarea) = localvaluesd(iarea) + aread
    localvalues(iarea) = localvalues(iarea) + area
    localvaluesd(imassrho) = localvaluesd(imassrho) + mass_rhod
    localvalues(imassrho) = localvalues(imassrho) + mass_rho
    localvaluesd(imassa) = localvaluesd(imassa) + mass_ad
    localvalues(imassa) = localvalues(imassa) + mass_a
    localvaluesd(imassptot) = localvaluesd(imassptot) + mass_ptotd
    localvalues(imassptot) = localvalues(imassptot) + mass_ptot
    localvaluesd(imassttot) = localvaluesd(imassttot) + mass_ttotd
    localvalues(imassttot) = localvalues(imassttot) + mass_ttot
    localvaluesd(imassps) = localvaluesd(imassps) + mass_psd
    localvalues(imassps) = localvalues(imassps) + mass_ps
    localvaluesd(imassmn) = localvaluesd(imassmn) + mass_mnd
    localvalues(imassmn) = localvalues(imassmn) + mass_mn
    localvaluesd(ifp:ifp+2) = localvaluesd(ifp:ifp+2) + fpd
    localvalues(ifp:ifp+2) = localvalues(ifp:ifp+2) + fp
    localvaluesd(iflowfm:iflowfm+2) = localvaluesd(iflowfm:iflowfm+2) + &
&     fmomd
    localvalues(iflowfm:iflowfm+2) = localvalues(iflowfm:iflowfm+2) + &
&     fmom
    localvaluesd(iflowmp:iflowmp+2) = localvaluesd(iflowmp:iflowmp+2) + &
&     mpd
    localvalues(iflowmp:iflowmp+2) = localvalues(iflowmp:iflowmp+2) + mp
    localvaluesd(iflowmm:iflowmm+2) = localvaluesd(iflowmm:iflowmm+2) + &
&     mmomd
    localvalues(iflowmm:iflowmm+2) = localvalues(iflowmm:iflowmm+2) + &
&     mmom
    localvaluesd(icoforcex:icoforcex+2) = localvaluesd(icoforcex:&
&     icoforcex+2) + cofsumfxd
    localvalues(icoforcex:icoforcex+2) = localvalues(icoforcex:icoforcex&
&     +2) + cofsumfx
    localvaluesd(icoforcey:icoforcey+2) = localvaluesd(icoforcey:&
&     icoforcey+2) + cofsumfyd
    localvalues(icoforcey:icoforcey+2) = localvalues(icoforcey:icoforcey&
&     +2) + cofsumfy
    localvaluesd(icoforcez:icoforcez+2) = localvaluesd(icoforcez:&
&     icoforcez+2) + cofsumfzd
    localvalues(icoforcez:icoforcez+2) = localvalues(icoforcez:icoforcez&
&     +2) + cofsumfz
    localvaluesd(iareaptot) = localvaluesd(iareaptot) + area_ptotd
    localvalues(iareaptot) = localvalues(iareaptot) + area_ptot
    localvaluesd(iareaps) = localvaluesd(iareaps) + area_psd
    localvalues(iareaps) = localvalues(iareaps) + area_ps
    localvaluesd(imassvx) = localvaluesd(imassvx) + mass_vxd
    localvalues(imassvx) = localvalues(imassvx) + mass_vx
    localvaluesd(imassvy) = localvaluesd(imassvy) + mass_vyd
    localvalues(imassvy) = localvalues(imassvy) + mass_vy
    localvaluesd(imassvz) = localvaluesd(imassvz) + mass_vzd
    localvalues(imassvz) = localvalues(imassvz) + mass_vz
    localvaluesd(imassnx) = localvaluesd(imassnx) + mass_nxd
    localvalues(imassnx) = localvalues(imassnx) + mass_nx
    localvaluesd(imassny) = localvaluesd(imassny) + mass_nyd
    localvalues(imassny) = localvalues(imassny) + mass_ny
    localvaluesd(imassnz) = localvaluesd(imassnz) + mass_nzd
    localvalues(imassnz) = localvalues(imassnz) + mass_nz
    localvaluesd(imassvi) = localvaluesd(imassvi) + mass_vid
    localvalues(imassvi) = localvalues(imassvi) + mass_vi
  end subroutine flowintegrationface_d
 
  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_d, only : computeptot, computettot
    use bcpointers_d, only : ssi, sface, ww1, ww2, pp1, pp2, xx, gamma1,&
&   gamma2
    use utils_d, 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
    real(kind=realtype) :: arg1
    real(kind=realtype) :: result1
    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
      arg1 = vxm**2 + vym**2 + vzm**2
      result1 = sqrt(arg1)
      vmag = result1 - sf
      arg1 = gammam*pm/rhom
      am = sqrt(arg1)
      mnm = vmag/am
      arg1 = ssi(i, j, 1)**2 + ssi(i, j, 2)**2 + ssi(i, j, 3)**2
      cellarea = sqrt(arg1)
      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
      arg1 = viconst*(one-pratio**gm1ovg)*ttot*tref
      vilocal = sqrt(arg1)
      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_d