surfaceIntegrations_fast_b.f90

Go to the documentation of this file.

!        generated by tapenade     (inria, ecuador team)
!  tapenade 3.16 (develop) - 22 aug 2023 15:51
!
module surfaceintegrations_fast_b
  implicit none
 
contains
  subroutine getcostfunctions(globalvals, funcvalues)
    use constants
    use inputtimespectral, only : ntimeintervalsspectral
    use flowvarrefstate, only : pref, rhoref, tref, lref, gammainf, pinf&
&   , uref, uinf
    use inputphysics, only : liftdirection, dragdirection, surfaceref, &
&   machcoef, lengthref, alpha, beta, liftindex, cpmin_family, cpmin_rho&
&   , sepsenmaxfamily, sepsenmaxrho
    use inputcostfunctions, only : computecavitation, computesepsensorks&
&   , sepsensorkssharpness, sepsensorksoffset
    use inputtsstabderiv, only : tsstability
    use utils_fast_b, only : computetsderivatives
    use flowutils_fast_b, only : getdirvector
    implicit none
! input/output
    real(kind=realtype), dimension(:, :), intent(in) :: globalvals
    real(kind=realtype), dimension(:), intent(out) :: funcvalues
! working
    real(kind=realtype) :: fact, factmoment, ovrnts
    real(kind=realtype), dimension(3, ntimeintervalsspectral) :: force, &
&   forcep, forcev, forcem, moment, cforce, cforcep, cforcev, cforcem, &
&   cmoment, cofx, cofy, cofz
    real(kind=realtype), dimension(3) :: vcoordref, vfreestreamref
    real(kind=realtype) :: mavgptot, mavgttot, mavgrho, mavgps, mflow, &
&   mavgmn, mavga, mavgvx, mavgvy, mavgvz, garea, mavgvi, fxlift, fylift&
&   , fzlift
    real(kind=realtype) :: vdotn, mag, u, v, w, ks_comp
    integer(kind=inttype) :: sps
    real(kind=realtype), dimension(8) :: dcdq, dcdqdot
    real(kind=realtype), dimension(8) :: dcdalpha, dcdalphadot
    real(kind=realtype), dimension(8) :: coef0
    intrinsic log
    intrinsic exp
    intrinsic sqrt
! factor used for time-averaged quantities.
    ovrnts = one/ntimeintervalsspectral
! sum pressure and viscous contributions
    force = globalvals(ifp:ifp+2, :) + globalvals(ifv:ifv+2, :) + &
&     globalvals(iflowfm:iflowfm+2, :)
    forcep = globalvals(ifp:ifp+2, :)
    forcev = globalvals(ifv:ifv+2, :)
    forcem = globalvals(iflowfm:iflowfm+2, :)
    cofx = globalvals(icoforcex:icoforcex+2, :)
    cofy = globalvals(icoforcey:icoforcey+2, :)
    cofz = globalvals(icoforcez:icoforcez+2, :)
    moment = globalvals(imp:imp+2, :) + globalvals(imv:imv+2, :) + &
&     globalvals(iflowmm:iflowmm+2, :)
    fact = two/(gammainf*machcoef*machcoef*surfaceref*lref*lref*pref)
    cforce = fact*force
    cforcep = fact*forcep
    cforcev = fact*forcev
    cforcem = fact*forcem
! moment factor has an extra lengthref
    fact = fact/(lengthref*lref)
    cmoment = fact*moment
! zero values since we are summing.
    funcvalues = zero
!$ad ii-loop
! here we finally assign the final function values
    do sps=1,ntimeintervalsspectral
      funcvalues(costfuncforcex) = funcvalues(costfuncforcex) + ovrnts*&
&       force(1, sps)
      funcvalues(costfuncforcey) = funcvalues(costfuncforcey) + ovrnts*&
&       force(2, sps)
      funcvalues(costfuncforcez) = funcvalues(costfuncforcez) + ovrnts*&
&       force(3, sps)
      funcvalues(costfuncforcexpressure) = funcvalues(&
&       costfuncforcexpressure) + ovrnts*forcep(1, sps)
      funcvalues(costfuncforceypressure) = funcvalues(&
&       costfuncforceypressure) + ovrnts*forcep(2, sps)
      funcvalues(costfuncforcezpressure) = funcvalues(&
&       costfuncforcezpressure) + ovrnts*forcep(3, sps)
      funcvalues(costfuncforcexviscous) = funcvalues(&
&       costfuncforcexviscous) + ovrnts*forcev(1, sps)
      funcvalues(costfuncforceyviscous) = funcvalues(&
&       costfuncforceyviscous) + ovrnts*forcev(2, sps)
      funcvalues(costfuncforcezviscous) = funcvalues(&
&       costfuncforcezviscous) + ovrnts*forcev(3, sps)
      funcvalues(costfuncforcexmomentum) = funcvalues(&
&       costfuncforcexmomentum) + ovrnts*forcem(1, sps)
      funcvalues(costfuncforceymomentum) = funcvalues(&
&       costfuncforceymomentum) + ovrnts*forcem(2, sps)
      funcvalues(costfuncforcezmomentum) = funcvalues(&
&       costfuncforcezmomentum) + ovrnts*forcem(3, sps)
! ------------
      funcvalues(costfuncforcexcoef) = funcvalues(costfuncforcexcoef) + &
&       ovrnts*cforce(1, sps)
      funcvalues(costfuncforceycoef) = funcvalues(costfuncforceycoef) + &
&       ovrnts*cforce(2, sps)
      funcvalues(costfuncforcezcoef) = funcvalues(costfuncforcezcoef) + &
&       ovrnts*cforce(3, sps)
      funcvalues(costfuncforcexcoefpressure) = funcvalues(&
&       costfuncforcexcoefpressure) + ovrnts*cforcep(1, sps)
      funcvalues(costfuncforceycoefpressure) = funcvalues(&
&       costfuncforceycoefpressure) + ovrnts*cforcep(2, sps)
      funcvalues(costfuncforcezcoefpressure) = funcvalues(&
&       costfuncforcezcoefpressure) + ovrnts*cforcep(3, sps)
      funcvalues(costfuncforcexcoefviscous) = funcvalues(&
&       costfuncforcexcoefviscous) + ovrnts*cforcev(1, sps)
      funcvalues(costfuncforceycoefviscous) = funcvalues(&
&       costfuncforceycoefviscous) + ovrnts*cforcev(2, sps)
      funcvalues(costfuncforcezcoefviscous) = funcvalues(&
&       costfuncforcezcoefviscous) + ovrnts*cforcev(3, sps)
      funcvalues(costfuncforcexcoefmomentum) = funcvalues(&
&       costfuncforcexcoefmomentum) + ovrnts*cforcem(1, sps)
      funcvalues(costfuncforceycoefmomentum) = funcvalues(&
&       costfuncforceycoefmomentum) + ovrnts*cforcem(2, sps)
      funcvalues(costfuncforcezcoefmomentum) = funcvalues(&
&       costfuncforcezcoefmomentum) + ovrnts*cforcem(3, sps)
! ------------
! center of pressure (these are actually center of all forces)
! protect the divisions against zero, and divide the weighed sum by the force magnitude
! for this time spectral instance before we add it to the sum
      if (force(1, sps) .ne. zero) then
        cofx(:, sps) = cofx(:, sps)/force(1, sps)
      else
        cofx(:, sps) = zero
      end if
      if (force(2, sps) .ne. zero) then
        cofy(:, sps) = cofy(:, sps)/force(2, sps)
      else
        cofy(:, sps) = zero
      end if
      if (force(3, sps) .ne. zero) then
        cofz(:, sps) = cofz(:, sps)/force(3, sps)
      else
        cofz(:, sps) = zero
      end if
