zipperIntegrations.F90

Go to the documentation of this file.

module zipperintegrations
contains
 
    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, only: computeptot, computettot
        use surfacefamilies, only: familyexchange, bcfamexchange
        use utils, only: mynorm2, cross_prod
        implicit none
 
        ! Input/output Variables
        logical, intent(in) :: isInflow
        integer(kind=intType), intent(in), dimension(:, :) :: conn
        integer(kind=intType), intent(in), dimension(:) :: 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
 
        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 == internalflow) then
            internalflowfact = -one
        end if
 
        inflowfact = one
        if (isinflow) then
            inflowfact = -one
        end if
 
        !$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.)) then
                        triisvalid = .false.
                    end if
                end if
 
                validtrianlge: 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)
                    call cross_prod(x2 - x1, x3 - x1, 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
                    ! the prefs in psinf / ptot cancel out so we can just take the ratio
                    ! we need to clip the ratio to stay under one. right next to the wall,
                    ! the pTot can go below the static free stream pressure. To prevent
                    ! nans from the sqrt, we just clip this. This does not affect the computation
                    ! because when pTot is this small, the velocities are also small, and the
                    ! mdot is almost zero, so the cells in this area don't contribute much
                    ! to the mass weighed sum.
                    pratio = min(one, one / ptot)
                    vilocal = sqrt(viconst * (one - (pratio)**gm1ovg) * ttot * tref)
                    mass_vi = mass_vi + vilocal * massflowratelocal
 
                    mass_nx = mass_nx + 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 validtrianlge
            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
 
    subroutine wallintegrationzipper(conn, fams, vars, localValues, famList, sps)
 
        use constants
        use sorting, only: faminlist
        use flowvarrefstate, only: lref
        use inputphysics, only: pointref
        use utils, only: mynorm2, cross_prod
        implicit none
 
        ! Input/Output
        integer(kind=intType), intent(in), dimension(:, :) :: conn
        integer(kind=intType), intent(in), dimension(:) :: fams
        real(kind=realtype), intent(in), dimension(:, :) :: 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
 
        ! 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)
                call cross_prod(x2 - x1, x3 - x1, norm)
                ss = half * norm
                ! The third here is to account for the summation of P1, p2
                ! and P3
                triarea = mynorm2(ss) * 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
 
#ifndef USE_TAPENADE
 
    subroutine integratezippers(localValues, famList, sps)
        !--------------------------------------------------------------
        ! Manual Differentiation Warning: Modifying this routine requires
        ! modifying the hand-written forward and reverse routines.
        ! --------------------------------------------------------------
 
        ! Integrate over the triangles formed by the zipper mesh. This
        ! will perform both all necesasry zipper integrations. Currently
        ! this includes the wall force integrations as well as the
        ! flow-though surface integration.
 
        use constants
        use oversetdata, only: zippermeshes, zippermesh
        use haloexchange, only: wallintegrationzippercomm, flowintegrationzippercomm
        implicit none
 
        ! Input Variables
        real(kind=realtype), dimension(nLocalValues), intent(inout) :: localvalues
        integer(kind=intType), dimension(:), intent(in) :: famList
        integer(kind=intType), intent(in) :: sps
        real(kind=realtype), dimension(:, :), allocatable :: vars
        type(zippermesh), pointer :: zipper
 
        ! Determine if we have a wall Zipper:
        zipper => zippermeshes(ibcgroupwalls)
        if (zipper%allocated) then
 
            ! Allocate space necessary to store variables. Only non-zero on
            ! root proc.
            allocate (vars(size(zipper%indices), nzippwallcomm))
 
            ! Gather up the required variables in "vars" on the root
            ! proc. This routine is differientated by hand.
            call wallintegrationzippercomm(vars, sps)
 
            ! Perform actual integration. Tapenade ADs this routine.
            call wallintegrationzipper(zipper%conn, zipper%fam, vars, localvalues, famlist, sps)
 
            ! Cleanup vars
            deallocate (vars)
        end if
 
        zipper => zippermeshes(ibcgroupinflow)
        ! Determine if we have a flowthrough Zipper:
        if (zipper%allocated) then
 
            ! Allocate space necessary to store variables. Only non-zero on
            ! root proc.
            allocate (vars(size(zipper%indices), nzippflowcomm))
 
            ! Gather up the required variables in "vars" on the root
            ! proc. This routine is differientated by hand.
            call flowintegrationzippercomm(.true., vars, sps)
 
            ! Perform actual integration. Tapenade ADs this routine.
            call flowintegrationzipper(.true., zipper%conn, zipper%fam, vars, &
                                       localvalues, famlist, sps)
 
            ! Cleanup vars
            deallocate (vars)
        end if
 
        zipper => zippermeshes(ibcgroupoutflow)
        ! Determine if we have a flowthrough Zipper:
        if (zipper%allocated) then
 
            ! Allocate space necessary to store variables. Only non-zero on
            ! root proc.
            allocate (vars(size(zipper%indices), nzippflowcomm))
 
            ! Gather up the required variables in "vars" on the root
            ! proc. This routine is differientated by hand.
            call flowintegrationzippercomm(.false., vars, sps)
 
            ! Perform actual integration. Tapenade ADs this routine.
            call flowintegrationzipper(.false., zipper%conn, zipper%fam, vars, &
                                       localvalues, famlist, sps)
 
            ! Cleanup vars
            deallocate (vars)
        end if
    end subroutine integratezippers
 
