zipperIntegrations_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 zipperintegrations_b
  implicit none
 
contains
!  differentiation of flowintegrationzipper in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: timeref tref pref rgas rhoref
!                pointref vars localvalues
!   with respect to varying inputs: timeref tref pref rgas rhoref
!                pointref vars localvalues
!   rw status of diff variables: timeref:incr tref:incr pref:incr
!                rgas:incr rhoref:incr pointref:incr vars:incr
!                localvalues:in-out
  subroutine flowintegrationzipper_b(isinflow, conn, fams, vars, varsd, &
&   localvalues, localvaluesd, famlist, sps, ptvalid)
! integrate over the trianges for the inflow/outflow conditions.
    use constants
    use blockpointers, only : bctype
    use sorting, only : faminlist
    use flowvarrefstate, only : pref, prefd, pinf, pinfd, rhoref, &
&   rhorefd, pref, prefd, timeref, timerefd, lref, tref, trefd, rgas, &
&   rgasd, uref, urefd, uinf, uinfd, rhoinf, rhoinfd, gammainf
    use inputphysics, only : pointref, pointrefd, flowtype, rgasdim
    use flowutils_b, only : computeptot, computeptot_b, computettot, &
&   computettot_b
    use surfacefamilies, only : familyexchange, bcfamexchange
    use utils_b, only : mynorm2, cross_prod, cross_prod_b
    implicit none
! input/output variables
    logical, intent(in) :: isinflow
    integer(kind=inttype), dimension(:, :), intent(in) :: conn
    integer(kind=inttype), dimension(:), intent(in) :: fams
    real(kind=realtype), dimension(:, :), intent(in) :: vars
    real(kind=realtype), dimension(:, :) :: varsd
    real(kind=realtype), dimension(nlocalvalues), intent(inout) :: &
&   localvalues
    real(kind=realtype), dimension(nlocalvalues), intent(inout) :: &
&   localvaluesd
    integer(kind=inttype), dimension(:), intent(in) :: famlist
    integer(kind=inttype), intent(in) :: sps
    logical(kind=inttype), dimension(:), optional, intent(in) :: ptvalid
! working variables
    integer(kind=inttype) :: i, j
    real(kind=realtype) :: sf, vmag, vnm, vxm, vym, vzm, fx, fy, fz, u, &
&   v, w, vnmfreestreamref
    real(kind=realtype) :: sfd, vmagd, vnmd, vxmd, vymd, vzmd, fxd, fyd&
&   , fzd
    real(kind=realtype), dimension(3) :: fp, mp, fmom, mmom, refpoint, &
&   ss, x1, x2, x3, norm, sfacecoordref
    real(kind=realtype), dimension(3) :: fpd, mpd, fmomd, mmomd, &
&   refpointd, ssd, x1d, x2d, x3d, normd, sfacecoordrefd
    real(kind=realtype) :: pm, ptot, ttot, rhom, gammam, mnm, &
&   massflowratelocal, am
    real(kind=realtype) :: pmd, ptotd, ttotd, rhomd, gammamd, mnmd, &
&   massflowratelocald, amd
    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, cellarea, overcellarea
    real(kind=realtype) :: aread, cellaread, overcellaread
    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
    real(kind=realtype) :: internalflowfact, inflowfact, xc, xco, yc, &
&   yco, zc, zco, mx, my, mz
    real(kind=realtype) :: xcd, xcod, ycd, ycod, zcd, zcod, mxd, myd, &
&   mzd
    real(kind=realtype), dimension(3) :: cofsumfx, cofsumfy, cofsumfz
    real(kind=realtype), dimension(3) :: cofsumfxd, cofsumfyd, cofsumfzd
    logical :: triisvalid
    intrinsic sqrt
    intrinsic size
    intrinsic present
    intrinsic min
    real(kind=realtype), dimension(3) :: arg1
    real(kind=realtype), dimension(3) :: arg1d
    real(kind=realtype), dimension(3) :: arg2
    real(kind=realtype), dimension(3) :: arg2d
    real(kind=realtype) :: temp
    real(kind=realtype) :: tempd
    real(kind=realtype) :: temp0
    real(kind=realtype) :: tempd0
    integer :: branch
    mredim = sqrt(pref*rhoref)
    refpoint(1) = lref*pointref(1)
    refpoint(2) = lref*pointref(2)
    refpoint(3) = lref*pointref(3)
    internalflowfact = one
    if (flowtype .eq. internalflow) internalflowfact = -one
    inflowfact = one
    if (isinflow) inflowfact = -one
    mass_vid = localvaluesd(imassvi)
    mass_nzd = localvaluesd(imassnz)
    mass_nyd = localvaluesd(imassny)
    mass_nxd = localvaluesd(imassnx)
    mass_vzd = localvaluesd(imassvz)
    mass_vyd = localvaluesd(imassvy)
    mass_vxd = localvaluesd(imassvx)
    area_psd = localvaluesd(iareaps)
    area_ptotd = localvaluesd(iareaptot)
    cofsumfzd = 0.0_8
    cofsumfzd = localvaluesd(icoforcez:icoforcez+2)
    cofsumfyd = 0.0_8
    cofsumfyd = localvaluesd(icoforcey:icoforcey+2)
    cofsumfxd = 0.0_8
    cofsumfxd = localvaluesd(icoforcex:icoforcex+2)
    mmomd = 0.0_8
    mmomd = localvaluesd(iflowmm:iflowmm+2)
    mpd = 0.0_8
    mpd = localvaluesd(iflowmp:iflowmp+2)
    fmomd = 0.0_8
    fmomd = localvaluesd(iflowfm:iflowfm+2)
    fpd = 0.0_8
    fpd = localvaluesd(ifp:ifp+2)
    mass_mnd = localvaluesd(imassmn)
    mass_psd = localvaluesd(imassps)
    mass_ttotd = localvaluesd(imassttot)
    mass_ptotd = localvaluesd(imassptot)
    mass_ad = localvaluesd(imassa)
    mass_rhod = localvaluesd(imassrho)
    aread = localvaluesd(iarea)
    massflowrated = localvaluesd(imassflow)
    mredimd = 0.0_8
    normd = 0.0_8
    ptotd = 0.0_8
    refpointd = 0.0_8
    sfacecoordrefd = 0.0_8
    ttotd = 0.0_8