! fx
      funcvalues(costfunccoforcexx) = funcvalues(costfunccoforcexx) + &
&       ovrnts*cofx(1, sps)
      funcvalues(costfunccoforcexy) = funcvalues(costfunccoforcexy) + &
&       ovrnts*cofx(2, sps)
      funcvalues(costfunccoforcexz) = funcvalues(costfunccoforcexz) + &
&       ovrnts*cofx(3, sps)
! fy
      funcvalues(costfunccoforceyx) = funcvalues(costfunccoforceyx) + &
&       ovrnts*cofy(1, sps)
      funcvalues(costfunccoforceyy) = funcvalues(costfunccoforceyy) + &
&       ovrnts*cofy(2, sps)
      funcvalues(costfunccoforceyz) = funcvalues(costfunccoforceyz) + &
&       ovrnts*cofy(3, sps)
! fz
      funcvalues(costfunccoforcezx) = funcvalues(costfunccoforcezx) + &
&       ovrnts*cofz(1, sps)
      funcvalues(costfunccoforcezy) = funcvalues(costfunccoforcezy) + &
&       ovrnts*cofz(2, sps)
      funcvalues(costfunccoforcezz) = funcvalues(costfunccoforcezz) + &
&       ovrnts*cofz(3, sps)
! ------------
      funcvalues(costfuncmomx) = funcvalues(costfuncmomx) + ovrnts*&
&       moment(1, sps)
      funcvalues(costfuncmomy) = funcvalues(costfuncmomy) + ovrnts*&
&       moment(2, sps)
      funcvalues(costfuncmomz) = funcvalues(costfuncmomz) + ovrnts*&
&       moment(3, sps)
      funcvalues(costfuncmomxcoef) = funcvalues(costfuncmomxcoef) + &
&       ovrnts*cmoment(1, sps)
      funcvalues(costfuncmomycoef) = funcvalues(costfuncmomycoef) + &
&       ovrnts*cmoment(2, sps)
      funcvalues(costfuncmomzcoef) = funcvalues(costfuncmomzcoef) + &
&       ovrnts*cmoment(3, sps)
! final part of the ks computation
      if (computesepsensorks) then
! only calculate the log part if we are actually computing for separation for ks method.
        ks_comp = ovrnts*(sepsenmaxfamily(sps)+log(globalvals(&
&         isepsensorks, sps))/sepsenmaxrho)
        funcvalues(costfuncsepsensorks) = funcvalues(costfuncsepsensorks&
&         ) + ks_comp
        funcvalues(costfuncsepsensorksarea) = funcvalues(&
&         costfuncsepsensorksarea) + ovrnts*globalvals(isepsensorksarea&
&         , sps)*ks_comp*one/(one+exp(2*sepsensorkssharpness*(ks_comp+&
&         sepsensorksoffset))) + ovrnts*globalvals(isepsensorarea, sps)
      end if
      funcvalues(costfuncsepsensor) = funcvalues(costfuncsepsensor) + &
&       ovrnts*globalvals(isepsensor, sps)
      funcvalues(costfunccavitation) = funcvalues(costfunccavitation) + &
&       ovrnts*globalvals(icavitation, sps)
! final part of the ks computation
      if (computecavitation) then
! only calculate the log part if we are actually computing for cavitation.
! if we are not computing cavitation, the icpmin in globalvals will be zero,
! which doesn't play well with log. we just want to return zero here.
        funcvalues(costfunccpmin) = funcvalues(costfunccpmin) + ovrnts*(&
&         cpmin_family(sps)-log(globalvals(icpmin, sps))/cpmin_rho)
      end if
      funcvalues(costfuncaxismoment) = funcvalues(costfuncaxismoment) + &
&       ovrnts*globalvals(iaxismoment, sps)
      funcvalues(costfuncsepsensoravgx) = funcvalues(&
&       costfuncsepsensoravgx) + ovrnts*globalvals(isepavg, sps)
      funcvalues(costfuncsepsensoravgy) = funcvalues(&
&       costfuncsepsensoravgy) + ovrnts*globalvals(isepavg+1, sps)
      funcvalues(costfuncsepsensoravgz) = funcvalues(&
&       costfuncsepsensoravgz) + ovrnts*globalvals(isepavg+2, sps)
      funcvalues(costfuncarea) = funcvalues(costfuncarea) + ovrnts*&
&       globalvals(iarea, sps)
      funcvalues(costfuncflowpower) = funcvalues(costfuncflowpower) + &
&       ovrnts*globalvals(ipower, sps)
      funcvalues(costfunccperror2) = funcvalues(costfunccperror2) + &
&       ovrnts*globalvals(icperror2, sps)
! mass flow like objective
      mflow = globalvals(imassflow, sps)
      if (mflow .ne. zero) then
        mavgptot = globalvals(imassptot, sps)/mflow
        mavgttot = globalvals(imassttot, sps)/mflow
        mavgrho = globalvals(imassrho, sps)/mflow
        mavgps = globalvals(imassps, sps)/mflow
        mavgmn = globalvals(imassmn, sps)/mflow
        mavga = globalvals(imassa, sps)/mflow
        mavgvx = globalvals(imassvx, sps)/mflow
        mavgvy = globalvals(imassvy, sps)/mflow
        mavgvz = globalvals(imassvz, sps)/mflow
        mavgvi = globalvals(imassvi, sps)/mflow
        mag = sqrt(globalvals(imassnx, sps)**2 + globalvals(imassny, sps&
&         )**2 + globalvals(imassnz, sps)**2)
      else
        mavgptot = zero
        mavgttot = zero
        mavgrho = zero
        mavgps = zero
        mavgmn = zero
        mavga = zero
        mavgvx = zero
        mavgvy = zero
        mavgvz = zero
        mavgvi = zero
      end if
! area averaged objectives
      garea = globalvals(iarea, sps)
      if (garea .ne. zero) then
! area averaged pressure
        funcvalues(costfuncaavgptot) = funcvalues(costfuncaavgptot) + &
&         ovrnts*globalvals(iareaptot, sps)/garea
        funcvalues(costfuncaavgps) = funcvalues(costfuncaavgps) + ovrnts&
&         *globalvals(iareaps, sps)/garea
      end if
      funcvalues(costfuncmdot) = funcvalues(costfuncmdot) + ovrnts*mflow
      funcvalues(costfuncmavgptot) = funcvalues(costfuncmavgptot) + &
&       ovrnts*mavgptot
      funcvalues(costfuncmavgttot) = funcvalues(costfuncmavgttot) + &
&       ovrnts*mavgttot
      funcvalues(costfuncmavgrho) = funcvalues(costfuncmavgrho) + ovrnts&
&       *mavgrho
      funcvalues(costfuncmavgps) = funcvalues(costfuncmavgps) + ovrnts*&
&       mavgps
      funcvalues(costfuncmavgmn) = funcvalues(costfuncmavgmn) + ovrnts*&
&       mavgmn
      funcvalues(costfuncmavga) = funcvalues(costfuncmavga) + ovrnts*&
&       mavga
      funcvalues(costfuncmavgvx) = funcvalues(costfuncmavgvx) + ovrnts*&