#ifndef USE_COMPLEX
    subroutine integratezippers_d(localValues, localValuesd, famList, sps)
        !------------------------------------------------------------------------
        ! Manual Differentiation Warning: This routine is differentiated by hand.
        ! -----------------------------------------------------------------------
 
        ! Forward mode linearization of the zipper integration.
 
        use constants
        use oversetdata, only: zippermeshes, zippermesh
        use haloexchange, only: wallintegrationzippercomm_d, flowintegrationzippercomm_d
        use zipperintegrations_d, only: flowintegrationzipper_d, wallintegrationzipper_d
        implicit none
 
        ! Input Variables
        real(kind=realtype), dimension(nLocalValues), intent(inout) :: localvalues, localvaluesd
        integer(kind=intType), dimension(:), intent(in) :: famlist
        integer(kind=intType), intent(in) :: sps
        real(kind=realtype), dimension(:, :), allocatable :: vars, varsd
        type(zippermesh), pointer :: zipper
 
        zipper => zippermeshes(ibcgroupwalls)
        if (zipper%allocated) then
 
            ! Allocate space necessary to store variables. Only non-zero on
            ! root proc.
            allocate (vars(size(zipper%indices), nzippwallcomm), varsd(size(zipper%indices), nzippwallcomm))
 
            ! Gather up the required variables in "vars" on the root
            ! proc. This routine is differientated by hand.
            call wallintegrationzippercomm_d(vars, varsd, sps)
 
            ! Perform actual integration. Tapenade ADs this routine.
            call wallintegrationzipper_d(zipper%conn, zipper%fam, vars, varsd, &
                                         localvalues, localvaluesd, famlist, sps)
 
            ! Cleanup vars
            deallocate (vars, varsd)
        end if
 
        zipper => zippermeshes(ibcgroupinflow)
        ! Determine if we have a flowthrough Zipper:
        if (zipper%allocated) then
 
            ! Allocate space necessary to store variables. Only non-zero on
            ! root proc.
            allocate (vars(size(zipper%indices), nzippflowcomm), varsd(size(zipper%indices), nzippflowcomm))
 
            ! Gather up the required variables in "vars" on the root
            ! proc. This routine is differientated by hand.
            call flowintegrationzippercomm_d(.true., vars, varsd, sps)
 
            ! Perform actual integration. Tapenade ADs this routine.
            call flowintegrationzipper_d(.true., zipper%conn, zipper%fam, vars, varsd, &
                                         localvalues, localvaluesd, famlist, sps)
 
            ! Cleanup vars
            deallocate (vars, varsd)
        end if
 
        zipper => zippermeshes(ibcgroupoutflow)
        ! Determine if we have a flowthrough Zipper:
        if (zipper%allocated) then
 
            ! Allocate space necessary to store variables. Only non-zero on
            ! root proc.
            allocate (vars(size(zipper%indices), nzippflowcomm), varsd(size(zipper%indices), nzippflowcomm))
 
            ! Gather up the required variables in "vars" on the root
            ! proc. This routine is differientated by hand.
            call flowintegrationzippercomm_d(.false., vars, varsd, sps)
 
            ! Perform actual integration. Tapenade ADs this routine.
            call flowintegrationzipper_d(.false., zipper%conn, zipper%fam, vars, varsd, &
                                         localvalues, localvaluesd, famlist, sps)
 
            ! Cleanup vars
            deallocate (vars, varsd)
        end if
    end subroutine integratezippers_d
 
    subroutine integratezippers_b(localValues, localValuesd, famList, sps)
        !------------------------------------------------------------------------
        ! Manual Differentiation Warning: This routine is differentiated by hand.
        ! -----------------------------------------------------------------------
 
        ! Reverse mode linearization of the zipper integration.
 
        use constants
        use oversetdata, only: zippermeshes, zippermesh
        use haloexchange, only: wallintegrationzippercomm_b, flowintegrationzippercomm_b, &
                                wallintegrationzippercomm, flowintegrationzippercomm
        use zipperintegrations_b, only: wallintegrationzipper_b, flowintegrationzipper_b
        implicit none
 
        ! Input Variables
        real(kind=realtype), dimension(nLocalValues), intent(inout) :: localvalues, localvaluesd
        integer(kind=intType), dimension(:), intent(in) :: famlist
        integer(kind=intType), intent(in) :: sps
        real(kind=realtype), dimension(:, :), allocatable :: vars, varsd
        type(zippermesh), pointer :: zipper
 
        ! Determine if we have a wall Zipper:
        zipper => zippermeshes(ibcgroupwalls)
        if (zipper%allocated) then
 
            ! Allocate space necessary to store variables. Only non-zero on
            ! root proc.
            allocate (vars(size(zipper%indices), nzippwallcomm), varsd(size(zipper%indices), nzippwallcomm))
 
            ! Set the forward variables
            call wallintegrationzippercomm(vars, sps)
 
            ! Perform actual integration. Tapenade ADs this routine.
            varsd = zero
            call wallintegrationzipper_b(zipper%conn, zipper%fam, vars, varsd, &
                                         localvalues, localvaluesd, famlist, sps)
 
            ! Scatter (becuase we are reverse) the values from the root
            ! back out to all necessary procs.  This routine is
            ! differientated by hand.
            call wallintegrationzippercomm_b(vars, varsd, sps)
 
            ! Cleanup vars
            deallocate (vars, varsd)
        end if
 
        zipper => zippermeshes(ibcgroupinflow)
        ! Determine if we have a flowthrough Zipper:
        if (zipper%allocated) then
 
            ! Allocate space necessary to store variables. Only non-zero on
            ! root proc.
            allocate (vars(size(zipper%indices), nzippflowcomm), varsd(size(zipper%indices), nzippflowcomm))
 
            ! Set the forward variables
            call flowintegrationzippercomm(.true., vars, sps)
 
            ! Perform actual integration. Tapenade ADs this routine.
            varsd = zero
            call flowintegrationzipper_b(.true., zipper%conn, zipper%fam, vars, varsd, &
                                         localvalues, localvaluesd, famlist, sps)
 
            ! Scatter (becuase we are reverse) the values from the root
            ! back out to all necessary procs.  This routine is
            ! differientated by hand.
            call flowintegrationzippercomm_b(.true., vars, varsd, sps)
 
            ! Cleanup vars
            deallocate (vars, varsd)
        end if
 
        zipper => zippermeshes(ibcgroupoutflow)
        ! Determine if we have a flowthrough Zipper:
        if (zipper%allocated) then
 
            ! Allocate space necessary to store variables. Only non-zero on
            ! root proc.
            allocate (vars(size(zipper%indices), nzippflowcomm), varsd(size(zipper%indices), nzippflowcomm))
 
            ! Set the forward variables
            call flowintegrationzippercomm(.false., vars, sps)
 
            ! Perform actual integration. Tapenade ADs this routine.
            varsd = zero
            call flowintegrationzipper_b(.false., zipper%conn, zipper%fam, vars, varsd, &
                                         localvalues, localvaluesd, famlist, sps)
 
            ! Scatter (becuase we are reverse) the values from the root
            ! back out to all necessary procs.  This routine is
            ! differientated by hand.
            call flowintegrationzippercomm_b(.false., vars, varsd, sps)
 
            ! Cleanup vars
            deallocate (vars, varsd)
        end if
    end subroutine integratezippers_b
#endif
#endif
 
end module zipperintegrations