!$bwd-of ii-loop 
    do i=1,size(conn, 2)
      if (faminlist(fams(i), famlist)) then
! if the ptvalid list is given, check if we should integrate
! this triangle.
        triisvalid = .true.
        if (present(ptvalid)) then
! check if each of the three nodes are valid
          if (((ptvalid(conn(1, i)) .eqv. .false.) .or. (ptvalid(conn(2&
&             , i)) .eqv. .false.)) .or. (ptvalid(conn(3, i)) .eqv. &
&             .false.)) then
            call pushcontrol1b(0)
            triisvalid = .false.
          else
            call pushcontrol1b(0)
          end if
        else
          call pushcontrol1b(1)
        end if
        if (triisvalid) then
! compute the averaged values for this triangle
          vxm = zero
          vym = zero
          vzm = zero
          rhom = zero
          pm = zero
          gammam = zero
          sf = zero
          do j=1,3
            rhom = rhom + vars(conn(j, i), irho)
            vxm = vxm + vars(conn(j, i), ivx)
            vym = vym + vars(conn(j, i), ivy)
            vzm = vzm + vars(conn(j, i), ivz)
            pm = pm + vars(conn(j, i), irhoe)
            gammam = gammam + vars(conn(j, i), izippflowgamma)
            sf = sf + vars(conn(j, i), izippflowsface)
          end do
! divide by 3 due to the summation above:
          rhom = third*rhom
          vxm = third*vxm
          vym = third*vym
          vzm = third*vzm
          pm = third*pm
          gammam = third*gammam
          sf = third*sf
! get the nodes of triangle.
          x1 = vars(conn(1, i), izippflowx:izippflowz)
          x2 = vars(conn(2, i), izippflowx:izippflowz)
          x3 = vars(conn(3, i), izippflowx:izippflowz)
          arg1(:) = x2 - x1
          arg2(:) = x3 - x1
          call cross_prod(arg1(:), arg2(:), norm)
          ss = half*norm
          call computeptot(rhom, vxm, vym, vzm, pm, ptot)
          call computettot(rhom, vxm, vym, vzm, pm, ttot)
          vnm = vxm*ss(1) + vym*ss(2) + vzm*ss(3) - sf
          vmag = sqrt(vxm**2 + vym**2 + vzm**2) - sf
          am = sqrt(gammam*pm/rhom)
          mnm = vmag/sqrt(gammam*pm/rhom)
          cellarea = sqrt(ss(1)**2 + ss(2)**2 + ss(3)**2)
          overcellarea = 1/cellarea
          massflowratelocal = rhom*vnm*mredim
          call pushreal8(pm)
          pm = pm*pref
          sfacecoordref(1) = sf*ss(1)*overcellarea
          sfacecoordref(2) = sf*ss(2)*overcellarea
          sfacecoordref(3) = sf*ss(3)*overcellarea
          govgm1 = gammainf/(gammainf-one)
          gm1ovg = one/govgm1
          viconst = two*govgm1*rgasdim
          if (one .gt. one/ptot) then
            pratio = one/ptot
            call pushcontrol1b(0)
          else
            call pushcontrol1b(1)
            pratio = one
          end if
          vilocal = sqrt(viconst*(one-pratio**gm1ovg)*ttot*tref)
! compute the average cell center.
          xco = zero
          yco = zero
          zco = zero
          do j=1,3
            xco = xco + vars(conn(1, i), izippflowx)
            yco = yco + vars(conn(2, i), izippflowy)
            zco = zco + vars(conn(3, i), izippflowz)
          end do
! finish average for cell center
          xco = third*xco
          yco = third*yco
          zco = third*zco
! x-y-zco is the cell center w.r.t. the origin, x-y-zc is w.r.t. the reference point
          xc = xco - refpoint(1)
          yc = yco - refpoint(2)
          zc = zco - refpoint(3)
          call pushreal8(pm)
          pm = -(pm-pinf*pref)
! update the pressure force and moment coefficients.
! center of force computations. here we accumulate in the sums.
! force-x
! force-y
! force-z
! momentum forces
! get unit normal vector.
          call pushreal8array(ss, 3)
          ss = ss/cellarea
          call pushreal8(massflowratelocal)
          massflowratelocal = massflowratelocal/timeref*internalflowfact&
&           *inflowfact
          fx = massflowratelocal*ss(1)*vxm
          fy = massflowratelocal*ss(2)*vym
          fz = massflowratelocal*ss(3)*vzm
! center of force computations. here we accumulate in the sums.
! force-x
! force-y
! force-z
          mxd = mmomd(1)
          myd = mmomd(2)
          mzd = mmomd(3)
          zcod = fz*cofsumfzd(3) + fy*cofsumfyd(3) + fx*cofsumfxd(3)
          fzd = zco*cofsumfzd(3) + yco*cofsumfzd(2) + xco*cofsumfzd(1) +&
&           yc*mxd - xc*myd - fmomd(3)
          ycod = fz*cofsumfzd(2) + fy*cofsumfyd(2) + fx*cofsumfxd(2)
          xcod = fz*cofsumfzd(1) + fy*cofsumfyd(1) + fx*cofsumfxd(1)
          fyd = zco*cofsumfyd(3) + yco*cofsumfyd(2) + xco*cofsumfyd(1) +&
&           xc*mzd - zc*mxd - fmomd(2)
          fxd = zco*cofsumfxd(3) + yco*cofsumfxd(2) + xco*cofsumfxd(1) +&