&       mavgvx
      funcvalues(costfuncmavgvy) = funcvalues(costfuncmavgvy) + ovrnts*&
&       mavgvy
      funcvalues(costfuncmavgvz) = funcvalues(costfuncmavgvz) + ovrnts*&
&       mavgvz
      funcvalues(costfuncmavgvi) = funcvalues(costfuncmavgvi) + ovrnts*&
&       mavgvi
! bending moment calc - also broken.
! call computerootbendingmoment(cforce, cmoment, liftindex, bendingmoment)
! funcvalues(costfuncbendingcoef) = funcvalues(costfuncbendingcoef) + ovrnts*bendingmoment
    end do
! lift and drag (coefficients): dot product with the lift/drag direction.
    funcvalues(costfunclift) = funcvalues(costfuncforcex)*liftdirection(&
&     1) + funcvalues(costfuncforcey)*liftdirection(2) + funcvalues(&
&     costfuncforcez)*liftdirection(3)
    funcvalues(costfuncliftpressure) = funcvalues(costfuncforcexpressure&
&     )*liftdirection(1) + funcvalues(costfuncforceypressure)*&
&     liftdirection(2) + funcvalues(costfuncforcezpressure)*&
&     liftdirection(3)
    funcvalues(costfuncliftviscous) = funcvalues(costfuncforcexviscous)*&
&     liftdirection(1) + funcvalues(costfuncforceyviscous)*liftdirection&
&     (2) + funcvalues(costfuncforcezviscous)*liftdirection(3)
    funcvalues(costfuncliftmomentum) = funcvalues(costfuncforcexmomentum&
&     )*liftdirection(1) + funcvalues(costfuncforceymomentum)*&
&     liftdirection(2) + funcvalues(costfuncforcezmomentum)*&
&     liftdirection(3)
!-----
    funcvalues(costfuncdrag) = funcvalues(costfuncforcex)*dragdirection(&
&     1) + funcvalues(costfuncforcey)*dragdirection(2) + funcvalues(&
&     costfuncforcez)*dragdirection(3)
    funcvalues(costfuncdragpressure) = funcvalues(costfuncforcexpressure&
&     )*dragdirection(1) + funcvalues(costfuncforceypressure)*&
&     dragdirection(2) + funcvalues(costfuncforcezpressure)*&
&     dragdirection(3)
    funcvalues(costfuncdragviscous) = funcvalues(costfuncforcexviscous)*&
&     dragdirection(1) + funcvalues(costfuncforceyviscous)*dragdirection&
&     (2) + funcvalues(costfuncforcezviscous)*dragdirection(3)
    funcvalues(costfuncdragmomentum) = funcvalues(costfuncforcexmomentum&
&     )*dragdirection(1) + funcvalues(costfuncforceymomentum)*&
&     dragdirection(2) + funcvalues(costfuncforcezmomentum)*&
&     dragdirection(3)
!-----
    funcvalues(costfuncliftcoef) = funcvalues(costfuncforcexcoef)*&
&     liftdirection(1) + funcvalues(costfuncforceycoef)*liftdirection(2)&
&     + funcvalues(costfuncforcezcoef)*liftdirection(3)
    funcvalues(costfuncliftcoefpressure) = funcvalues(&
&     costfuncforcexcoefpressure)*liftdirection(1) + funcvalues(&
&     costfuncforceycoefpressure)*liftdirection(2) + funcvalues(&
&     costfuncforcezcoefpressure)*liftdirection(3)
    funcvalues(costfuncliftcoefviscous) = funcvalues(&
&     costfuncforcexcoefviscous)*liftdirection(1) + funcvalues(&
&     costfuncforceycoefviscous)*liftdirection(2) + funcvalues(&
&     costfuncforcezcoefviscous)*liftdirection(3)
    funcvalues(costfuncliftcoefmomentum) = funcvalues(&
&     costfuncforcexcoefmomentum)*liftdirection(1) + funcvalues(&
&     costfuncforceycoefmomentum)*liftdirection(2) + funcvalues(&
&     costfuncforcezcoefmomentum)*liftdirection(3)
!-----
    funcvalues(costfuncdragcoef) = funcvalues(costfuncforcexcoef)*&
&     dragdirection(1) + funcvalues(costfuncforceycoef)*dragdirection(2)&
&     + funcvalues(costfuncforcezcoef)*dragdirection(3)
    funcvalues(costfuncdragcoefpressure) = funcvalues(&
&     costfuncforcexcoefpressure)*dragdirection(1) + funcvalues(&
&     costfuncforceycoefpressure)*dragdirection(2) + funcvalues(&
&     costfuncforcezcoefpressure)*dragdirection(3)
    funcvalues(costfuncdragcoefviscous) = funcvalues(&
&     costfuncforcexcoefviscous)*dragdirection(1) + funcvalues(&
&     costfuncforceycoefviscous)*dragdirection(2) + funcvalues(&
&     costfuncforcezcoefviscous)*dragdirection(3)
    funcvalues(costfuncdragcoefmomentum) = funcvalues(&
&     costfuncforcexcoefmomentum)*dragdirection(1) + funcvalues(&
&     costfuncforceycoefmomentum)*dragdirection(2) + funcvalues(&
&     costfuncforcezcoefmomentum)*dragdirection(3)
! ----- center of lift
! dot product the 3 forces with the lift direction separately
    fxlift = funcvalues(costfuncforcex)*liftdirection(1)
    fylift = funcvalues(costfuncforcey)*liftdirection(2)
    fzlift = funcvalues(costfuncforcez)*liftdirection(3)
! run the weighed average for the 3 components of center of lift
! protect against division by zero
    if (fxlift + fylift + fzlift .ne. zero) then
      funcvalues(costfunccofliftx) = (fxlift*funcvalues(&
&       costfunccoforcexx)+fylift*funcvalues(costfunccoforceyx)+fzlift*&
&       funcvalues(costfunccoforcezx))/(fxlift+fylift+fzlift)
      funcvalues(costfunccoflifty) = (fxlift*funcvalues(&
&       costfunccoforcexy)+fylift*funcvalues(costfunccoforceyy)+fzlift*&
&       funcvalues(costfunccoforcezy))/(fxlift+fylift+fzlift)
      funcvalues(costfunccofliftz) = (fxlift*funcvalues(&
&       costfunccoforcexz)+fylift*funcvalues(costfunccoforceyz)+fzlift*&
&       funcvalues(costfunccoforcezz))/(fxlift+fylift+fzlift)
    else
      funcvalues(costfunccofliftx) = zero
      funcvalues(costfunccoflifty) = zero
      funcvalues(costfunccofliftz) = zero
    end if
! -------------------- time spectral objectives ------------------
    if (tsstability) then
      print*, &
&     'error: tsstabilityderivatives are *broken*. they need to be ', &
&     'completely verifed from scratch'
      stop
    end if
  end subroutine getcostfunctions
 
  subroutine wallintegrationface(localvalues, mm)
!
!       wallintegrations computes the contribution of the block
!       given by the pointers in blockpointers to the force and
!       moment of the geometry. a distinction is made
!       between the inviscid and viscous parts. in case the maximum
!       yplus value must be monitored (only possible for rans), this
!       value is also computed. the separation sensor and the cavita-
!       tion sensor is also computed
!       here.
!
    use constants
    use communication
    use blockpointers
    use flowvarrefstate
    use inputcostfunctions
    use inputphysics, only : machcoef, pointref, veldirfreestream, &
&   equations, momentaxis, cpmin_family, cpmin_rho, cavitationnumber, &
&   sepsenmaxfamily, sepsenmaxrho
    use bcpointers
    implicit none
! input/output variables
    real(kind=realtype), dimension(nlocalvalues), intent(inout) :: &
&   localvalues
    integer(kind=inttype) :: mm
! local variables.
    real(kind=realtype), dimension(3) :: fp, fv, mp, mv
    real(kind=realtype), dimension(3) :: cofsumfx, cofsumfy, cofsumfz
    real(kind=realtype) :: yplusmax, sepsensorks, sepsensor, &
&   sepsensoravg(3), sepsensorarea, cavitation, cpmin_ks_sum, &
&   sepsensorksarea
    integer(kind=inttype) :: i, j, ii, blk
    real(kind=realtype) :: pm1, fx, fy, fz, fn
    real(kind=realtype) :: vecttangential(3)
    real(kind=realtype) :: vectdotproductfsnormal
    real(kind=realtype) :: xc, xco, yc, yco, zc, zco, qf(3), r(3), n(3)&
&   , l
    real(kind=realtype) :: fact, rho, mul, yplus, dwall
    real(kind=realtype) :: v(3), sensor, sensor1, cp, tmp, plocal, &
&   ks_exponent
    real(kind=realtype) :: tauxx, tauyy, tauzz
    real(kind=realtype) :: tauxy, tauxz, tauyz
    real(kind=realtype), dimension(3) :: refpoint
    real(kind=realtype), dimension(3, 2) :: axispoints
    real(kind=realtype) :: mx, my, mz, cellarea, m0x, m0y, m0z, mvaxis, &
&   mpaxis
    real(kind=realtype) :: cperror, cperror2
    intrinsic mod
    intrinsic max
    intrinsic sqrt
    intrinsic cos
    intrinsic exp
    select case  (bcfaceid(mm)) 
    case (imin, jmin, kmin) 
      fact = -one
    case (imax, jmax, kmax) 
      fact = one
    end select
! determine the reference point for the moment computation in
! meters.
    refpoint(1) = lref*pointref(1)
    refpoint(2) = lref*pointref(2)
    refpoint(3) = lref*pointref(3)
! determine the points defining the axis about which to compute a moment
    axispoints(1, 1) = lref*momentaxis(1, 1)
    axispoints(1, 2) = lref*momentaxis(1, 2)
    axispoints(2, 1) = lref*momentaxis(2, 1)
    axispoints(2, 2) = lref*momentaxis(2, 2)
    axispoints(3, 1) = lref*momentaxis(3, 1)
    axispoints(3, 2) = lref*momentaxis(3, 2)
! initialize the force and moment coefficients to 0 as well as
! yplusmax.
    fp = zero
    fv = zero
    mp = zero
    mv = zero
    cofsumfx = zero
    cofsumfy = zero
    cofsumfz = zero
    yplusmax = zero
    sepsensor = zero
    sepsensorks = zero
    sepsensorarea = zero
    sepsensorksarea = zero
    cavitation = zero
    cpmin_ks_sum = zero
    sepsensoravg = zero
    mpaxis = zero
    mvaxis = zero
    cperror2 = zero
!$ad ii-loop
!
!         integrate the inviscid contribution over the solid walls,
!         either inviscid or viscous. the integration is done with
!         cp. for closed contours this is equal to the integration
!         of p; for open contours this is not the case anymore.
!         question is whether a force for an open contour is
!         meaningful anyway.
!
! loop over the quadrilateral faces of the subface. note that
! the nodal range of bcdata must be used and not the cell
! range, because the latter may include the halo's in i and
! j-direction. the offset +1 is there, because inbeg and jnbeg
! refer to nodal ranges and not to cell ranges. the loop
! (without the ad stuff) would look like:
!
! do j=(bcdata(mm)%jnbeg+1),bcdata(mm)%jnend
!    do i=(bcdata(mm)%inbeg+1),bcdata(mm)%inend
    do ii=0,(bcdata(mm)%jnend-bcdata(mm)%jnbeg)*(bcdata(mm)%inend-bcdata&
&       (mm)%inbeg)-1
      i = mod(ii, bcdata(mm)%inend - bcdata(mm)%inbeg) + bcdata(mm)%&
&       inbeg + 1
      j = ii/(bcdata(mm)%inend-bcdata(mm)%inbeg) + bcdata(mm)%jnbeg + 1
! compute the average pressure minus 1 and the coordinates
! of the centroid of the face relative from from the
! moment reference point. due to the usage of pointers for
! the coordinates, whose original array starts at 0, an
! offset of 1 must be used. the pressure is multipled by
! fact to account for the possibility of an inward or
! outward pointing normal.
      pm1 = fact*(half*(pp2(i, j)+pp1(i, j))-pinf)*pref
      tmp = two/(gammainf*pinf*machcoef*machcoef)
      cp = (half*(pp2(i, j)+pp1(i, j))-pinf)*tmp
      cperror = cp - bcdata(mm)%cptarget(i, j)
      cperror2 = cperror2 + cperror*cperror
      xc = fourth*(xx(i, j, 1)+xx(i+1, j, 1)+xx(i, j+1, 1)+xx(i+1, j+1, &
&       1)) - refpoint(1)
      yc = fourth*(xx(i, j, 2)+xx(i+1, j, 2)+xx(i, j+1, 2)+xx(i+1, j+1, &
&       2)) - refpoint(2)
      zc = fourth*(xx(i, j, 3)+xx(i+1, j, 3)+xx(i, j+1, 3)+xx(i+1, j+1, &
&       3)) - refpoint(3)
      if (bcdata(mm)%iblank(i, j) .lt. 0) then
        blk = 0
      else
        blk = bcdata(mm)%iblank(i, j)
      end if
      fx = pm1*ssi(i, j, 1)
      fy = pm1*ssi(i, j, 2)
      fz = pm1*ssi(i, j, 3)
! note from ay: technically, we can just compute the moments using the center of force
! terms. however, the moment computations coded here distinguish pressure,
! viscous, and momentum contributions to moment. even though these individual
! contributions are not exposed to python, i still wanted to keep how it's done in the
! code in case its useful in the future. this is also true for the face integrations
! update the inviscid force and moment coefficients. iblank as we sum
      fp(1) = fp(1) + fx*blk
      fp(2) = fp(2) + fy*blk
      fp(3) = fp(3) + fz*blk
      mx = yc*fz - zc*fy
      my = zc*fx - xc*fz
      mz = xc*fy - yc*fx
      mp(1) = mp(1) + mx*blk
      mp(2) = mp(2) + my*blk
      mp(3) = mp(3) + mz*blk
! the force integral for the center of pressure computation.
! we need the cell centers wrt origin
      xco = fourth*(xx(i, j, 1)+xx(i+1, j, 1)+xx(i, j+1, 1)+xx(i+1, j+1&
&       , 1))
      yco = fourth*(xx(i, j, 2)+xx(i+1, j, 2)+xx(i, j+1, 2)+xx(i+1, j+1&
&       , 2))
      zco = fourth*(xx(i, j, 3)+xx(i+1, j, 3)+xx(i, j+1, 3)+xx(i+1, j+1&
&       , 3))
! accumulate in the sums. each force component is tracked separately
! force-x
      cofsumfx(1) = cofsumfx(1) + xco*fx*blk
      cofsumfx(2) = cofsumfx(2) + yco*fx*blk
      cofsumfx(3) = cofsumfx(3) + zco*fx*blk