&           zc*myd - yc*mzd - fmomd(1)
          xcd = fy*mzd - fz*myd
          ycd = fz*mxd - fx*mzd
          zcd = fx*myd - fy*mxd
          ssd = 0.0_8
          ssd(3) = ssd(3) + massflowratelocal*vzm*fzd
          massflowratelocald = vzm*ss(3)*fzd + vym*ss(2)*fyd + vxm*ss(1)&
&           *fxd
          vzmd = massflowratelocal*ss(3)*fzd
          ssd(2) = ssd(2) + massflowratelocal*vym*fyd
          vymd = massflowratelocal*ss(2)*fyd
          ssd(1) = ssd(1) + massflowratelocal*vxm*fxd
          vxmd = massflowratelocal*ss(1)*fxd
          call popreal8(massflowratelocal)
          tempd0 = internalflowfact*inflowfact*massflowratelocald/&
&           timeref
          massflowratelocald = tempd0
          timerefd = timerefd - massflowratelocal*tempd0/timeref
          call popreal8array(ss, 3)
          cellaread = -(sum(ss*ssd)/cellarea**2)
          ssd = ssd/cellarea
          fz = pm*ss(3)
          fy = pm*ss(2)
          fx = pm*ss(1)
          mzd = mpd(3)
          myd = mpd(2)
          fxd = zco*cofsumfxd(3) + yco*cofsumfxd(2) + xco*cofsumfxd(1) +&
&           zc*myd - yc*mzd + fpd(1)
          mxd = mpd(1)
          fzd = zco*cofsumfzd(3) + yco*cofsumfzd(2) + xco*cofsumfzd(1) +&
&           yc*mxd - xc*myd + fpd(3)
          fyd = zco*cofsumfyd(3) + yco*cofsumfyd(2) + xco*cofsumfyd(1) +&
&           xc*mzd + fpd(2) - zc*mxd
          xcd = xcd + fy*mzd - fz*myd
          xcod = xcod + fz*cofsumfzd(1) + fy*cofsumfyd(1) + fx*cofsumfxd&
&           (1) + xcd
          ycd = ycd + fz*mxd - fx*mzd
          ycod = ycod + fz*cofsumfzd(2) + fy*cofsumfyd(2) + fx*cofsumfxd&
&           (2) + ycd
          zcd = zcd + fx*myd - fy*mxd
          zcod = zcod + fz*cofsumfzd(3) + fy*cofsumfyd(3) + fx*cofsumfxd&
&           (3) + zcd
          pmd = ss(3)*fzd + ss(2)*fyd + ss(1)*fxd
          ssd(3) = ssd(3) + pm*fzd
          ssd(2) = ssd(2) + pm*fyd
          ssd(1) = ssd(1) + pm*fxd
          call popreal8(pm)
          prefd = prefd + pinf*pmd
          pmd = -pmd
          refpointd(3) = refpointd(3) - zcd
          refpointd(2) = refpointd(2) - ycd
          refpointd(1) = refpointd(1) - xcd
          zcod = third*zcod
          ycod = third*ycod
          xcod = third*xcod
          do j=3,1,-1
            varsd(conn(3, i), izippflowz) = varsd(conn(3, i), izippflowz&
&             ) + zcod
            varsd(conn(2, i), izippflowy) = varsd(conn(2, i), izippflowy&
&             ) + ycod
            varsd(conn(1, i), izippflowx) = varsd(conn(1, i), izippflowx&
&             ) + xcod
          end do
          ssd(3) = ssd(3) + overcellarea*massflowratelocal*mass_nzd
          overcellaread = massflowratelocal*ss(3)*mass_nzd + &
&           massflowratelocal*ss(2)*mass_nyd + massflowratelocal*ss(1)*&
&           mass_nxd
          massflowratelocald = massflowratelocald + overcellarea*ss(3)*&
&           mass_nzd + overcellarea*ss(2)*mass_nyd + overcellarea*ss(1)*&
&           mass_nxd + vilocal*mass_vid
          ssd(2) = ssd(2) + overcellarea*massflowratelocal*mass_nyd
          ssd(1) = ssd(1) + overcellarea*massflowratelocal*mass_nxd
          vilocald = massflowratelocal*mass_vid
          temp0 = one - pratio**gm1ovg
          if (viconst*(temp0*(ttot*tref)) .eq. 0.0_8) then
            tempd = 0.0_8
          else
            tempd = viconst*vilocald/(2.0*sqrt(viconst*(temp0*(ttot*tref&
&             ))))
          end if
          if (pratio .le. 0.0_8 .and. (gm1ovg .eq. 0.0_8 .or. gm1ovg &
&             .ne. int(gm1ovg))) then
            pratiod = 0.0_8
          else
            pratiod = -(gm1ovg*pratio**(gm1ovg-1)*ttot*tref*tempd)
          end if
          ttotd = ttotd + tref*temp0*tempd
          trefd = trefd + ttot*temp0*tempd
          call popcontrol1b(branch)
          if (branch .eq. 0) ptotd = ptotd - one*pratiod/ptot**2
          sfacecoordrefd(3) = sfacecoordrefd(3) - massflowratelocal*&
&           mass_vzd
          massflowratelocald = massflowratelocald + (uref*vzm-&
&           sfacecoordref(3))*mass_vzd + (uref*vym-sfacecoordref(2))*&
&           mass_vyd + (uref*vxm-sfacecoordref(1))*mass_vxd + mnm*&
&           mass_mnd + pm*mass_psd + am*uref*mass_ad + rhom*rhoref*&
&           mass_rhod + ttot*tref*mass_ttotd + ptot*pref*mass_ptotd + &
&           massflowrated
          sfacecoordrefd(2) = sfacecoordrefd(2) - massflowratelocal*&
&           mass_vyd
          sfacecoordrefd(1) = sfacecoordrefd(1) - massflowratelocal*&
&           mass_vxd
          overcellaread = overcellaread + sf*ss(3)*sfacecoordrefd(3) + &
&           sf*ss(2)*sfacecoordrefd(2) + sf*ss(1)*sfacecoordrefd(1)
          pmd = pmd + cellarea*area_psd + massflowratelocal*mass_psd
          cellaread = cellaread + pm*area_psd + ptot*pref*area_ptotd + &
&           aread - overcellaread/cellarea**2
          mnmd = massflowratelocal*mass_mnd
          amd = massflowratelocal*uref*mass_ad
          call popreal8(pm)
          prefd = prefd + ptot*cellarea*area_ptotd + ptot*&
&           massflowratelocal*mass_ptotd + pm*pmd
          vnmd = rhom*mredim*massflowratelocald
          if (ss(1)**2 + ss(2)**2 + ss(3)**2 .eq. 0.0_8) then
            tempd0 = 0.0_8
          else
            tempd0 = cellaread/(2.0*sqrt(ss(1)**2+ss(2)**2+ss(3)**2))
          end if
          temp = gammam*pm/rhom
          temp0 = sqrt(temp)
          vmagd = mnmd/temp0
          if (temp .eq. 0.0_8) then
            tempd = 0.0_8
          else
            tempd = -(vmag*mnmd/(rhom*2.0*temp0**3))
          end if
          rhomd = massflowratelocal*rhoref*mass_rhod + vnm*mredim*&
&           massflowratelocald - temp*tempd
          pmd = pref*pmd + gammam*tempd
          gammamd = pm*tempd
          temp = gammam*pm/rhom
          if (temp .eq. 0.0_8) then
            tempd = 0.0_8
          else
            tempd = amd/(rhom*2.0*sqrt(temp))
          end if
          gammamd = gammamd + pm*tempd
          pmd = pmd + gammam*tempd
          rhomd = rhomd - temp*tempd
          if (vxm**2 + vym**2 + vzm**2 .eq. 0.0_8) then
            tempd = 0.0_8
          else
            tempd = vmagd/(2.0*sqrt(vxm**2+vym**2+vzm**2))
          end if
          vzmd = vzmd + uref*massflowratelocal*mass_vzd + 2*vzm*tempd + &
&           ss(3)*vnmd
          vymd = vymd + uref*massflowratelocal*mass_vyd + 2*vym*tempd + &
&           ss(2)*vnmd
          vxmd = vxmd + uref*massflowratelocal*mass_vxd + 2*vxm*tempd + &
&           ss(1)*vnmd
          ssd(3) = ssd(3) + sf*overcellarea*sfacecoordrefd(3) + 2*ss(3)*&
&           tempd0 + vzm*vnmd
          sfd = overcellarea*ss(3)*sfacecoordrefd(3) + overcellarea*ss(2&
&           )*sfacecoordrefd(2) + overcellarea*ss(1)*sfacecoordrefd(1) -&
&           vmagd - vnmd
          sfacecoordrefd(3) = 0.0_8
          ssd(2) = ssd(2) + sf*overcellarea*sfacecoordrefd(2) + 2*ss(2)*&
&           tempd0 + vym*vnmd
          sfacecoordrefd(2) = 0.0_8
          ssd(1) = ssd(1) + sf*overcellarea*sfacecoordrefd(1) + 2*ss(1)*&
&           tempd0 + vxm*vnmd
          sfacecoordrefd(1) = 0.0_8
          ptotd = ptotd + pref*cellarea*area_ptotd + massflowratelocal*&
&           pref*mass_ptotd
          rhorefd = rhorefd + rhom*massflowratelocal*mass_rhod
          ttotd = ttotd + massflowratelocal*tref*mass_ttotd
          trefd = trefd + ttot*massflowratelocal*mass_ttotd
          mredimd = mredimd + rhom*vnm*massflowratelocald
          call computettot_b(rhom, rhomd, vxm, vxmd, vym, vymd, vzm, &
&                      vzmd, pm, pmd, ttot, ttotd)
          call computeptot_b(rhom, rhomd, vxm, vxmd, vym, vymd, vzm, &
&                      vzmd, pm, pmd, ptot, ptotd)
          normd = normd + half*ssd
          call cross_prod_b(arg1(:), arg1d(:), arg2(:), arg2d(:), norm, &
&                     normd)
          x1d = 0.0_8
          x3d = 0.0_8
          x3d = arg2d
          x1d = -arg2d - arg1d
          x2d = 0.0_8
          x2d = arg1d
          varsd(conn(3, i), izippflowx:izippflowz) = varsd(conn(3, i), &
&           izippflowx:izippflowz) + x3d
          varsd(conn(2, i), izippflowx:izippflowz) = varsd(conn(2, i), &
&           izippflowx:izippflowz) + x2d
          varsd(conn(1, i), izippflowx:izippflowz) = varsd(conn(1, i), &
&           izippflowx:izippflowz) + x1d
          sfd = third*sfd
          gammamd = third*gammamd
          pmd = third*pmd
          vzmd = third*vzmd
          vymd = third*vymd
          vxmd = third*vxmd
          rhomd = third*rhomd
          do j=3,1,-1
            varsd(conn(j, i), izippflowsface) = varsd(conn(j, i), &
&             izippflowsface) + sfd
            varsd(conn(j, i), izippflowgamma) = varsd(conn(j, i), &
&             izippflowgamma) + gammamd
            varsd(conn(j, i), irhoe) = varsd(conn(j, i), irhoe) + pmd
            varsd(conn(j, i), ivz) = varsd(conn(j, i), ivz) + vzmd
            varsd(conn(j, i), ivy) = varsd(conn(j, i), ivy) + vymd
            varsd(conn(j, i), ivx) = varsd(conn(j, i), ivx) + vxmd
            varsd(conn(j, i), irho) = varsd(conn(j, i), irho) + rhomd
          end do
        end if
        call popcontrol1b(branch)
      end if
    end do
    pointrefd(3) = pointrefd(3) + lref*refpointd(3)
    refpointd(3) = 0.0_8
    pointrefd(2) = pointrefd(2) + lref*refpointd(2)
    refpointd(2) = 0.0_8
    pointrefd(1) = pointrefd(1) + lref*refpointd(1)
    if (pref*rhoref .eq. 0.0_8) then
      tempd = 0.0_8
    else
      tempd = mredimd/(2.0*sqrt(pref*rhoref))
    end if
    prefd = prefd + rhoref*tempd
    rhorefd = rhorefd + pref*tempd
  end subroutine flowintegrationzipper_b
 
  subroutine flowintegrationzipper(isinflow, conn, fams, vars, &
&   localvalues, famlist, sps, ptvalid)
! integrate over the trianges for the inflow/outflow conditions.
    use constants
    use blockpointers, only : bctype
    use sorting, only : faminlist
    use flowvarrefstate, only : pref, pinf, rhoref, pref, timeref, &