! force-y
      cofsumfy(1) = cofsumfy(1) + xco*fy*blk
      cofsumfy(2) = cofsumfy(2) + yco*fy*blk
      cofsumfy(3) = cofsumfy(3) + zco*fy*blk
! force-z
      cofsumfz(1) = cofsumfz(1) + xco*fz*blk
      cofsumfz(2) = cofsumfz(2) + yco*fz*blk
      cofsumfz(3) = cofsumfz(3) + zco*fz*blk
! compute the r and n vectores for the moment around an
! axis computation where r is the distance from the
! force to the first point on the axis and n is a unit
! normal in the direction of the axis
      r(1) = fourth*(xx(i, j, 1)+xx(i+1, j, 1)+xx(i, j+1, 1)+xx(i+1, j+1&
&       , 1)) - axispoints(1, 1)
      r(2) = fourth*(xx(i, j, 2)+xx(i+1, j, 2)+xx(i, j+1, 2)+xx(i+1, j+1&
&       , 2)) - axispoints(2, 1)
      r(3) = fourth*(xx(i, j, 3)+xx(i+1, j, 3)+xx(i, j+1, 3)+xx(i+1, j+1&
&       , 3)) - axispoints(3, 1)
      l = sqrt((axispoints(1, 2)-axispoints(1, 1))**2 + (axispoints(2, 2&
&       )-axispoints(2, 1))**2 + (axispoints(3, 2)-axispoints(3, 1))**2)
      n(1) = (axispoints(1, 2)-axispoints(1, 1))/l
      n(2) = (axispoints(2, 2)-axispoints(2, 1))/l
      n(3) = (axispoints(3, 2)-axispoints(3, 1))/l
! compute the moment of the force about the first point
! used to define the axis, and the project that axis in
! the n direction
      m0x = r(2)*fz - r(3)*fy
      m0y = r(3)*fx - r(1)*fz
      m0z = r(1)*fy - r(2)*fx
      mpaxis = mpaxis + (m0x*n(1)+m0y*n(2)+m0z*n(3))*blk
! save the face-based forces and area
      bcdata(mm)%fp(i, j, 1) = fx
      bcdata(mm)%fp(i, j, 2) = fy
      bcdata(mm)%fp(i, j, 3) = fz
      cellarea = sqrt(ssi(i, j, 1)**2 + ssi(i, j, 2)**2 + ssi(i, j, 3)**&
&       2)
      bcdata(mm)%area(i, j) = cellarea
! get normalized surface velocity:
      v(1) = ww2(i, j, ivx)
      v(2) = ww2(i, j, ivy)
      v(3) = ww2(i, j, ivz)
      v = v/(sqrt(v(1)**2+v(2)**2+v(3)**2)+1e-16)
      if (computesepsensorks) then
! freestream projection over the surface.
        vectdotproductfsnormal = veldirfreestream(1)*bcdata(mm)%norm(i, &
&         j, 1) + veldirfreestream(2)*bcdata(mm)%norm(i, j, 2) + &
&         veldirfreestream(3)*bcdata(mm)%norm(i, j, 3)
! tangential vector on the surface, which is the freestream projected vector
        vecttangential(1) = veldirfreestream(1) - vectdotproductfsnormal&
&         *bcdata(mm)%norm(i, j, 1)
        vecttangential(2) = veldirfreestream(2) - vectdotproductfsnormal&
&         *bcdata(mm)%norm(i, j, 2)
        vecttangential(3) = veldirfreestream(3) - vectdotproductfsnormal&
&         *bcdata(mm)%norm(i, j, 3)
        vecttangential = vecttangential/(sqrt(vecttangential(1)**2+&
&         vecttangential(2)**2+vecttangential(3)**2)+1e-16)
! computing separation sensor
! velocity dot products
        sensor = v(1)*vecttangential(1) + v(2)*vecttangential(2) + v(3)*&
&         vecttangential(3)
! sepsensor value
        sensor = (cos(degtorad*sepsensorksphi)-sensor)/(-cos(degtorad*&
&         sepsensorksphi)+cos(zero)+1e-16)
! also do the ks-based spensenor max computation
        call ksaggregationfunction(sensor, sepsenmaxfamily(spectralsol)&
&                            , sepsenmaxrho, ks_exponent)
        sepsensorksarea = sepsensorksarea + blk*cellarea
        sepsensorarea = cellarea*blk*one/(one+exp(-(2*&
&         sepsensorkssharpness*(sensor+sepsensorksoffset)))) + &
&         sepsensorarea
        sepsensorks = sepsensorks + ks_exponent*blk
      end if
! dot product with free stream
      sensor = -(v(1)*veldirfreestream(1)+v(2)*veldirfreestream(2)+v(3)*&
&       veldirfreestream(3))
!now run through a smooth heaviside function:
      sensor = one/(one+exp(-(2*sepsensorsharpness*(sensor-&
&       sepsensoroffset))))
! and integrate over the area of this cell and save, blanking as we go.
      sensor = sensor*cellarea*blk
      sepsensor = sepsensor + sensor
! also accumulate into the sepsensoravg
! x-y-zco are computed above for center of force computations
      sepsensoravg(1) = sepsensoravg(1) + sensor*xco
      sepsensoravg(2) = sepsensoravg(2) + sensor*yco
      sepsensoravg(3) = sepsensoravg(3) + sensor*zco
      if (computecavitation) then
        plocal = pp2(i, j)
        tmp = two/(gammainf*machcoef*machcoef)
        cp = tmp*(plocal-pinf)
        sensor1 = -cp - cavitationnumber
        sensor1 = sensor1**cavexponent/(one+exp(2*cavsensorsharpness*(-&
&         sensor1+cavsensoroffset)))
        sensor1 = sensor1*cellarea*blk
        cavitation = cavitation + sensor1
! also do the ks-based cpmin computation
        ks_exponent = exp(cpmin_rho*(-cp+cpmin_family(spectralsol)))
        cpmin_ks_sum = cpmin_ks_sum + ks_exponent*blk
      end if
    end do
!
! integration of the viscous forces.
! only for viscous boundaries.
!
    if (bctype(mm) .eq. nswalladiabatic .or. bctype(mm) .eq. &
&       nswallisothermal) then
! initialize dwall for the laminar case and set the pointer
! for the unit normals.
      dwall = zero
!$ad ii-loop
! loop over the quadrilateral faces of the subface and
! compute the viscous contribution to the force and
! moment and update the maximum value of y+.
      do ii=0,(bcdata(mm)%jnend-bcdata(mm)%jnbeg)*(bcdata(mm)%inend-&
&         bcdata(mm)%inbeg)-1
        i = mod(ii, bcdata(mm)%inend - bcdata(mm)%inbeg) + bcdata(mm)%&
&         inbeg + 1
        j = ii/(bcdata(mm)%inend-bcdata(mm)%inbeg) + bcdata(mm)%jnbeg + &
&         1
        if (bcdata(mm)%iblank(i, j) .lt. 0) then
          blk = 0
        else
          blk = bcdata(mm)%iblank(i, j)
        end if
        tauxx = viscsubface(mm)%tau(i, j, 1)
        tauyy = viscsubface(mm)%tau(i, j, 2)
        tauzz = viscsubface(mm)%tau(i, j, 3)
        tauxy = viscsubface(mm)%tau(i, j, 4)
        tauxz = viscsubface(mm)%tau(i, j, 5)
        tauyz = viscsubface(mm)%tau(i, j, 6)