&   lref, tref, rgas, uref, uinf, rhoinf, gammainf
    use inputphysics, only : pointref, flowtype, rgasdim
    use flowutils_b, only : computeptot, computettot
    use surfacefamilies, only : familyexchange, bcfamexchange
    use utils_b, only : mynorm2, cross_prod
    implicit none
! input/output variables
    logical, intent(in) :: isinflow
    integer(kind=inttype), dimension(:, :), intent(in) :: conn
    integer(kind=inttype), dimension(:), intent(in) :: fams
    real(kind=realtype), dimension(:, :), intent(in) :: vars
    real(kind=realtype), dimension(nlocalvalues), intent(inout) :: &
&   localvalues
    integer(kind=inttype), dimension(:), intent(in) :: famlist
    integer(kind=inttype), intent(in) :: sps
    logical(kind=inttype), dimension(:), optional, intent(in) :: ptvalid
! working variables
    integer(kind=inttype) :: i, j
    real(kind=realtype) :: sf, vmag, vnm, vxm, vym, vzm, fx, fy, fz, u, &
&   v, w, vnmfreestreamref
    real(kind=realtype), dimension(3) :: fp, mp, fmom, mmom, refpoint, &
&   ss, x1, x2, x3, norm, sfacecoordref
    real(kind=realtype) :: pm, ptot, ttot, rhom, gammam, mnm, &
&   massflowratelocal, am
    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, cellarea, overcellarea
    real(kind=realtype) :: area_ptot, area_ps
    real(kind=realtype) :: govgm1, gm1ovg, viconst, vilocal, pratio
    real(kind=realtype) :: mredim
    real(kind=realtype) :: internalflowfact, inflowfact, xc, xco, yc, &
&   yco, zc, zco, mx, my, mz
    real(kind=realtype), dimension(3) :: cofsumfx, cofsumfy, cofsumfz
    logical :: triisvalid
    intrinsic sqrt
    intrinsic size
    intrinsic present
    intrinsic min
    real(kind=realtype), dimension(3) :: arg1
    real(kind=realtype), dimension(3) :: arg2
    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
    refpoint(1) = lref*pointref(1)
    refpoint(2) = lref*pointref(2)
    refpoint(3) = lref*pointref(3)
    internalflowfact = one
    if (flowtype .eq. internalflow) internalflowfact = -one
    inflowfact = one
    if (isinflow) inflowfact = -one
!$ad ii-loop
    do i=1,size(conn, 2)
      if (faminlist(fams(i), famlist)) then
! if the ptvalid list is given, check if we should integrate
! this triangle.
        triisvalid = .true.
        if (present(ptvalid)) then
! check if each of the three nodes are valid
          if (((ptvalid(conn(1, i)) .eqv. .false.) .or. (ptvalid(conn(2&
&             , i)) .eqv. .false.)) .or. (ptvalid(conn(3, i)) .eqv. &
&             .false.)) triisvalid = .false.
        end if
        if (triisvalid) then
! compute the averaged values for this triangle
          vxm = zero
          vym = zero
          vzm = zero
          rhom = zero
          pm = zero
          mnm = zero
          gammam = zero
          sf = zero
          do j=1,3
            rhom = rhom + vars(conn(j, i), irho)
            vxm = vxm + vars(conn(j, i), ivx)
            vym = vym + vars(conn(j, i), ivy)
            vzm = vzm + vars(conn(j, i), ivz)
            pm = pm + vars(conn(j, i), irhoe)
            gammam = gammam + vars(conn(j, i), izippflowgamma)
            sf = sf + vars(conn(j, i), izippflowsface)
          end do
! divide by 3 due to the summation above:
          rhom = third*rhom
          vxm = third*vxm
          vym = third*vym
          vzm = third*vzm
          pm = third*pm
          gammam = third*gammam
          sf = third*sf
! get the nodes of triangle.
          x1 = vars(conn(1, i), izippflowx:izippflowz)
          x2 = vars(conn(2, i), izippflowx:izippflowz)
          x3 = vars(conn(3, i), izippflowx:izippflowz)
          arg1(:) = x2 - x1
          arg2(:) = x3 - x1
          call cross_prod(arg1(:), arg2(:), norm)
          ss = half*norm
          call computeptot(rhom, vxm, vym, vzm, pm, ptot)
          call computettot(rhom, vxm, vym, vzm, pm, ttot)
          vnm = vxm*ss(1) + vym*ss(2) + vzm*ss(3) - sf
          vmag = sqrt(vxm**2 + vym**2 + vzm**2) - sf
          am = sqrt(gammam*pm/rhom)
          mnm = vmag/sqrt(gammam*pm/rhom)
          cellarea = sqrt(ss(1)**2 + ss(2)**2 + ss(3)**2)
          area = area + cellarea
          overcellarea = 1/cellarea
          massflowratelocal = rhom*vnm*mredim
          massflowrate = massflowrate + massflowratelocal
          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
          area_ps = area_ps + pm*cellarea
          sfacecoordref(1) = sf*ss(1)*overcellarea
          sfacecoordref(2) = sf*ss(2)*overcellarea
          sfacecoordref(3) = sf*ss(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 + ss(1)*overcellarea*massflowratelocal
          mass_ny = mass_ny + ss(2)*overcellarea*massflowratelocal
          mass_nz = mass_nz + ss(3)*overcellarea*massflowratelocal
! compute the average cell center.
          xco = zero
          yco = zero
          zco = zero
          do j=1,3
            xco = xco + vars(conn(1, i), izippflowx)
            yco = yco + vars(conn(2, i), izippflowy)
            zco = zco + vars(conn(3, i), izippflowz)
          end do
! finish average for cell center
          xco = third*xco
          yco = third*yco
          zco = third*zco
! x-y-zco is the cell center w.r.t. the origin, x-y-zc is w.r.t. the reference point
          xc = xco - refpoint(1)
          yc = yco - refpoint(2)
          zc = zco - refpoint(3)
          pm = -(pm-pinf*pref)
          fx = pm*ss(1)
          fy = pm*ss(2)
          fz = pm*ss(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
! center of force computations. here we accumulate in the sums.
! 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
! get unit normal vector.
          ss = ss/cellarea
          massflowratelocal = massflowratelocal/timeref*internalflowfact&
&           *inflowfact
          fx = massflowratelocal*ss(1)*vxm
          fy = massflowratelocal*ss(2)*vym
          fz = massflowratelocal*ss(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.
! 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 if
      end if
    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 flowintegrationzipper
 
!  differentiation of wallintegrationzipper in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: pointref vars localvalues
!   with respect to varying inputs: pointref vars localvalues
!   rw status of diff variables: pointref:incr vars:incr localvalues:in-out
  subroutine wallintegrationzipper_b(conn, fams, vars, varsd, &
&   localvalues, localvaluesd, famlist, sps)
    use constants
    use sorting, only : faminlist
    use flowvarrefstate, only : lref
    use inputphysics, only : pointref, pointrefd
    use utils_b, only : mynorm2, mynorm2_b, cross_prod, cross_prod_b
    implicit none
! input/output
    integer(kind=inttype), dimension(:, :), intent(in) :: conn
    integer(kind=inttype), dimension(:), intent(in) :: fams
    real(kind=realtype), dimension(:, :), intent(in) :: vars
    real(kind=realtype), dimension(:, :) :: varsd
    real(kind=realtype), intent(inout) :: localvalues(nlocalvalues)
    real(kind=realtype), intent(inout) :: localvaluesd(nlocalvalues)
    integer(kind=inttype), dimension(:), intent(in) :: famlist
    integer(kind=inttype), intent(in) :: sps
! working
    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
    integer(kind=inttype) :: i, j
    real(kind=realtype), dimension(3) :: ss, norm, refpoint
    real(kind=realtype), dimension(3) :: ssd, normd, refpointd
    real(kind=realtype), dimension(3) :: p1, p2, p3, v1, v2, v3, x1, x2&
&   , x3
    real(kind=realtype), dimension(3) :: p1d, p2d, p3d, v1d, v2d, v3d, &
&   x1d, x2d, x3d
    real(kind=realtype) :: fact, triarea, fx, fy, fz, mx, my, mz, xc, &
&   xco, yc, yco, zc, zco
    real(kind=realtype) :: triaread, fxd, fyd, fzd, mxd, myd, mzd, xcd, &
&   xcod, ycd, ycod, zcd, zcod
    intrinsic size
    real(kind=realtype), dimension(3) :: arg1
    real(kind=realtype), dimension(3) :: arg1d
    real(kind=realtype), dimension(3) :: arg2
    real(kind=realtype), dimension(3) :: arg2d
    real(kind=realtype) :: result1
    real(kind=realtype) :: result1d
    real(kind=realtype) :: tempd
! 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)
    cofsumfzd = 0.0_8
    cofsumfzd = localvaluesd(icoforcez:icoforcez+2)
    cofsumfyd = 0.0_8
    cofsumfyd = localvaluesd(icoforcey:icoforcey+2)
    cofsumfxd = 0.0_8
    cofsumfxd = localvaluesd(icoforcex:icoforcex+2)
    mvd = 0.0_8
    mvd = localvaluesd(imv:imv+2)
    mpd = 0.0_8
    mpd = localvaluesd(imp:imp+2)
    fvd = 0.0_8
    fvd = localvaluesd(ifv:ifv+2)
    fpd = 0.0_8
    fpd = localvaluesd(ifp:ifp+2)
    normd = 0.0_8
    refpointd = 0.0_8
!$bwd-of ii-loop 
    do i=1,size(conn, 2)
      if (faminlist(fams(i), famlist)) then
! get the nodes of triangle.
        x1 = vars(conn(1, i), izippwallx:izippwallz)
        x2 = vars(conn(2, i), izippwallx:izippwallz)
        x3 = vars(conn(3, i), izippwallx:izippwallz)
        arg1(:) = x2 - x1
        arg2(:) = x3 - x1
        call cross_prod(arg1(:), arg2(:), norm)
        ss = half*norm
! the third here is to account for the summation of p1, p2
! and p3
        result1 = mynorm2(ss)
        triarea = result1*third
! compute the average cell center.
        xco = third*(x1(1)+x2(1)+x3(1))
        yco = third*(x1(2)+x2(2)+x3(2))
        zco = third*(x1(3)+x2(3)+x3(3))
        xc = xco - refpoint(1)
        yc = yco - refpoint(2)
        zc = zco - refpoint(3)
! update the pressure force and moment coefficients.
        p1 = vars(conn(1, i), izippwalltpx:izippwalltpz)
        p2 = vars(conn(2, i), izippwalltpx:izippwalltpz)
        p3 = vars(conn(3, i), izippwalltpx:izippwalltpz)
! accumulate in the sums. each force component is tracked separately
! force-x
! force-y
! force-z
! update the viscous force and moment coefficients
        v1 = vars(conn(1, i), izippwalltvx:izippwalltvz)
        v2 = vars(conn(2, i), izippwalltvx:izippwalltvz)
        v3 = vars(conn(3, i), izippwalltvx:izippwalltvz)
        fx = (v1(1)+v2(1)+v3(1))*triarea
        fy = (v1(2)+v2(2)+v3(2))*triarea
        fz = (v1(3)+v2(3)+v3(3))*triarea
! note: momentum forces have opposite sign to pressure forces
! accumulate in the sums. each force component is tracked separately
! force-x
! force-y
! force-z
        mxd = mvd(1)
        myd = mvd(2)
        mzd = mvd(3)
        zcod = fz*cofsumfzd(3) + fy*cofsumfyd(3) + fx*cofsumfxd(3)
        fzd = zco*cofsumfzd(3) + yco*cofsumfzd(2) + xco*cofsumfzd(1) + &
&         yc*mxd - xc*myd + fvd(3)
        ycod = fz*cofsumfzd(2) + fy*cofsumfyd(2) + fx*cofsumfxd(2)
        xcod = fz*cofsumfzd(1) + fy*cofsumfyd(1) + fx*cofsumfxd(1)
        fyd = zco*cofsumfyd(3) + yco*cofsumfyd(2) + xco*cofsumfyd(1) + &
&         xc*mzd + fvd(2) - zc*mxd
        fxd = zco*cofsumfxd(3) + yco*cofsumfxd(2) + xco*cofsumfxd(1) + &
&         zc*myd - yc*mzd + fvd(1)
        xcd = fy*mzd - fz*myd
        ycd = fz*mxd - fx*mzd
        zcd = fx*myd - fy*mxd
        v1d = 0.0_8
        v2d = 0.0_8
        v3d = 0.0_8
        tempd = triarea*fzd
        triaread = (v1(3)+v2(3)+v3(3))*fzd + (v1(2)+v2(2)+v3(2))*fyd + (&
&         v1(1)+v2(1)+v3(1))*fxd
        v1d(3) = v1d(3) + tempd
        v2d(3) = v2d(3) + tempd
        v3d(3) = v3d(3) + tempd
        tempd = triarea*fyd
        v1d(2) = v1d(2) + tempd
        v2d(2) = v2d(2) + tempd
        v3d(2) = v3d(2) + tempd
        tempd = triarea*fxd
        v1d(1) = v1d(1) + tempd
        v2d(1) = v2d(1) + tempd
        v3d(1) = v3d(1) + tempd
        varsd(conn(3, i), izippwalltvx:izippwalltvz) = varsd(conn(3, i)&
&         , izippwalltvx:izippwalltvz) + v3d
        varsd(conn(2, i), izippwalltvx:izippwalltvz) = varsd(conn(2, i)&
&         , izippwalltvx:izippwalltvz) + v2d
        varsd(conn(1, i), izippwalltvx:izippwalltvz) = varsd(conn(1, i)&
&         , izippwalltvx:izippwalltvz) + v1d
        fz = (p1(3)+p2(3)+p3(3))*triarea
        fy = (p1(2)+p2(2)+p3(2))*triarea
        fx = (p1(1)+p2(1)+p3(1))*triarea
        mzd = mpd(3)
        myd = mpd(2)
        fxd = zco*cofsumfxd(3) + yco*cofsumfxd(2) + xco*cofsumfxd(1) + &
&         zc*myd - yc*mzd + fpd(1)
        mxd = mpd(1)
        fzd = zco*cofsumfzd(3) + yco*cofsumfzd(2) + xco*cofsumfzd(1) + &
&         yc*mxd - xc*myd + fpd(3)
        fyd = zco*cofsumfyd(3) + yco*cofsumfyd(2) + xco*cofsumfyd(1) + &
&         xc*mzd + fpd(2) - zc*mxd
        xcd = xcd + fy*mzd - fz*myd
        xcod = xcod + fz*cofsumfzd(1) + fy*cofsumfyd(1) + fx*cofsumfxd(1&
&         ) + xcd
        ycd = ycd + fz*mxd - fx*mzd
        ycod = ycod + fz*cofsumfzd(2) + fy*cofsumfyd(2) + fx*cofsumfxd(2&
&         ) + ycd
        zcd = zcd + fx*myd - fy*mxd
        zcod = zcod + fz*cofsumfzd(3) + fy*cofsumfyd(3) + fx*cofsumfxd(3&
&         ) + zcd
        p1d = 0.0_8
        p2d = 0.0_8
        p3d = 0.0_8
        tempd = triarea*fzd
        triaread = triaread + (p1(3)+p2(3)+p3(3))*fzd + (p1(2)+p2(2)+p3(&
&         2))*fyd + (p1(1)+p2(1)+p3(1))*fxd
        p1d(3) = p1d(3) + tempd
        p2d(3) = p2d(3) + tempd
        p3d(3) = p3d(3) + tempd
        tempd = triarea*fyd
        p1d(2) = p1d(2) + tempd
        p2d(2) = p2d(2) + tempd
        p3d(2) = p3d(2) + tempd
        tempd = triarea*fxd
        p1d(1) = p1d(1) + tempd
        p2d(1) = p2d(1) + tempd
        p3d(1) = p3d(1) + tempd
        varsd(conn(3, i), izippwalltpx:izippwalltpz) = varsd(conn(3, i)&
&         , izippwalltpx:izippwalltpz) + p3d
        varsd(conn(2, i), izippwalltpx:izippwalltpz) = varsd(conn(2, i)&
&         , izippwalltpx:izippwalltpz) + p2d
        varsd(conn(1, i), izippwalltpx:izippwalltpz) = varsd(conn(1, i)&
&         , izippwalltpx:izippwalltpz) + p1d
        refpointd(3) = refpointd(3) - zcd
        refpointd(2) = refpointd(2) - ycd
        refpointd(1) = refpointd(1) - xcd
        x1d = 0.0_8
        x2d = 0.0_8
        x3d = 0.0_8
        tempd = third*zcod
        x1d(3) = x1d(3) + tempd
        x2d(3) = x2d(3) + tempd
        x3d(3) = x3d(3) + tempd
        tempd = third*ycod
        x1d(2) = x1d(2) + tempd
        x2d(2) = x2d(2) + tempd
        x3d(2) = x3d(2) + tempd
        tempd = third*xcod
        x1d(1) = x1d(1) + tempd
        x2d(1) = x2d(1) + tempd
        x3d(1) = x3d(1) + tempd
        result1d = third*triaread
        call mynorm2_b(ss, ssd, result1d)
        normd = normd + half*ssd
        call cross_prod_b(arg1(:), arg1d(:), arg2(:), arg2d(:), norm, &
&                   normd)
        x3d = x3d + arg2d
        x1d = x1d - arg2d - arg1d
        x2d = x2d + arg1d
        varsd(conn(3, i), izippwallx:izippwallz) = varsd(conn(3, i), &
&         izippwallx:izippwallz) + x3d
        varsd(conn(2, i), izippwallx:izippwallz) = varsd(conn(2, i), &
&         izippwallx:izippwallz) + x2d
        varsd(conn(1, i), izippwallx:izippwallz) = varsd(conn(1, i), &
&         izippwallx:izippwallz) + x1d
      end if
    end do
    pointrefd(3) = pointrefd(3) + lref*refpointd(3)
    refpointd(3) = 0.0_8
    pointrefd(2) = pointrefd(2) + lref*refpointd(2)
    refpointd(2) = 0.0_8
    pointrefd(1) = pointrefd(1) + lref*refpointd(1)
  end subroutine wallintegrationzipper_b
 
  subroutine wallintegrationzipper(conn, fams, vars, localvalues, &
&   famlist, sps)
    use constants
    use sorting, only : faminlist
    use flowvarrefstate, only : lref
    use inputphysics, only : pointref
    use utils_b, only : mynorm2, cross_prod
    implicit none
! input/output
    integer(kind=inttype), dimension(:, :), intent(in) :: conn
    integer(kind=inttype), dimension(:), intent(in) :: fams
    real(kind=realtype), dimension(:, :), intent(in) :: vars
    real(kind=realtype), intent(inout) :: localvalues(nlocalvalues)
    integer(kind=inttype), dimension(:), intent(in) :: famlist
    integer(kind=inttype), intent(in) :: sps
! working
    real(kind=realtype), dimension(3) :: fp, fv, mp, mv
    real(kind=realtype), dimension(3) :: cofsumfx, cofsumfy, cofsumfz
    integer(kind=inttype) :: i, j
    real(kind=realtype), dimension(3) :: ss, norm, refpoint
    real(kind=realtype), dimension(3) :: p1, p2, p3, v1, v2, v3, x1, x2&
&   , x3
    real(kind=realtype) :: fact, triarea, fx, fy, fz, mx, my, mz, xc, &
&   xco, yc, yco, zc, zco
    intrinsic size
    real(kind=realtype), dimension(3) :: arg1
    real(kind=realtype), dimension(3) :: arg2
    real(kind=realtype) :: result1
! 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)
    fp = zero
    fv = zero
    mp = zero
    mv = zero
    cofsumfx = zero
    cofsumfy = zero
    cofsumfz = zero
!$ad ii-loop
    do i=1,size(conn, 2)
      if (faminlist(fams(i), famlist)) then
! get the nodes of triangle.
        x1 = vars(conn(1, i), izippwallx:izippwallz)
        x2 = vars(conn(2, i), izippwallx:izippwallz)
        x3 = vars(conn(3, i), izippwallx:izippwallz)
        arg1(:) = x2 - x1
        arg2(:) = x3 - x1
        call cross_prod(arg1(:), arg2(:), norm)
        ss = half*norm
! the third here is to account for the summation of p1, p2
! and p3
        result1 = mynorm2(ss)
        triarea = result1*third
! compute the average cell center.
        xco = third*(x1(1)+x2(1)+x3(1))
        yco = third*(x1(2)+x2(2)+x3(2))
        zco = third*(x1(3)+x2(3)+x3(3))
        xc = xco - refpoint(1)
        yc = yco - refpoint(2)
        zc = zco - refpoint(3)
! update the pressure force and moment coefficients.
        p1 = vars(conn(1, i), izippwalltpx:izippwalltpz)
        p2 = vars(conn(2, i), izippwalltpx:izippwalltpz)
        p3 = vars(conn(3, i), izippwalltpx:izippwalltpz)
        fx = (p1(1)+p2(1)+p3(1))*triarea
        fy = (p1(2)+p2(2)+p3(2))*triarea
        fz = (p1(3)+p2(3)+p3(3))*triarea
        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
! accumulate in the sums. each force component is tracked separately
! 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
! update the viscous force and moment coefficients
        v1 = vars(conn(1, i), izippwalltvx:izippwalltvz)
        v2 = vars(conn(2, i), izippwalltvx:izippwalltvz)
        v3 = vars(conn(3, i), izippwalltvx:izippwalltvz)
        fx = (v1(1)+v2(1)+v3(1))*triarea
        fy = (v1(2)+v2(2)+v3(2))*triarea
        fz = (v1(3)+v2(3)+v3(3))*triarea
! note: momentum forces have opposite sign to pressure forces
        fv(1) = fv(1) + fx
        fv(2) = fv(2) + fy
        fv(3) = fv(3) + fz
        mx = yc*fz - zc*fy
        my = zc*fx - xc*fz
        mz = xc*fy - yc*fx
        mv(1) = mv(1) + mx
        mv(2) = mv(2) + my
        mv(3) = mv(3) + mz
! accumulate in the sums. each force component is tracked separately
! 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 if
    end do
! increment into the local vector
    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
  end subroutine wallintegrationzipper
 
end module zipperintegrations_b