! compute the viscous force on the face. a minus sign
! is now present, due to the definition of this force.
        fx = -(fact*(tauxx*ssi(i, j, 1)+tauxy*ssi(i, j, 2)+tauxz*ssi(i, &
&         j, 3))*pref)
        fy = -(fact*(tauxy*ssi(i, j, 1)+tauyy*ssi(i, j, 2)+tauyz*ssi(i, &
&         j, 3))*pref)
        fz = -(fact*(tauxz*ssi(i, j, 1)+tauyz*ssi(i, j, 2)+tauzz*ssi(i, &
&         j, 3))*pref)
! compute the coordinates of the centroid of the face
! relative from the moment reference point. due to the
! usage of pointers for xx and offset of 1 is present,
! because x originally starts at 0.
        xc = fourth*(xx(i, j, 1)+xx(i+1, j, 1)+xx(i, j+1, 1)+xx(i+1, j+1&
&         , 1)) - refpoint(1)
        yc = fourth*(xx(i, j, 2)+xx(i+1, j, 2)+xx(i, j+1, 2)+xx(i+1, j+1&
&         , 2)) - refpoint(2)
        zc = fourth*(xx(i, j, 3)+xx(i+1, j, 3)+xx(i, j+1, 3)+xx(i+1, j+1&
&         , 3)) - refpoint(3)
! update the viscous force and moment coefficients, blanking as we go.
        fv(1) = fv(1) + fx*blk
        fv(2) = fv(2) + fy*blk
        fv(3) = fv(3) + fz*blk
        mx = yc*fz - zc*fy
        my = zc*fx - xc*fz
        mz = xc*fy - yc*fx
        mv(1) = mv(1) + mx*blk
        mv(2) = mv(2) + my*blk
        mv(3) = mv(3) + mz*blk
! the force integral for the center of pressure computation.
! we need the cell centers wrt origin
        xco = fourth*(xx(i, j, 1)+xx(i+1, j, 1)+xx(i, j+1, 1)+xx(i+1, j+&
&         1, 1))
        yco = fourth*(xx(i, j, 2)+xx(i+1, j, 2)+xx(i, j+1, 2)+xx(i+1, j+&
&         1, 2))
        zco = fourth*(xx(i, j, 3)+xx(i+1, j, 3)+xx(i, j+1, 3)+xx(i+1, j+&
&         1, 3))
! accumulate in the sums. each force component is tracked separately
! force-x
        cofsumfx(1) = cofsumfx(1) + xco*fx*blk
        cofsumfx(2) = cofsumfx(2) + yco*fx*blk
        cofsumfx(3) = cofsumfx(3) + zco*fx*blk
! force-y
        cofsumfy(1) = cofsumfy(1) + xco*fy*blk
        cofsumfy(2) = cofsumfy(2) + yco*fy*blk
        cofsumfy(3) = cofsumfy(3) + zco*fy*blk
! force-z
        cofsumfz(1) = cofsumfz(1) + xco*fz*blk
        cofsumfz(2) = cofsumfz(2) + yco*fz*blk
        cofsumfz(3) = cofsumfz(3) + zco*fz*blk
! compute the r and n vectors for the moment around an
! axis computation where r is the distance from the
! force to the first point on the axis and n is a unit
! normal in the direction of the axis
        r(1) = fourth*(xx(i, j, 1)+xx(i+1, j, 1)+xx(i, j+1, 1)+xx(i+1, j&
&         +1, 1)) - axispoints(1, 1)
        r(2) = fourth*(xx(i, j, 2)+xx(i+1, j, 2)+xx(i, j+1, 2)+xx(i+1, j&
&         +1, 2)) - axispoints(2, 1)
        r(3) = fourth*(xx(i, j, 3)+xx(i+1, j, 3)+xx(i, j+1, 3)+xx(i+1, j&
&         +1, 3)) - axispoints(3, 1)
        l = sqrt((axispoints(1, 2)-axispoints(1, 1))**2 + (axispoints(2&
&         , 2)-axispoints(2, 1))**2 + (axispoints(3, 2)-axispoints(3, 1)&
&         )**2)
        n(1) = (axispoints(1, 2)-axispoints(1, 1))/l
        n(2) = (axispoints(2, 2)-axispoints(2, 1))/l
        n(3) = (axispoints(3, 2)-axispoints(3, 1))/l
! compute the moment of the force about the first point
! used to define the axis, and then project that axis in
! the n direction
        m0x = r(2)*fz - r(3)*fy
        m0y = r(3)*fx - r(1)*fz
        m0z = r(1)*fy - r(2)*fx
        mvaxis = mvaxis + (m0x*n(1)+m0y*n(2)+m0z*n(3))*blk
! save the face based forces for the slice operations
        bcdata(mm)%fv(i, j, 1) = fx
        bcdata(mm)%fv(i, j, 2) = fy
        bcdata(mm)%fv(i, j, 3) = fz
! compute the tangential component of the stress tensor,
! which is needed to monitor y+. the result is stored
! in fx, fy, fz, although it is not really a force.
! as later on only the magnitude of the tangential
! component is important, there is no need to take the
! sign into account (it should be a minus sign).
        fx = tauxx*bcdata(mm)%norm(i, j, 1) + tauxy*bcdata(mm)%norm(i, j&
&         , 2) + tauxz*bcdata(mm)%norm(i, j, 3)
        fy = tauxy*bcdata(mm)%norm(i, j, 1) + tauyy*bcdata(mm)%norm(i, j&
&         , 2) + tauyz*bcdata(mm)%norm(i, j, 3)
        fz = tauxz*bcdata(mm)%norm(i, j, 1) + tauyz*bcdata(mm)%norm(i, j&
&         , 2) + tauzz*bcdata(mm)%norm(i, j, 3)
        fn = fx*bcdata(mm)%norm(i, j, 1) + fy*bcdata(mm)%norm(i, j, 2) +&
&         fz*bcdata(mm)%norm(i, j, 3)
        fx = fx - fn*bcdata(mm)%norm(i, j, 1)
        fy = fy - fn*bcdata(mm)%norm(i, j, 2)
        fz = fz - fn*bcdata(mm)%norm(i, j, 3)
! compute the local value of y+. due to the usage
! of pointers there is on offset of -1 in dd2wall..
      end do
    else
! if we had no viscous force, set the viscous component to zero
      bcdata(mm)%fv = zero
    end if
! increment the local values array with the values we computed here.
    localvalues(ifp:ifp+2) = localvalues(ifp:ifp+2) + fp
    localvalues(ifv:ifv+2) = localvalues(ifv:ifv+2) + fv
    localvalues(imp:imp+2) = localvalues(imp:imp+2) + mp
    localvalues(imv:imv+2) = localvalues(imv:imv+2) + mv
    localvalues(icoforcex:icoforcex+2) = localvalues(icoforcex:icoforcex&
&     +2) + cofsumfx
    localvalues(icoforcey:icoforcey+2) = localvalues(icoforcey:icoforcey&
&     +2) + cofsumfy
    localvalues(icoforcez:icoforcez+2) = localvalues(icoforcez:icoforcez&
&     +2) + cofsumfz
    localvalues(isepsensor) = localvalues(isepsensor) + sepsensor
    localvalues(isepsensorks) = localvalues(isepsensorks) + sepsensorks
    localvalues(isepsensorksarea) = localvalues(isepsensorksarea) + &
&     sepsensorksarea
    localvalues(isepsensorarea) = localvalues(isepsensorarea) + &
&     sepsensorarea
    localvalues(icavitation) = localvalues(icavitation) + cavitation
    localvalues(icpmin) = localvalues(icpmin) + cpmin_ks_sum
    localvalues(isepavg:isepavg+2) = localvalues(isepavg:isepavg+2) + &
&     sepsensoravg
    localvalues(iaxismoment) = localvalues(iaxismoment) + mpaxis + &
&     mvaxis
    localvalues(icperror2) = localvalues(icperror2) + cperror2
  end subroutine wallintegrationface
 
  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
 
  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_fast_b, only : computeptot, computettot
    use bcpointers, only : ssi, sface, ww1, ww2, pp1, pp2, xx, gamma1, &
&   gamma2
    use utils_fast_b, only : mynorm2
    implicit none
! input/output variables
    logical, intent(in) :: isinflow
    real(kind=realtype), dimension(nlocalvalues), intent(inout) :: &
&   localvalues
    integer(kind=inttype), intent(in) :: mm
! local variables
    real(kind=realtype) :: massflowrate, mass_ptot, mass_ttot, mass_ps, &
&   mass_mn, mass_a, mass_rho, mass_vx, mass_vy, mass_vz, mass_nx, &
&   mass_ny, mass_nz, mass_vi
    real(kind=realtype) :: area_ptot, area_ps
    real(kind=realtype) :: govgm1, gm1ovg, viconst, vilocal, pratio
    real(kind=realtype) :: mredim
    integer(kind=inttype) :: i, j, ii, blk
    real(kind=realtype) :: internalflowfact, inflowfact, fact, xc, xco, &
&   yc, yco, zc, zco, mx, my, mz
    real(kind=realtype) :: sf, vmag, vnm, vnmfreestreamref, vxm, vym, &
&   vzm, fx, fy, fz, u, v, w
    real(kind=realtype) :: pm, ptot, ttot, rhom, gammam, am
    real(kind=realtype) :: area, cellarea, overcellarea
    real(kind=realtype), dimension(3) :: fp, mp, fmom, mmom, refpoint, &
&   sfacecoordref
    real(kind=realtype), dimension(3) :: cofsumfx, cofsumfy, cofsumfz
    real(kind=realtype) :: mnm, massflowratelocal
    intrinsic sqrt
    intrinsic mod
    intrinsic max
    intrinsic min
    refpoint(1) = lref*pointref(1)
    refpoint(2) = lref*pointref(2)
    refpoint(3) = lref*pointref(3)
! note that these are *opposite* of force integrations. the reason
! is that we want positive mass flow into the domain and negative
! mass flow out of the domain. since the low faces have ssi
! vectors pointining into the domain, this is correct. the high
! end faces need to flip this.
    select case  (bcfaceid(mm)) 
    case (imin, jmin, kmin) 
      fact = one
    case (imax, jmax, kmax) 
      fact = -one
    end select
! the sign of momentum forces are flipped for internal flows
    internalflowfact = one
    if (flowtype .eq. internalflow) internalflowfact = -one
    inflowfact = one
    if (isinflow) inflowfact = -one
! loop over the quadrilateral faces of the subface. note that
! the nodal range of bcdata must be used and not the cell
! range, because the latter may include the halo's in i and
! j-direction. the offset +1 is there, because inbeg and jnbeg
! refer to nodal ranges and not to cell ranges. the loop
! (without the ad stuff) would look like:
!
! do j=(bcdata(mm)%jnbeg+1),bcdata(mm)%jnend
!    do i=(bcdata(mm)%inbeg+1),bcdata(mm)%inend
    mredim = sqrt(pref*rhoref)
    fp = zero
    mp = zero
    fmom = zero
    mmom = zero
    cofsumfx = zero
    cofsumfy = zero
    cofsumfz = zero
    massflowrate = zero
    area = zero
    mass_ptot = zero
    mass_ttot = zero
    mass_ps = zero
    mass_mn = zero
    mass_a = zero
    mass_rho = zero
    mass_vx = zero
    mass_vy = zero
    mass_vz = zero
    mass_nx = zero
    mass_ny = zero
    mass_nz = zero
    mass_vi = zero
    area_ptot = zero
    area_ps = zero
!$ad ii-loop
    do ii=0,(bcdata(mm)%jnend-bcdata(mm)%jnbeg)*(bcdata(mm)%inend-bcdata&
&       (mm)%inbeg)-1
      i = mod(ii, bcdata(mm)%inend - bcdata(mm)%inbeg) + bcdata(mm)%&
&       inbeg + 1
      j = ii/(bcdata(mm)%inend-bcdata(mm)%inbeg) + bcdata(mm)%jnbeg + 1
      if (addgridvelocities) then
        sf = sface(i, j)
      else
        sf = zero
      end if
      if (bcdata(mm)%iblank(i, j) .lt. 0) then
        blk = 0
      else
        blk = bcdata(mm)%iblank(i, j)
      end if
      vxm = half*(ww1(i, j, ivx)+ww2(i, j, ivx))
      vym = half*(ww1(i, j, ivy)+ww2(i, j, ivy))
      vzm = half*(ww1(i, j, ivz)+ww2(i, j, ivz))
      rhom = half*(ww1(i, j, irho)+ww2(i, j, irho))
      pm = half*(pp1(i, j)+pp2(i, j))
      gammam = half*(gamma1(i, j)+gamma2(i, j))
      vnm = vxm*ssi(i, j, 1) + vym*ssi(i, j, 2) + vzm*ssi(i, j, 3) - sf
      vmag = sqrt(vxm**2 + vym**2 + vzm**2) - sf
      am = sqrt(gammam*pm/rhom)
      mnm = vmag/am
      cellarea = sqrt(ssi(i, j, 1)**2 + ssi(i, j, 2)**2 + ssi(i, j, 3)**&
&       2)
      area = area + cellarea*blk
      overcellarea = 1/cellarea
      call computeptot(rhom, vxm, vym, vzm, pm, ptot)
      call computettot(rhom, vxm, vym, vzm, pm, ttot)
      massflowratelocal = rhom*vnm*blk*fact*mredim
      massflowrate = massflowrate + massflowratelocal
! re-dimentionalize quantities
      pm = pm*pref
      mass_ptot = mass_ptot + ptot*massflowratelocal*pref
      mass_ttot = mass_ttot + ttot*massflowratelocal*tref
      mass_rho = mass_rho + rhom*massflowratelocal*rhoref
      mass_a = mass_a + am*massflowratelocal*uref
      mass_ps = mass_ps + pm*massflowratelocal
      mass_mn = mass_mn + mnm*massflowratelocal
      area_ptot = area_ptot + ptot*pref*cellarea*blk
      area_ps = area_ps + pm*cellarea*blk
      sfacecoordref(1) = sf*ssi(i, j, 1)*overcellarea
      sfacecoordref(2) = sf*ssi(i, j, 2)*overcellarea
      sfacecoordref(3) = sf*ssi(i, j, 3)*overcellarea
      mass_vx = mass_vx + (vxm*uref-sfacecoordref(1))*massflowratelocal
      mass_vy = mass_vy + (vym*uref-sfacecoordref(2))*massflowratelocal
      mass_vz = mass_vz + (vzm*uref-sfacecoordref(3))*massflowratelocal
      govgm1 = gammainf/(gammainf-one)
      gm1ovg = one/govgm1
      viconst = two*govgm1*rgasdim
      if (one .gt. one/ptot) then
        pratio = one/ptot
      else
        pratio = one
      end if
      vilocal = sqrt(viconst*(one-pratio**gm1ovg)*ttot*tref)
      mass_vi = mass_vi + vilocal*massflowratelocal
      mass_nx = mass_nx + ssi(i, j, 1)*overcellarea*massflowratelocal
      mass_ny = mass_ny + ssi(i, j, 2)*overcellarea*massflowratelocal
      mass_nz = mass_nz + ssi(i, j, 3)*overcellarea*massflowratelocal
      xc = fourth*(xx(i, j, 1)+xx(i+1, j, 1)+xx(i, j+1, 1)+xx(i+1, j+1, &
&       1)) - refpoint(1)
      yc = fourth*(xx(i, j, 2)+xx(i+1, j, 2)+xx(i, j+1, 2)+xx(i+1, j+1, &
&       2)) - refpoint(2)
      zc = fourth*(xx(i, j, 3)+xx(i+1, j, 3)+xx(i, j+1, 3)+xx(i+1, j+1, &
&       3)) - refpoint(3)
! pressure forces. note that these need a *negative* and to subtract
! the reference pressure sign to be consistent with the force
! computation on the walls.
      pm = -((pm-pinf*pref)*fact*blk)
      fx = pm*ssi(i, j, 1)
      fy = pm*ssi(i, j, 2)
      fz = pm*ssi(i, j, 3)
! update the pressure force and moment coefficients.
      fp(1) = fp(1) + fx
      fp(2) = fp(2) + fy
      fp(3) = fp(3) + fz
      mx = yc*fz - zc*fy
      my = zc*fx - xc*fz
      mz = xc*fy - yc*fx
      mp(1) = mp(1) + mx
      mp(2) = mp(2) + my
      mp(3) = mp(3) + mz
! the force integral for the center of pressure computation.
! we need the cell centers wrt origin
      xco = fourth*(xx(i, j, 1)+xx(i+1, j, 1)+xx(i, j+1, 1)+xx(i+1, j+1&
&       , 1))
      yco = fourth*(xx(i, j, 2)+xx(i+1, j, 2)+xx(i, j+1, 2)+xx(i+1, j+1&
&       , 2))
      zco = fourth*(xx(i, j, 3)+xx(i+1, j, 3)+xx(i, j+1, 3)+xx(i+1, j+1&
&       , 3))
! center of force computations. here we accumulate in the sums.
! accumulate in the sums. each force component is tracked separately
! blanking is included in the mdot multiplier for the force.
! force-x
      cofsumfx(1) = cofsumfx(1) + xco*fx
      cofsumfx(2) = cofsumfx(2) + yco*fx
      cofsumfx(3) = cofsumfx(3) + zco*fx
! force-y
      cofsumfy(1) = cofsumfy(1) + xco*fy
      cofsumfy(2) = cofsumfy(2) + yco*fy
      cofsumfy(3) = cofsumfy(3) + zco*fy
! force-z
      cofsumfz(1) = cofsumfz(1) + xco*fz
      cofsumfz(2) = cofsumfz(2) + yco*fz
      cofsumfz(3) = cofsumfz(3) + zco*fz
! momentum forces are a little tricky.  we negate because
! have to re-apply fact to massflowratelocal to undoo it, because
! we need the signed behavior of ssi to get the momentum forces correct.
! also, the sign is flipped between inflow and outflow types
      massflowratelocal = massflowratelocal*fact/timeref*blk/cellarea*&
&       internalflowfact*inflowfact
      fx = massflowratelocal*ssi(i, j, 1)*vxm
      fy = massflowratelocal*ssi(i, j, 2)*vym
      fz = massflowratelocal*ssi(i, j, 3)*vzm
      fmom(1) = fmom(1) + fx
      fmom(2) = fmom(2) + fy
      fmom(3) = fmom(3) + fz
      mx = yc*fz - zc*fy
      my = zc*fx - xc*fz
      mz = xc*fy - yc*fx
      mmom(1) = mmom(1) + mx
      mmom(2) = mmom(2) + my
      mmom(3) = mmom(3) + mz
! center of force computations. here we accumulate in the sums.
! each force component is tracked separately
! blanking is included in the mdot multiplier for the force.
! force-x
      cofsumfx(1) = cofsumfx(1) + xco*fx
      cofsumfx(2) = cofsumfx(2) + yco*fx
      cofsumfx(3) = cofsumfx(3) + zco*fx
! force-y
      cofsumfy(1) = cofsumfy(1) + xco*fy
      cofsumfy(2) = cofsumfy(2) + yco*fy
      cofsumfy(3) = cofsumfy(3) + zco*fy
! force-z
      cofsumfz(1) = cofsumfz(1) + xco*fz
      cofsumfz(2) = cofsumfz(2) + yco*fz
      cofsumfz(3) = cofsumfz(3) + zco*fz
    end do
! increment the local values array with what we computed here
    localvalues(imassflow) = localvalues(imassflow) + massflowrate
    localvalues(iarea) = localvalues(iarea) + area
    localvalues(imassrho) = localvalues(imassrho) + mass_rho
    localvalues(imassa) = localvalues(imassa) + mass_a
    localvalues(imassptot) = localvalues(imassptot) + mass_ptot
    localvalues(imassttot) = localvalues(imassttot) + mass_ttot
    localvalues(imassps) = localvalues(imassps) + mass_ps
    localvalues(imassmn) = localvalues(imassmn) + mass_mn
    localvalues(ifp:ifp+2) = localvalues(ifp:ifp+2) + fp
    localvalues(iflowfm:iflowfm+2) = localvalues(iflowfm:iflowfm+2) + &
&     fmom
    localvalues(iflowmp:iflowmp+2) = localvalues(iflowmp:iflowmp+2) + mp
    localvalues(iflowmm:iflowmm+2) = localvalues(iflowmm:iflowmm+2) + &
&     mmom
    localvalues(icoforcex:icoforcex+2) = localvalues(icoforcex:icoforcex&
&     +2) + cofsumfx
    localvalues(icoforcey:icoforcey+2) = localvalues(icoforcey:icoforcey&
&     +2) + cofsumfy
    localvalues(icoforcez:icoforcez+2) = localvalues(icoforcez:icoforcez&
&     +2) + cofsumfz
    localvalues(iareaptot) = localvalues(iareaptot) + area_ptot
    localvalues(iareaps) = localvalues(iareaps) + area_ps
    localvalues(imassvx) = localvalues(imassvx) + mass_vx
    localvalues(imassvy) = localvalues(imassvy) + mass_vy
    localvalues(imassvz) = localvalues(imassvz) + mass_vz
    localvalues(imassnx) = localvalues(imassnx) + mass_nx
    localvalues(imassny) = localvalues(imassny) + mass_ny
    localvalues(imassnz) = localvalues(imassnz) + mass_nz
    localvalues(imassvi) = localvalues(imassvi) + mass_vi
  end subroutine flowintegrationface
! ----------------------------------------------------------------------
!                                                                      |
!                    no tapenade routine below this line               |
!                                                                      |
! ----------------------------------------------------------------------
 
end module surfaceintegrations_fast_b