wallDistance.F90

Go to the documentation of this file.

module walldistance
 
    use constants, only: inttype, realtype
    use walldistancedata
    implicit none
    save
 
#ifndef USE_TAPENADE
    ! nquadVisc:             Number of local quads on the viscous
    !                        bodies.
    ! nNodeVisc:             Number of local nodes on the viscous
    !                        bodies.
    ! nquadViscGlob:         Global number of viscous quads.
    ! connVisc(4,nquadVisc): Connectivity of the local viscous
    !                        quads.
    ! coorVisc(3,nNodeVisc): The coordinates of the local nodes.
 
    integer(kind=intType) :: nquadvisc, nnodevisc
    integer(kind=intType) :: nquadviscglob
 
    integer(kind=intType), dimension(:, :), allocatable :: connvisc
    real(kind=realtype), dimension(:, :), allocatable :: coorvisc
 
    ! rotMatrixSections(nSections,3,3): Rotation matrices needed
    !                                   for the alignment of the
    !                                   sections. The rotation
    !                                   matrix is a**n, where a is
    !                                   periodic transformation
    !                                   matrix; n is an integer.
 
    real(kind=realtype), dimension(:, :, :), allocatable :: rotmatrixsections
#endif
 
contains
 
    subroutine updatewalldistancesquickly(nn, level, sps)
 
        ! This is the actual update routine that uses xSurf. It is done on
        ! block-level-sps basis.  This is the used to update the wall
        ! distance. Most importantly, this routine is included in the
        ! reverse mode AD routines, but NOT the forward mode. Since it is
        ! done on a per-block basis, it is assumed that the required block
        ! pointers are already set.
 
        use constants
        use blockpointers, only: nx, ny, nz, il, jl, kl, x, flowdoms, d2wall
        implicit none
 
        ! Subroutine arguments
        integer(kind=intType) :: nn, level, sps
 
        ! Local Variables
        integer(kind=intType) :: i, j, k, ii, ind(4)
        real(kind=realtype) :: xp(3), xc(3), u, v
 
#ifdef TAPENADE_REVERSE
        !$AD II-LOOP
        do ii = 0, nx * ny * nz - 1
            i = mod(ii, nx) + 2
            j = mod(ii / nx, ny) + 2
            k = ii / (nx * ny) + 2
#else
            do k = 2, kl
                do j = 2, jl
                    do i = 2, il
#endif
 
                        if (flowdoms(nn, level, sps)%surfNodeIndices(1, i, j, k) == 0) then
                            ! This node is too far away and has no
                            ! association. Set the distance to a large constant.
                            d2wall(i, j, k) = large
                            cycle
                        end if
 
                        ! Extract elemID and u-v position for the association of
                        ! this cell:
 
                        ind = flowdoms(nn, level, sps)%surfNodeIndices(:, i, j, k)
                        u = flowdoms(nn, level, sps)%uv(1, i, j, k)
                        v = flowdoms(nn, level, sps)%uv(2, i, j, k)
 
                        ! Now we have the 4 corners, use bi-linear shape
                        ! functions o to get target: (CCW ordering remember!)
 
                        xp(:) = &
                            (one - u) * (one - v) * xsurf(3 * (ind(1) - 1) + 1:3 * ind(1)) + &
                            (u) * (one - v) * xsurf(3 * (ind(2) - 1) + 1:3 * ind(2)) + &
                            (u) * (v) * xsurf(3 * (ind(3) - 1) + 1:3 * ind(3)) + &
                            (one - u) * (v) * xsurf(3 * (ind(4) - 1) + 1:3 * ind(4))
 
                        ! Get the cell center
                        xc(1) = eighth * (x(i - 1, j - 1, k - 1, 1) + x(i, j - 1, k - 1, 1) &
                                          + x(i - 1, j, k - 1, 1) + x(i, j, k - 1, 1) &
                                          + x(i - 1, j - 1, k, 1) + x(i, j - 1, k, 1) &
                                          + x(i - 1, j, k, 1) + x(i, j, k, 1))
 
                        xc(2) = eighth * (x(i - 1, j - 1, k - 1, 2) + x(i, j - 1, k - 1, 2) &
                                          + x(i - 1, j, k - 1, 2) + x(i, j, k - 1, 2) &
                                          + x(i - 1, j - 1, k, 2) + x(i, j - 1, k, 2) &
                                          + x(i - 1, j, k, 2) + x(i, j, k, 2))
 
                        xc(3) = eighth * (x(i - 1, j - 1, k - 1, 3) + x(i, j - 1, k - 1, 3) &
                                          + x(i - 1, j, k - 1, 3) + x(i, j, k - 1, 3) &
                                          + x(i - 1, j - 1, k, 3) + x(i, j - 1, k, 3) &
                                          + x(i - 1, j, k, 3) + x(i, j, k, 3))
 
                        ! Now we have the two points...just take the norm of the
                        ! distance between them
 
                        d2wall(i, j, k) = sqrt( &
                                          (xc(1) - xp(1))**2 + (xc(2) - xp(2))**2 + (xc(3) - xp(3))**2)
#ifdef TAPENADE_REVERSE
                    end do
#else
                end do
            end do
        end do
#endif
 
    end subroutine updatewalldistancesquickly
 
    ! ----------------------------------------------------------------------
    !                                                                      |
    !                    No Tapenade Routine below this line               |
    !                                                                      |
    ! ----------------------------------------------------------------------
 
#ifndef USE_TAPENADE
    subroutine computewalldistance(level, allocMem)
        !
        !       wallDistance computes the distances of the cell centers to
        !       the nearest viscous wall. An adt type of method is used, which
        !       guarantees to find the minimum distance to the wall. Possible
        !       periodic transformations are taken into account, such that
        !       also in case of a periodic problem the correct distance is
        !       computed; the nearest wall point may lie in a periodic domain.
        use constants
        use blockpointers, only: ndom
        use communication, only: sendbuffer, recvbuffer, myid, adflow_comm_world, &
                                 sendbuffersize, recvbuffersize
        use inputphysics, only: equations, walldistanceneeded
        use inputtimespectral, only: ntimeintervalsspectral
        use inputdiscretization, only: useapproxwalldistance, updatewallassociations
        use oversetdata, only: oversetpresent
        use utils, only: setpointers, echk, terminate, &
                         deallocatetempmemory, allocatetempmemory
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: level
        logical, intent(in) :: allocMem
        !
        !      Local variables.
        !
        integer :: ierr, i, j, k, nn, ii, l
 
        integer(kind=intType) :: sps, sps2, ll, nLevels
        logical :: tempLogical
        real(kind=alwaysrealtype) :: t0
        character(len=3) :: integerString
 
        ! Check if the RANS equations are solved. If not, the wall
        ! distance is not needed and a return can be made.
 
        if (equations /= ransequations) return
 
        ! If the turbulence model is wall distance free just compute the
        ! normal spacing of the first cell and store this in the wall
        ! distance. It may be needed for the boundary conditions and
        ! the monitoring of the y+. Return afterwards.
 
        if (.not. walldistanceneeded) then
 
            ! Loop over the number of spectral solutions, initialize the
            ! distance and compute the initial normal spacing.
 
            do sps = 1, ntimeintervalsspectral
                call initwalldistance(level, sps, allocmem)
                call computenormalspacing(level, sps)
            end do
 
            ! And return.
 
            return
        end if
 
        ! Write a message to stdout to indicate that the wall distance
        ! computation starts for the given level.
 
        write (integerstring, "(i2)") level
        integerstring = adjustl(integerstring)
 
        ! Store the start time.
 
        t0 = mpi_wtime()
 
        ! Release temporarily some memory such that the overall memory
        ! requirement is not dictated by this routine. What memory is
        ! released depends on allocMem. If this is .True. It means that
        ! this is the first time the wall distances are computed, i.e. in
        ! the preprocessing phase. Then only send and receive buffers are
        ! released. If allocMem is .False., this means that the wall
        ! distances are computed in a moving mesh computation and
        ! consequently some more memory can be released.
 
        if (allocmem) then
            deallocate (sendbuffer, recvbuffer, stat=ierr)
            if (ierr /= 0) &
                call terminate("wallDistance", &
                               "Deallocation error for communication buffers")
        else
            call deallocatetempmemory(.true.)
        end if
 
        ! There are two different searches we can do: the original code
        ! always works and it capable to dealing with rotating/periodic
        ! geometries. It uses constant memory and is slow. The
        ! alternative method uses memory that scales with the size of
        ! the surface grid per processor and only works for
        ! steady/unsteady simulations without periodic/rotating
        ! components. But it is fast. It is designed to be used for
        ! updating the wall distances between iterations of
        ! aerostructural solutions.
 
        ! Normal, original wall distance calc. Cannot be used when
        ! overset is present due to possibility of overlapping walls.
        if ((.not. useapproxwalldistance) .and. (.not. oversetpresent)) then
            ! Loop over the number of spectral solutions.
            spectralloop: do sps = 1, ntimeintervalsspectral
 
                ! Initialize the wall distances.
 
                call initwalldistance(level, sps, allocmem)
 
                ! Build the viscous surface mesh.
 
                call viscoussurfacemesh(level, sps)
 
                ! If there are no viscous faces, processor 0 prints a warning
                ! and the wall distances are not computed.
 
                if (nquadviscglob == 0) then
 
                    if (myid == 0) then
                        print "(a)", "#"
                        print "(a)", "#               Warning!!!!"
                        print "(a)", "# No viscous boundary found. Wall &
                             &distances are set to infinity"
                        print "(a)", "#"
                    end if
 
                else
                    ! Determine the wall distances for the owned cells.
                    call determinedistance(level, sps)
                end if
            end do spectralloop
        else ! The user wants to use approx wall distance calcs OR we
            ! have overset mesh. :
 
            if (updatelevelwallassociation(level)) then
 
                ! Initialize the wall distance
                spectralloop2: do sps = 1, ntimeintervalsspectral
                    call initwalldistance(level, sps, allocmem)
                end do spectralloop2
 
                ! Destroy the PETSc wall distance data if necessary
                call destroywalldistancedatalevel(level)
 
                ! Do the associtaion. This allocates the data destroyed
                ! in the destroyWallDistanceData call
 
                do sps = 1, ntimeintervalsspectral
                    call determinewallassociation(level, sps)
                end do
 
                if (.not. updatewallassociations) then
                    ! don't re-compute the wall associations after the first call
                    updatelevelwallassociation(level) = .false.
                end if
            end if
 
            ! Update the xsurf vector from X
            call updatexsurf(level)
 
            ! Call the actual update routine, on each of the sps instances and blocks
            do sps = 1, ntimeintervalsspectral
 
                ! Now extract the vector of the surface data we need
                call vecgetarrayf90(xsurfvec(level, sps), xsurf, ierr)
                call echk(ierr, __file__, __line__)
 
                do nn = 1, ndom
                    call setpointers(nn, level, sps)
                    call updatewalldistancesquickly(nn, level, sps)
                end do
 
                call vecrestorearrayf90(xsurfvec(level, sps), xsurf, ierr)
                call echk(ierr, __file__, __line__)
 
            end do
        end if
 
        ! Allocate the temporarily released memory again. For more info
        ! see the comments at the beginning of this routine.
 
        if (allocmem) then
            allocate (sendbuffer(sendbuffersize), &
                      recvbuffer(recvbuffersize), stat=ierr)
            if (ierr /= 0) &
                call terminate("wallDistance", &
                               "Memory allocation failure for comm buffers")
        else
            call allocatetempmemory(.true.)
        end if
 
        ! Synchronize the processors.
 
        call mpi_barrier(adflow_comm_world, ierr)
 
        ! Write a message to stdout with the amount of time it
        ! took to compute the distances.
 
        !        if(myID == 0) then
        !          print "(*(A))", "# End wall distances level ", trim(integerString)
        !          print "(*(A, ES12.5))", "# Wall clock time:", mpi_wtime() - t0," sec."
        !          print "(a)", "#"
        !        endif
 
    end subroutine computewalldistance
 
    subroutine computenormalspacing(level, sps)
        !
        !       computeNormalSpacing computes the normal spacing of the first
        !       cell center from the viscous wall for the given multigrid
        !       level and spectral solution. This routine is called for
        !       turbulence models, which do not need the wall distance.
        !       However, they do need info of the first normal spacing for the
        !       monitoring of y+ and possibly for the boundary conditions.
        !       This is computed in this routine.
        !
        use constants
        use blockpointers, only: x, d2wall, nviscbocos, bcfaceid, bcdata, &
                                 nx, ny, nz, il, jl, kl, ndom
        use utils, only: setpointers
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: level, sps
        !
        !      Local variables.
        !
        integer(kind=intType) :: nn, mm, i, j
        integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd
 
        real(kind=realtype) :: nnx, nny, nnz, vecx, vecy, vecz, dot
 
        real(kind=realtype), dimension(:, :, :), pointer :: xface, xint
        real(kind=realtype), dimension(:, :), pointer :: dd2wall
 
        ! Loop over the domains.
 
        domain: do nn = 1, ndom
 
            ! Set the pointers for this block.
 
            call setpointers(nn, level, sps)
 
            ! Loop over the viscous subfaces of this block. Note that
            ! these are numbered first.
 
            bocos: do mm = 1, nviscbocos
 
                ! Set the pointers for the plane on the surface, one plane
                ! into the computational domain and the wall distance.
                ! This depends on the block face on which the subface is
                ! located. Note that the starting index of d2Wall is 2 and
                ! therefore a pointer offset will be needed later on.
 
                select case (bcfaceid(mm))
 
                case (imin)
                    xface => x(1, 1:, 1:, :); xint => x(2, 1:, 1:, :)
                    dd2wall => d2wall(2, :, :)
 
                case (imax)
                    xface => x(il, 1:, 1:, :); xint => x(nx, 1:, 1:, :)
                    dd2wall => d2wall(il, :, :)
 
                case (jmin)
                    xface => x(1:, 1, 1:, :); xint => x(1:, 2, 1:, :)
                    dd2wall => d2wall(:, 2, :)
 
                case (jmax)
                    xface => x(1:, jl, 1:, :); xint => x(1:, ny, 1:, :)
                    dd2wall => d2wall(:, jl, :)
 
                case (kmin)
                    xface => x(1:, 1:, 1, :); xint => x(1:, 1:, 2, :)
                    dd2wall => d2wall(:, :, 2)
 
                case (kmax)
                    xface => x(1:, 1:, kl, :); xint => x(1:, 1:, nz, :)
                    dd2wall => d2wall(:, :, kl)
 
                end select
 
                ! Store the face range of this subface a bit easier.
 
                jbeg = bcdata(mm)%jnBeg + 1; jend = bcdata(mm)%jnEnd
                ibeg = bcdata(mm)%inBeg + 1; iend = bcdata(mm)%inEnd
 
                ! Loop over the faces of the subfaces.
 
                do j = jbeg, jend
                    do i = ibeg, iend
 
                        ! Store the three components of the unit normal a
                        ! bit easier.
 
                        nnx = bcdata(mm)%norm(i, j, 1)
                        nny = bcdata(mm)%norm(i, j, 2)
                        nnz = bcdata(mm)%norm(i, j, 3)
 
                        ! Compute the vector from centroid of the adjacent cell
                        ! to the centroid of the face.
 
                        vecx = eighth * (xface(i - 1, j - 1, 1) + xface(i - 1, j, 1) &
                                         + xface(i, j - 1, 1) + xface(i, j, 1) &
                                         - xint(i - 1, j - 1, 1) - xint(i - 1, j, 1) &
                                         - xint(i, j - 1, 1) - xint(i, j, 1))
 
                        vecy = eighth * (xface(i - 1, j - 1, 2) + xface(i - 1, j, 2) &
                                         + xface(i, j - 1, 2) + xface(i, j, 2) &
                                         - xint(i - 1, j - 1, 2) - xint(i - 1, j, 2) &
                                         - xint(i, j - 1, 2) - xint(i, j, 2))
 
                        vecz = eighth * (xface(i - 1, j - 1, 3) + xface(i - 1, j, 3) &
                                         + xface(i, j - 1, 3) + xface(i, j, 3) &
                                         - xint(i - 1, j - 1, 3) - xint(i - 1, j, 3) &
                                         - xint(i, j - 1, 3) - xint(i, j, 3))
 
                        ! Compute the projection of this vector onto the normal
                        ! vector of the face. For a decent mesh there will not be
                        ! much of a difference between the projection and the
                        ! original mesh, but it does not hurt to do it.
 
                        dot = nnx * vecx + nny * vecy + nnz * vecz
 
                        ! As (nnx,nny,nnz) is a unit vector the distance to the
                        ! wall of the first cell center is given by the absolute
                        ! value of dot. Due to the use of pointers and the fact
                        ! that the original d2Wall array starts at 2 and offset
                        ! of -1 must be used to store the data at the correct
                        ! location.
 
                        dd2wall(i - 1, j - 1) = abs(dot)
 
                    end do
                end do
 
            end do bocos
 
        end do domain
 
    end subroutine computenormalspacing
 
    subroutine initwalldistance(level, sps, allocMem)
        !
        !       initWallDistance allocates the memory for the wall distance,
        !       if needed, and initializes the wall distance to a large value.
        !
        use constants
        use blockpointers, only: ndom, flowdoms
        use utils, only: terminate
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: level, sps
        logical, intent(in) :: allocMem
        !
        !      Local variables.
        !
        integer :: ierr
 
        integer(kind=intType) :: nn, il, jl, kl
 
        ! Loop over the domains.
 
        domain: do nn = 1, ndom
 
            ! Allocate the memory for d2Wall, if desired.
 
            if (allocmem) then
 
                il = flowdoms(nn, level, sps)%il
                jl = flowdoms(nn, level, sps)%jl
                kl = flowdoms(nn, level, sps)%kl
 
                allocate (flowdoms(nn, level, sps)%d2Wall(2:il, 2:jl, 2:kl), &
                          stat=ierr)
                if (ierr /= 0) &
                    call terminate("initWallDistance", &
                                   "Memory allocation failure for d2Wall")
            end if
 
            ! Initialize the wall distances to a large value.
 
            flowdoms(nn, level, sps)%d2Wall = large
 
        end do domain
 
    end subroutine initwalldistance
 
    subroutine determinedistance(level, sps)
        !
        !       determineDistance determines the distance from the center
        !       of the cell to the nearest viscous wall for owned cells.
        !
        use constants
        use adtapi, only: adtbuildsurfaceadt, adtmindistancesearch, adtdeallocateadts
        use blockpointers, only: x, flowdoms, kl, jl, il, ndom, nx, ny, nz, &
                                 sectionid, d2wall
        use communication, only: adflow_comm_world
        use section, only: sections
        use inputphysics, only: walloffset
        use utils, only: setpointers, terminate
        implicit none
        !
        !      Subroutine arguments
        !
        integer(kind=intType), intent(in) :: level, sps
        !
        !      Local parameter, which defines the name of the adt to be create
        !      and used here. Only needed because of the api of the adt library
        !
        character(len=10), parameter :: viscAdt = "ViscousADT"
        !
        !      Local variables.
        !
        integer :: ierr
 
        integer, dimension(:), allocatable :: procID
 
        integer(kind=intType) :: nCell, nCellPer, nTria
 
        integer(kind=intType), dimension(1, 1) :: connTria
        real(kind=realtype), dimension(3, 2) :: dummy
 
        integer(kind=intType), dimension(:), allocatable :: elementID
 
        real(kind=realtype), dimension(:), allocatable :: dist2
        real(kind=realtype), dimension(:), allocatable :: dist2per
        real(kind=realtype), dimension(:, :), allocatable :: coor, uvw
        real(kind=realtype), dimension(:, :), allocatable :: coorper
 
        integer(kind=adtElementType), dimension(:), allocatable :: &
            elementType
 
        integer(kind=intType) :: nn, mm, ll, ii, jj, i, j, k
 
        real(kind=realtype), dimension(3) :: xc
 
        ! Build the adt of the surface grid. As the api requires to
        ! specify both the quadrilateral and the triangular connectivity,
        ! some dummy variables must be passed.
 
        ntria = 0
        conntria = 0
        dummy = zero
 
        call adtbuildsurfaceadt(ntria, nquadvisc, nnodevisc, &
                                coorvisc, conntria, connvisc, &
                                dummy, .false., adflow_comm_world, &
                                viscadt)
 
        ! Determine the number of cell centers for which the distance
        ! must be computed. Also determine how many of them are periodic.
 
        ncell = 0
        ncellper = 0
 
        do nn = 1, ndom
            ll = flowdoms(nn, level, sps)%nx * flowdoms(nn, level, sps)%ny &
                 * flowdoms(nn, level, sps)%nz
            ncell = ncell + ll
 
            mm = flowdoms(nn, level, sps)%sectionID
            if (sections(mm)%periodic) ncellper = ncellper + ll
        end do
 
        ! Allocate the memory for the arrays needed by the ADT.
 
        allocate (coor(3, ncell), procid(ncell), elementtype(ncell), &
                  elementid(ncell), uvw(3, ncell), dist2(ncell), &
                  coorper(3, ncellper), dist2per(ncellper), stat=ierr)
        if (ierr /= 0) &
             call terminate("determineDistance", &
             "Memory allocation failure for the variables &
             &needed by the adt.")
        !
        !       Step 1: The search of the original coordinates; possibly a
        !               rotational periodic transformation is applied to align
        !               the sections.
        !
        ! Loop over the domains to store the coordinates of the cell
        ! centers of the owned cells. Apply the transformation such that
        ! the sections are aligned.
 
        mm = 0
        domains: do nn = 1, ndom
 
            ! Set the pointers for this block, store the section id a bit
            ! easier in ll and loop over the cell centers.
 
            call setpointers(nn, level, sps)
            ll = sectionid
 
            do k = 2, kl
                do j = 2, jl
                    do i = 2, il
 
                        ! Compute the coordinates of the cell center relative
                        ! to the rotation center of this section.
 
                        xc(1) = eighth * (x(i - 1, j - 1, k - 1, 1) + x(i, j - 1, k - 1, 1) &
                                          + x(i - 1, j, k - 1, 1) + x(i, j, k - 1, 1) &
                                          + x(i - 1, j - 1, k, 1) + x(i, j - 1, k, 1) &
                                          + x(i - 1, j, k, 1) + x(i, j, k, 1)) &
                                - sections(ll)%rotCenter(1)
 
                        xc(2) = eighth * (x(i - 1, j - 1, k - 1, 2) + x(i, j - 1, k - 1, 2) &
                                          + x(i - 1, j, k - 1, 2) + x(i, j, k - 1, 2) &
                                          + x(i - 1, j - 1, k, 2) + x(i, j - 1, k, 2) &
                                          + x(i - 1, j, k, 2) + x(i, j, k, 2)) &
                                - sections(ll)%rotCenter(2)
 
                        xc(3) = eighth * (x(i - 1, j - 1, k - 1, 3) + x(i, j - 1, k - 1, 3) &
                                          + x(i - 1, j, k - 1, 3) + x(i, j, k - 1, 3) &
                                          + x(i - 1, j - 1, k, 3) + x(i, j - 1, k, 3) &
                                          + x(i - 1, j, k, 3) + x(i, j, k, 3)) &
                                - sections(ll)%rotCenter(3)
 
                        ! Apply the periodic transformation for this section to
                        ! align it with other sections and store this coordinate
                        ! in the appropriate place in coor.
 
                        mm = mm + 1
                        coor(1, mm) = rotmatrixsections(ll, 1, 1) * xc(1) &
                                      + rotmatrixsections(ll, 1, 2) * xc(2) &
                                      + rotmatrixsections(ll, 1, 3) * xc(3) &
                                      + sections(ll)%rotCenter(1)
 
                        coor(2, mm) = rotmatrixsections(ll, 2, 1) * xc(1) &
                                      + rotmatrixsections(ll, 2, 2) * xc(2) &
                                      + rotmatrixsections(ll, 2, 3) * xc(3) &
                                      + sections(ll)%rotCenter(2)
 
                        coor(3, mm) = rotmatrixsections(ll, 3, 1) * xc(1) &
                                      + rotmatrixsections(ll, 3, 2) * xc(2) &
                                      + rotmatrixsections(ll, 3, 3) * xc(3) &
                                      + sections(ll)%rotCenter(3)
 
                        ! Initialize the distance squared, because this is an
                        ! inout argument in the call to adtMinDistanceSearch.
 
                        dist2(mm) = d2wall(i, j, k)
 
                    end do
                end do
            end do
        end do domains
 
        ! Perform the search. As no no interpolations are required,
        ! some dummies are passed.
 
        call adtmindistancesearch(ncell, coor, viscadt, &
                                  procid, elementtype, elementid, &
                                  uvw, dist2, 0_inttype, &
                                  dummy, dummy)
        !        if (myid == 0) then
        !           print *,'procID = '
        !           print *,procID
        !           print *,'elemID:',elementID
        !        end if
 
        !
        !       Step 2: For periodic sections the nearest wall may be in the
        !               periodic part of the grid that is not stored.
        !               Therefore apply the periodic transformation to the
        !               node and compute the minimum distance for this
        !               coordinate.
        !
        ! Initialize the counters mm and ii. Mm is the counter for coor;
        ! ii is the counter for coorPer.
 
        mm = 0
        ii = 0
 
        ! Loop over the domains and find the periodic ones.
 
        domainsper1: do nn = 1, ndom
            jj = flowdoms(nn, level, sps)%nx * flowdoms(nn, level, sps)%ny &
                 * flowdoms(nn, level, sps)%nz
 
            ll = flowdoms(nn, level, sps)%sectionID
 
            ! Check if the section is periodic.
 
            if (sections(ll)%periodic) then
 
                ! Loop over the corresponding entries in coor of this block
                ! and apply the periodic transformation. The transformed
                ! coordinates are stored in coorPer. Also initialize
                ! the wall distance squared to the value just computed.
 
                do i = 1, jj
                    mm = mm + 1
                    ii = ii + 1
 
                    xc(1) = coor(1, mm) - sections(ll)%rotCenter(1)
                    xc(2) = coor(2, mm) - sections(ll)%rotCenter(2)
                    xc(3) = coor(3, mm) - sections(ll)%rotCenter(3)
 
                    coorper(1, ii) = sections(ll)%rotMatrix(1, 1) * xc(1) &
                                     + sections(ll)%rotMatrix(1, 2) * xc(2) &
                                     + sections(ll)%rotMatrix(1, 3) * xc(3) &
                                     + sections(ll)%rotCenter(1) &
                                     + sections(ll)%translation(1)
 
                    coorper(2, ii) = sections(ll)%rotMatrix(2, 1) * xc(1) &
                                     + sections(ll)%rotMatrix(2, 2) * xc(2) &
                                     + sections(ll)%rotMatrix(2, 3) * xc(3) &
                                     + sections(ll)%rotCenter(2) &
                                     + sections(ll)%translation(2)
 
                    coorper(3, ii) = sections(ll)%rotMatrix(3, 1) * xc(1) &
                                     + sections(ll)%rotMatrix(3, 2) * xc(2) &
                                     + sections(ll)%rotMatrix(3, 3) * xc(3) &
                                     + sections(ll)%rotCenter(3) &
                                     + sections(ll)%translation(3)
 
                    dist2per(ii) = dist2(mm)
                end do
 
            else
                ! Section is not periodic. Update the counter mm.
 
                mm = mm + jj
            end if
 
        end do domainsper1
 
        ! Perform the adt search of this set of periodic coordinates.
        ! As no no interpolations are required, some dummies are passed.
 
        call adtmindistancesearch(ncellper, coorper, viscadt, &
                                  procid, elementtype, elementid, &
                                  uvw, dist2per, 0_inttype, &
                                  dummy, dummy)
        !
        !       Step 3: Also apply the inverse periodic transformation.
        !
        ! Initialize the counters mm and ii. Mm is the counter for coor;
        ! ii is the counter for coorPer.
 
        mm = 0
        ii = 0
 
        ! Loop over the domains and find the periodic ones.
 
        domainsper2: do nn = 1, ndom
            jj = flowdoms(nn, level, sps)%nx * flowdoms(nn, level, sps)%ny &
                 * flowdoms(nn, level, sps)%nz
 
            ll = flowdoms(nn, level, sps)%sectionID
 
            ! Check if the section is periodic.
 
            if (sections(ll)%periodic) then
 
                ! Loop over the corresponding entries in coor of this block
                ! and apply the inverse periodic transformation. Again the
                ! transformed coordinates are stored in coorPer. Note that
                ! the inverse of the rotation matrix is the transpose and
                ! that the translation vector should be multiplied by the
                ! inverse of the rotation matrix. Note that the wall distance
                ! has already been initialized in the previous periodic
                ! search.
 
                do i = 1, jj
                    mm = mm + 1
                    ii = ii + 1
 
                    xc(1) = coor(1, mm) - sections(ll)%rotCenter(1) &
                            - sections(ll)%translation(1)
                    xc(2) = coor(2, mm) - sections(ll)%rotCenter(2) &
                            - sections(ll)%translation(2)
                    xc(3) = coor(3, mm) - sections(ll)%rotCenter(3) &
                            - sections(ll)%translation(3)
 
                    coorper(1, ii) = sections(ll)%rotMatrix(1, 1) * xc(1) &
                                     + sections(ll)%rotMatrix(2, 1) * xc(2) &
                                     + sections(ll)%rotMatrix(3, 1) * xc(3) &
                                     + sections(ll)%rotCenter(1)
 
                    coorper(2, ii) = sections(ll)%rotMatrix(1, 2) * xc(1) &
                                     + sections(ll)%rotMatrix(2, 2) * xc(2) &
                                     + sections(ll)%rotMatrix(3, 2) * xc(3) &
                                     + sections(ll)%rotCenter(2)
 
                    coorper(3, ii) = sections(ll)%rotMatrix(1, 3) * xc(1) &
                                     + sections(ll)%rotMatrix(2, 3) * xc(2) &
                                     + sections(ll)%rotMatrix(3, 3) * xc(3) &
                                     + sections(ll)%rotCenter(3)
                end do
 
            else
                ! Section is not periodic. Update the counter mm.
 
                mm = mm + jj
            end if
 
        end do domainsper2
 
        ! Perform the adt search of this set of periodic coordinates.
        ! As no no interpolations are required, some dummies are passed.
 
        call adtmindistancesearch(ncellper, coorper, viscadt, &
                                  procid, elementtype, elementid, &
                                  uvw, dist2per, 0_inttype, &
                                  dummy, dummy)
        !
        !       Step 4: Store the minimum distance in the block type.
        !
        mm = 0
        ii = 0
 
        domainsstore: do nn = 1, ndom
 
            ! Set the pointers for this block and store the section id a
            ! bit easier in ll.
 
            call setpointers(nn, level, sps)
            ll = sectionid
 
            ! Check if the section is periodic. If so the distance is set
            ! to the minimum of the value computed in step 1 and the
            ! periodic values. Note that mm should not be updated in this
            ! loop, because it is updated in the loop over the cell centers
            ! of this block. Instead the counter j is used.
 
            if (sections(ll)%periodic) then
 
                jj = nx * ny * nz
 
                j = mm
                do i = 1, jj
                    j = j + 1
                    ii = ii + 1
 
                    dist2(j) = min(dist2(j), dist2per(ii))
                end do
            end if
 
            ! Loop over the cell centers of the block to store the wall
            ! distance. Note that dist2 stores the distance squared and
            ! thus a square root must be taken.
            ! Add a possible offset for the debugging of wall functions.
 
            do k = 2, kl
                do j = 2, jl
                    do i = 2, il
                        mm = mm + 1
                        d2wall(i, j, k) = sqrt(dist2(mm)) + walloffset
                    end do
                end do
            end do
 
        end do domainsstore
 
        ! Release the memory of the ADT and the arrays of the module
        ! viscSurface.
 
        call adtdeallocateadts(viscadt)
 
        deallocate (connvisc, coorvisc, rotmatrixsections, stat=ierr)
        if (ierr /= 0) &
             call terminate("determineDistance", &
             "Deallocation error for the arrays &
             &of viscSurface")
 
        ! Release the variables needed by the ADT.
 
        deallocate (coor, procid, elementtype, elementid, uvw, dist2, &
                    coorper, dist2per, stat=ierr)
        if (ierr /= 0) &
             call terminate("determineDistance", &
             "Deallocation failure for the variables &
             &needed by the adt.")
 
    end subroutine determinedistance
 
    subroutine localviscoussurfacemesh(multSections, level, sps)
        !
        !       localViscousSurfaceMesh stores the local viscous surface
        !       mesh (with possible periodic extensions in conn and coor.
        !
        use constants
        use blockpointers, only: bcdata, x, il, jl, kl, bcfaceid, sectionid, &
                                 flowdoms, nbocos, ndom, bctype
        use communication, only: myid, adflow_comm_world
        use section, only: sections, nsections
        use utils, only: setpointers, terminate
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: level, sps
        integer(kind=intType), dimension(*), intent(in) :: multSections
        !
        !      Local variables.
        !
        integer :: size, ierr
 
        integer(kind=intType) :: nn, mm, i, j, k
        integer(kind=intType) :: np, nq, npOld, np1, nqOld, mp, mq
        integer(kind=intType) :: nq1, nq2, nq3, nq4, sec, row, col
        integer(kind=intType) :: iBeg, jBeg, iEnd, jEnd
 
        integer(kind=intType), dimension(3) :: ind
 
        real(kind=realtype) :: length, dot, xx, yy, zz, r1, r2, aaa, bbb
        real(kind=realtype) :: theta, costheta, sintheta
 
        real(kind=realtype), dimension(3, 3) :: a
        real(kind=realtype), dimension(nSections) :: thetanmin, &
                                                     thetanmax, &
                                                     thetapmin, &
                                                     thetapmax, tmp
        real(kind=realtype), dimension(nSections, 3) :: rad1, rad2, axis
 
        real(kind=realtype), dimension(:, :, :), pointer :: xface
 
        !       Determine the unit vectors of the local coordinate system
        !       aligned with the rotation axis of the possible rotational
        !       periodic section.
        !
        do nn = 1, nsections
            if (sections(nn)%nSlices == 1) cycle
 
            ! Section is rotational periodic. First determine the rotation
            ! axis. This is the eigenvector which corresponds to the
            ! eigenvalue 1 of the transformation matrix.
 
            ! Store rot - i in a and initialize ind.
 
            ind(1) = 1; ind(2) = 2; ind(3) = 3
 
            a(1, 1) = sections(nn)%rotMatrix(1, 1) - one
            a(1, 2) = sections(nn)%rotMatrix(1, 2)
            a(1, 3) = sections(nn)%rotMatrix(1, 3)
 
            a(2, 1) = sections(nn)%rotMatrix(2, 1)
            a(2, 2) = sections(nn)%rotMatrix(2, 2) - one
            a(2, 3) = sections(nn)%rotMatrix(2, 3)
 
            a(3, 1) = sections(nn)%rotMatrix(3, 1)
            a(3, 2) = sections(nn)%rotMatrix(3, 2)
            a(3, 3) = sections(nn)%rotMatrix(3, 3) - one
 
            ! Loop over the two times that Gaussian elimination must be
            ! applied.
 
            loopgauss: do k = 1, 2
 
                ! Find the largest value in the sub-matrix.
 
                aaa = abs(a(k, k)); row = k; col = k
                do j = k, 3
                    do i = k, 3
                        bbb = abs(a(i, j))
                        if (bbb > aaa) then
                            aaa = bbb
                            row = i
                            col = j
                        end if
                    end do
                end do
 
                ! Swap the rows k and row.
 
                do j = 1, 3
                    aaa = a(k, j)
                    a(k, j) = a(row, j)
                    a(row, j) = aaa
                end do
 
                ! Swap the colums k and col; also swap ind(k) and ind(col).
 
                i = ind(k)
                ind(k) = ind(col)
                ind(col) = i
                do i = 1, 3
                    aaa = a(i, k)
                    a(i, k) = a(i, col)
                    a(i, col) = aaa
                end do
 
                ! Perform gaussian eliMination, because now it's sure that
                ! the element (k,k) is non-zero.
 
                aaa = one / a(k, k)
                do i = (k + 1), 3
                    bbb = a(i, k) * aaa
                    do j = k, 3
                        a(i, j) = a(i, j) - bbb * a(k, j)
                    end do
                end do
 
            end do loopgauss
 
            ! Due to the full pivoting it is now guaranteed that the elements
            ! a(ind(1),ind(1)) and a(ind(2),ind(2)) are nonzero and
            ! a(ind(3),ind(3)) == zero. Remember that the rotation matrix
            ! only has 1 eigenvalue of one. Set axis(ind(3)) to one and
            ! determine the other two elements of the eigen vector.
 
            axis(nn, ind(3)) = one
            axis(nn, ind(2)) = -(a(2, 3) * axis(nn, ind(3))) / a(2, 2)
            axis(nn, ind(1)) = -(a(1, 3) * axis(nn, ind(3)) &
                                 + a(1, 2) * axis(nn, ind(2))) / a(1, 1)
 
            ! Create a unit vector.
 
            length = one / sqrt(axis(nn, 1)**2 + axis(nn, 2)**2 + axis(nn, 3)**2)
            axis(nn, 1) = axis(nn, 1) * length
            axis(nn, 2) = axis(nn, 2) * length
            axis(nn, 3) = axis(nn, 3) * length
 
            ! Make sure that the largest component of this vector is
            ! positive, such that a unique definition of the rotation
            ! axis is obtained. Use dot and length as a temporary
            ! storage.
 
            dot = axis(nn, 1); length = abs(dot)
            if (abs(axis(nn, 2)) > length) then
                dot = axis(nn, 2); length = abs(dot)
            end if
            if (abs(axis(nn, 3)) > length) then
                dot = axis(nn, 3); length = abs(dot)
            end if
 
            if (dot < zero) then
                axis(nn, 1) = -axis(nn, 1)
                axis(nn, 2) = -axis(nn, 2)
                axis(nn, 3) = -axis(nn, 3)
            end if
 
            ! Determine the two vectors which determine the plane normal
            ! to the axis of rotation.
 
            ! Initial guess of rad1. First try the y-axis. If not good
            ! enough try the z-axis.
 
            if (abs(axis(nn, 2)) < 0.707107_realtype) then
                rad1(nn, 1) = zero
                rad1(nn, 2) = one
                rad1(nn, 3) = zero
            else
                rad1(nn, 1) = zero
                rad1(nn, 2) = zero
                rad1(nn, 3) = one
            end if
 
            ! Make sure that rad1 is normal to axis. Create a unit
            ! vector again.
 
            dot = rad1(nn, 1) * axis(nn, 1) + rad1(nn, 2) * axis(nn, 2) &
                  + rad1(nn, 3) * axis(nn, 3)
            rad1(nn, 1) = rad1(nn, 1) - dot * axis(nn, 1)
            rad1(nn, 2) = rad1(nn, 2) - dot * axis(nn, 2)
            rad1(nn, 3) = rad1(nn, 3) - dot * axis(nn, 3)
 
            length = one / (rad1(nn, 1)**2 + rad1(nn, 2)**2 + rad1(nn, 3)**2)
            rad1(nn, 1) = rad1(nn, 1) * length
            rad1(nn, 2) = rad1(nn, 2) * length
            rad1(nn, 3) = rad1(nn, 3) * length
 
            ! Create the second vector which spans the radIal plane. This
            ! must be normal to both axis and rad1, i.e. the cross-product.
 
            rad2(nn, 1) = axis(nn, 2) * rad1(nn, 3) - axis(nn, 3) * rad1(nn, 2)
            rad2(nn, 2) = axis(nn, 3) * rad1(nn, 1) - axis(nn, 1) * rad1(nn, 3)
            rad2(nn, 3) = axis(nn, 1) * rad1(nn, 2) - axis(nn, 2) * rad1(nn, 1)
 
        end do
 
        ! Initialize the values of thetaNMin, etc.
 
        thetanmin = zero
        thetanmax = -pi
        thetapmin = pi
        thetapmax = zero
        !
        !       Determine the local values of thetaNMin, etc. for the
        !       different sections.
        !
        do nn = 1, ndom
 
            ! Store the section id of the block a bit easier. Continue with
            ! the next block if this section consist of only one slice.
 
            sec = flowdoms(nn, level, sps)%sectionID
            if (sections(sec)%nSlices == 1) cycle
 
            ! Set the pointers for this block.
 
            call setpointers(nn, level, sps)
 
            ! Initialize nq1, nq2, nq3 and nq4 to 0. These integers store
            ! the number of nodes in the first, second, third and fourth
            ! quadrant respectivily.
 
            nq1 = 0; nq2 = 0; nq3 = 0; nq4 = 0
 
            ! Loop over the nodes of this block.
 
            do k = 1, kl
                do j = 1, jl
                    do i = 1, il
 
                        ! Determine the coordinates relative to the
                        ! center of rotation.
 
                        xx = x(i, j, k, 1) - sections(sec)%rotCenter(1)
                        yy = x(i, j, k, 2) - sections(sec)%rotCenter(2)
                        zz = x(i, j, k, 3) - sections(sec)%rotCenter(3)
 
                        ! Determine the radIal components in the local
                        ! cylindrical coordinate system of the section.
 
                        r1 = xx * rad1(sec, 1) + yy * rad1(sec, 2) + zz * rad1(sec, 3)
                        r2 = xx * rad2(sec, 1) + yy * rad2(sec, 2) + zz * rad2(sec, 3)
 
                        ! Determine the angle if r1 or r2 is nonzero.
 
                        if ((abs(r1) >= eps) .or. (abs(r2) >= eps)) then
 
                            theta = atan2(r2, r1)
 
                            ! Update the minimum and maximum angle for this
                            ! section, depending on the sign of theta.
 
                            if (theta >= zero) then
                                thetapmin(sec) = min(thetapmin(sec), theta)
                                thetapmax(sec) = max(thetapmax(sec), theta)
                            end if
 
                            if (theta <= zero) then
                                thetanmin(sec) = min(thetanmin(sec), theta)
                                thetanmax(sec) = max(thetanmax(sec), theta)
                            end if
 
                            ! Determine the quadrant in which this node is located
                            ! and update the corresponding counter.
 
                            if (theta <= -half * pi) then
                                nq3 = nq3 + 1
                            else if (theta <= zero) then
                                nq4 = nq4 + 1
                            else if (theta <= half * pi) then
                                nq1 = nq1 + 1
                            else
                                nq2 = nq2 + 1
                            end if
 
                        end if
 
                    end do
                end do
            end do
 
            ! Modify the minimum and maximum angles if nodes are present
            ! in multiple quadrants.
 
            if (nq1 > 0 .and. nq4 > 0) then
 
                ! Nodes in both the 1st and 4th quadrant. Update the
                ! corresponding minimum and maximum angle.
 
                thetanmax(sec) = zero
                thetapmin(sec) = zero
 
            end if
 
            if (nq2 > 0 .and. nq3 > 0) then
 
                ! Nodes in both the 2nd and 3rd quadrant. Update the
                ! corresponding minimum and maximum angle.
 
                thetanmin(sec) = -pi
                thetapmax(sec) = pi
 
            end if
 
        end do
 
        ! Determine the minimum of the minimum angles and the maximum of
        ! the maximum angles for all sections.
 
        size = nsections
        call mpi_allreduce(thetanmax, tmp, size, adflow_real, mpi_max, &
                           adflow_comm_world, ierr)
        thetanmax = tmp
 
        call mpi_allreduce(thetapmax, tmp, size, adflow_real, mpi_max, &
                           adflow_comm_world, ierr)
        thetapmax = tmp
 
        call mpi_allreduce(thetanmin, tmp, size, adflow_real, mpi_min, &
                           adflow_comm_world, ierr)
        thetanmin = tmp
 
        call mpi_allreduce(thetapmin, tmp, size, adflow_real, mpi_min, &
                           adflow_comm_world, ierr)
        thetapmin = tmp
 
        ! Allocate the memory for rotMatrixSections, the rotation
        ! matrices of the sections needed for the alignment.
 
        allocate (rotmatrixsections(nsections, 3, 3), stat=ierr)
        if (ierr /= 0) &
             call terminate("localViscousSurfaceMesh", &
             "Memory allocation failure for &
             &rotMatrixSections")
        !
        !       Determine the rotation matrix for each section, which aligns
        !       the rotational periodic sections with other sections.
        !
        do nn = 1, nsections
 
            ! Test if a rotation is actually needed.
 
            testrot: if (sections(nn)%nSlices == 1 .or. &
                         thetapmin(nn) == zero) then
 
                ! Section consist out of 1 slice or the slice crosses the
                ! line theta == 0. For both cases the rotation matrix is
                ! the identity matrix.
 
                rotmatrixsections(nn, 1, 1) = one
                rotmatrixsections(nn, 1, 2) = zero
                rotmatrixsections(nn, 1, 3) = zero
 
                rotmatrixsections(nn, 2, 1) = zero
                rotmatrixsections(nn, 2, 2) = one
                rotmatrixsections(nn, 2, 3) = zero
 
                rotmatrixsections(nn, 3, 1) = zero
                rotmatrixsections(nn, 3, 2) = zero
                rotmatrixsections(nn, 3, 3) = one
 
                else testrot
 
                ! Section consist out of multiple slices and the current slice
                ! does not cross the line theta == 0. The rotation matrix
                ! for alignment must be computed.
 
                theta = two * pi / sections(nn)%nSlices
 
                ! Determine the number of rotations needed to align the mesh.
 
                if (thetanmin(nn) < zero) then
 
                    ! The section lies (at least partially) in the third and
                    ! fourth quadrant. Determine the number of rotations for
                    ! alignment; this is a positive number.
 
                    mm = -thetanmax(nn) / theta + 1
 
                else
 
                    ! The section lies (completely) in the first and second
                    ! quadrant. The number of rotations will be a negative
                    ! number now.
 
                    mm = -thetapmin(nn) / theta - 1
 
                end if
 
                ! Compute the rotation angle in the local cylindrical frame
                ! and its sine and cosine.
 
                theta = mm * theta
                costheta = cos(theta)
                sintheta = sin(theta)
 
                ! Apply the transformation to obtain the matrix in the
                ! original cartesian frame.
 
                rotmatrixsections(nn, 1, 1) = axis(nn, 1) * axis(nn, 1) &
                                              + costheta * (rad1(nn, 1) * rad1(nn, 1) + rad2(nn, 1) * rad2(nn, 1))
                rotmatrixsections(nn, 1, 2) = axis(nn, 1) * axis(nn, 2) &
                                              + costheta * (rad1(nn, 1) * rad1(nn, 2) + rad2(nn, 1) * rad2(nn, 2)) &
                                              + sintheta * (rad1(nn, 2) * rad2(nn, 1) - rad1(nn, 1) * rad2(nn, 2))
                rotmatrixsections(nn, 1, 3) = axis(nn, 1) * axis(nn, 3) &
                                              + costheta * (rad1(nn, 1) * rad1(nn, 3) + rad2(nn, 1) * rad2(nn, 3)) &
                                              + sintheta * (rad1(nn, 3) * rad2(nn, 1) - rad1(nn, 1) * rad2(nn, 3))
 
                rotmatrixsections(nn, 2, 1) = axis(nn, 1) * axis(nn, 2) &
                                              + costheta * (rad1(nn, 1) * rad1(nn, 2) + rad2(nn, 1) * rad2(nn, 2)) &
                                              - sintheta * (rad1(nn, 2) * rad2(nn, 1) - rad1(nn, 1) * rad2(nn, 2))
                rotmatrixsections(nn, 2, 2) = axis(nn, 2) * axis(nn, 2) &
                                              + costheta * (rad1(nn, 2) * rad1(nn, 2) + rad2(nn, 2) * rad2(nn, 2))
                rotmatrixsections(nn, 2, 3) = axis(nn, 2) * axis(nn, 3) &
                                              + costheta * (rad1(nn, 2) * rad1(nn, 3) + rad2(nn, 2) * rad2(nn, 3)) &
                                              + sintheta * (rad1(nn, 3) * rad2(nn, 2) - rad1(nn, 2) * rad2(nn, 3))
 
                rotmatrixsections(nn, 3, 1) = axis(nn, 1) * axis(nn, 3) &
                                              + costheta * (rad1(nn, 1) * rad1(nn, 3) + rad2(nn, 1) * rad2(nn, 3)) &
                                              - sintheta * (rad1(nn, 3) * rad2(nn, 1) - rad1(nn, 1) * rad2(nn, 3))
                rotmatrixsections(nn, 3, 2) = axis(nn, 2) * axis(nn, 3) &
                                              + costheta * (rad1(nn, 2) * rad1(nn, 3) + rad2(nn, 2) * rad2(nn, 3)) &
                                              - sintheta * (rad1(nn, 3) * rad2(nn, 2) - rad1(nn, 2) * rad2(nn, 3))
                rotmatrixsections(nn, 3, 3) = axis(nn, 3) * axis(nn, 3) &
                                              + costheta * (rad1(nn, 3) * rad1(nn, 3) + rad2(nn, 3) * rad2(nn, 3))
 
            end if testrot
 
        end do
        !
        !       Determine the local viscous surface grid.
        !
        np = 0
        nq = 0
 
        loopdomains: do nn = 1, ndom
 
            ! Set the pointers for this block and store the section id
            ! a bit easier.
 
            call setpointers(nn, level, sps)
            sec = sectionid
 
            ! Loop over the subfaces of this block and test if this is
            ! a viscous subface.
 
            loopbocos: do mm = 1, nbocos
                testviscous: if (bctype(mm) == nswalladiabatic .or. &
                                 bctype(mm) == nswallisothermal) then
 
                    ! Viscous subface. Set the pointer for the coordinates of
                    ! the face.
 
                    select case (bcfaceid(mm))
 
                    case (imin)
                        xface => x(1, 1:, 1:, :)
 
                    case (imax)
                        xface => x(il, 1:, 1:, :)
 
                    case (jmin)
                        xface => x(1:, 1, 1:, :)
 
                    case (jmax)
                        xface => x(1:, jl, 1:, :)
 
                    case (kmin)
                        xface => x(1:, 1:, 1, :)
 
                    case (kmax)
                        xface => x(1:, 1:, kl, :)
 
                    end select
 
                    ! Store the nodal range of this subface a bit easier.
 
                    jbeg = bcdata(mm)%jnBeg; jend = bcdata(mm)%jnEnd
                    ibeg = bcdata(mm)%inBeg; iend = bcdata(mm)%inEnd
 
                    ! Store the old value of the number of points stored and
                    ! determine the new coordinates.
 
                    npold = np
                    do j = jbeg, jend
                        do i = ibeg, iend
 
                            ! Determine the coordinates relative to the rotation
                            ! center of this section.
 
                            xx = xface(i, j, 1) - sections(sec)%rotCenter(1)
                            yy = xface(i, j, 2) - sections(sec)%rotCenter(2)
                            zz = xface(i, j, 3) - sections(sec)%rotCenter(3)
 
                            ! Update the counter and determine the surface mesh
                            ! coordinates.
 
                            np = np + 1
                            coorvisc(1, np) = rotmatrixsections(sec, 1, 1) * xx &
                                              + rotmatrixsections(sec, 1, 2) * yy &
                                              + rotmatrixsections(sec, 1, 3) * zz &
                                              + sections(sec)%rotCenter(1)
 
                            coorvisc(2, np) = rotmatrixsections(sec, 2, 1) * xx &
                                              + rotmatrixsections(sec, 2, 2) * yy &
                                              + rotmatrixsections(sec, 2, 3) * zz &
                                              + sections(sec)%rotCenter(2)
 
                            coorvisc(3, np) = rotmatrixsections(sec, 3, 1) * xx &
                                              + rotmatrixsections(sec, 3, 2) * yy &
                                              + rotmatrixsections(sec, 3, 3) * zz &
                                              + sections(sec)%rotCenter(3)
                        end do
                    end do
 
                    ! Determine and store the connectivity of this subface.
 
                    np1 = iend - ibeg + 1
                    nqold = nq
 
                    do j = (jbeg + 1), jend
                        do i = (ibeg + 1), iend
 
                            ! Update the counter nq and determine the 4 indices
                            ! of the surface quad.
 
                            nq = nq + 1
 
                            connvisc(1, nq) = npold + (j - jbeg - 1) * np1 + i - ibeg
                            connvisc(2, nq) = connvisc(1, nq) + 1
                            connvisc(3, nq) = connvisc(2, nq) + np1
                            connvisc(4, nq) = connvisc(3, nq) - 1
 
                        end do
                    end do
 
                    ! Loop over the number of times the subface must be stored.
                    ! This happens when the rotational periodicity differs from
                    ! section to section. Note that this loop starts at k == 2.
 
                    loopmultiplicity: do k = 2, multsections(sec)
 
                        ! Store the current number of nodes and quads
                        ! in mp and mq respectivily.
 
                        mp = np
                        mq = nq
 
                        ! Loop over the of points on this subface.
 
                        do i = (npold + 1), mp
 
                            ! Determine the coordinates relative to the center
                            ! of rotation.
 
                            np = np + 1
 
                            xx = coorvisc(1, i) - sections(sec)%rotCenter(1)
                            yy = coorvisc(2, i) - sections(sec)%rotCenter(2)
                            zz = coorvisc(3, i) - sections(sec)%rotCenter(3)
 
                            ! Update the counter np and determine the new
                            ! coordinates after the transformation.
 
                            coorvisc(1, np) = sections(sec)%rotMatrix(1, 1) * xx &
                                              + sections(sec)%rotMatrix(1, 2) * yy &
                                              + sections(sec)%rotMatrix(1, 3) * zz &
                                              + sections(sec)%rotCenter(1) &
                                              + sections(sec)%translation(1)
 
                            coorvisc(2, np) = sections(sec)%rotMatrix(2, 1) * xx &
                                              + sections(sec)%rotMatrix(2, 2) * yy &
                                              + sections(sec)%rotMatrix(2, 3) * zz &
                                              + sections(sec)%rotCenter(2) &
                                              + sections(sec)%translation(2)
 
                            coorvisc(3, np) = sections(sec)%rotMatrix(3, 1) * xx &
                                              + sections(sec)%rotMatrix(3, 2) * yy &
                                              + sections(sec)%rotMatrix(3, 3) * zz &
                                              + sections(sec)%rotCenter(3) &
                                              + sections(sec)%translation(3)
                        end do
 
                        ! Store the number of nodes in this subface in j
                        ! and determine the connectivity of this rotated part.
 
                        j = np - mp
                        do i = (nqold + 1), mq
 
                            ! Update the counter nq and set the new connectivity,
                            ! which is the old connectivity plus an offset.
 
                            nq = nq + 1
                            connvisc(1, nq) = connvisc(1, i) + j
                            connvisc(2, nq) = connvisc(2, i) + j
                            connvisc(3, nq) = connvisc(3, i) + j
                            connvisc(4, nq) = connvisc(4, i) + j
 
                        end do
 
                        ! Copy the values of mp and mq into npOld and nqOld for
                        ! the next multiple of the slice. Idem for np in mp, etc.
 
                        npold = mp
                        nqold = mq
 
                        mp = np
                        mq = nq
 
                    end do loopmultiplicity
 
                end if testviscous
            end do loopbocos
        end do loopdomains
 
    end subroutine localviscoussurfacemesh
 
    subroutine updatewalldistancealllevels
        !
        !       updateWallDistanceAllLevels updates the wall distances for
        !       the cell centers on all grid levels. This routine is typically
        !       called when grid parts have been moved, either due to a
        !       physical motion of some parts or due to deformation.
        !
        use constants
        use block, only: flowdoms
        use inputphysics, only: equations
        use iteration, only: groundlevel
        implicit none
        !
        !      Local variables.
        !
        integer(kind=intType) :: nLevels, nn
 
        ! Return immediately if the rans equations are not solved.
 
        if (equations /= ransequations) return
 
        ! Loop over the grid levels and call wallDistance.
 
        nlevels = ubound(flowdoms, 2)
        do nn = groundlevel, nlevels
            call computewalldistance(nn, .false.)
        end do
 
    end subroutine updatewalldistancealllevels
 
    subroutine viscoussurfacemesh(level, sps)
        !
        !       viscousSurfaceMesh determines and stores the entire viscous
        !       surface possibly extended by periodic parts.
        !
        use constants
        use block, only: flowdoms, ndom
        use communication, only: adflow_comm_world, myid
        use section, only: nsections, sections
        use utils, only: terminate
        implicit none
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: level, sps
        !
        !      Local variables.
        !
        integer :: ierr
 
        integer(kind=intType) :: nn, mm, ii
        integer(kind=intType) :: ni, nj, nk
 
        integer(kind=intType), dimension(nSections) :: multSections
 
        ! Determine the minimum number of slices present in a certain
        ! section. For time accurate computations the number of slices
        ! is identical for all sections, but for steady flow using the
        ! mixing plane assumption this is not necessarily the case.
 
        mm = sections(1)%nSlices
        do nn = 2, nsections
            mm = min(mm, sections(nn)%nSlices)
        end do
 
        ! Determine the multiplicity of every section needed in the
        ! surface mesh, such that every part covers an angle which is
        ! at least equal to the angle of the largest section. Again note
        ! that this multiplicity is 1 for all sections if a time
        ! accurate computation is performed.
 
        do nn = 1, nsections
            multsections(nn) = sections(nn)%nSlices / mm
            if (sections(nn)%nSlices > mm * multsections(nn)) &
                multsections(nn) = multsections(nn) + 1
        end do
 
        ! Determine the local number of viscous nodes and quads.
        ! Note that these numbers are identical for all spectral
        ! solutions and thus it is okay to take the 1st one.
 
        nnodevisc = 0
        nquadvisc = 0
 
        do nn = 1, ndom
            do mm = 1, flowdoms(nn, level, 1)%nBocos
                if (flowdoms(nn, level, 1)%BCType(mm) == nswalladiabatic .or. &
                    flowdoms(nn, level, 1)%BCType(mm) == nswallisothermal) then
 
                    ! Determine the number of nodes of the subface in the
                    ! three directions.
 
                    ni = flowdoms(nn, level, 1)%inEnd(mm) &
                         - flowdoms(nn, level, 1)%inBeg(mm)
                    nj = flowdoms(nn, level, 1)%jnEnd(mm) &
                         - flowdoms(nn, level, 1)%jnBeg(mm)
                    nk = flowdoms(nn, level, 1)%knEnd(mm) &
                         - flowdoms(nn, level, 1)%knBeg(mm)
 
                    ! Determine the multiplication factor, because of the
                    ! possible multiple sections.
 
                    ii = flowdoms(nn, level, 1)%sectionId
                    ii = multsections(ii)
 
                    ! Update the number of nodes and quads. Take the
                    ! multiplicity into account.
 
                    nnodevisc = nnodevisc + ii * (ni + 1) * (nj + 1) * (nk + 1)
                    nquadvisc = nquadvisc + ii * max(ni, 1_inttype) &
                                * max(nj, 1_inttype) &
                                * max(nk, 1_inttype)
                end if
            end do
        end do
 
        ! Determine the global number of elements on the viscous
        ! surfaces. Return if there are no viscous quads present.
 
        call mpi_allreduce(nquadvisc, nquadviscglob, 1, adflow_integer, &
                           mpi_sum, adflow_comm_world, ierr)
 
        if (nquadviscglob == 0) return
 
        ! Allocate the memory for the local connectivity and coordinates.
 
        allocate (connvisc(4, nquadvisc), coorvisc(3, nnodevisc), &
                  stat=ierr)
        if (ierr /= 0) &
             call terminate("viscousSurfaceMesh", &
             "Memory allocation failure for connVisc &
             &and coorVisc.")
 
        ! Determine the local viscous surface mesh, possibly rotated
        ! to align the other sections.
 
        call localviscoussurfacemesh(multsections, level, sps)
 
    end subroutine viscoussurfacemesh
 
    subroutine determinewallassociation(level, sps)
 
        ! This routine will determine the closest surface point for every
        ! field cell. Special treatment is required for overlapping surfaces.
 
        use constants
        use adtapi, only: mindistancesearch
        use adtdata, only: adtbboxtargettype
        use adtlocalsearch, only: mindistancetreesearchsinglepoint
        use adtutils, only: stack
        use adtbuild, only: buildserialquad, destroyserialquad
        use blockpointers
        use communication
        use inputphysics
        use inputtimespectral
        use oversetdata, only: oversetpresent, oversetwall, nclusters, clusters, cumdomproc
        use inputoverset
        use adjointvars
        use surfacefamilies, only: bcfamgroups
        use utils, only: setpointers, echk
        use sorting, only: unique
        implicit none
 
        ! Input Variables
        integer(kind=intType), intent(in) :: level, sps
 
        ! Local Variables
        integer(kind=intType) :: i, j, k, l, ii, jj, kk, nn, mm, iNode, iCell, c
        integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, ni, nj, nUnique, cellID, cellID2
        integer(kind=intType) :: ierr, iDim
 
        ! Data for local surface
        integer(kind=intType) :: nNodes, nCells
        logical :: gridHasOverset
 
        ! Overset Walls for storing the surface ADT's
        type(oversetwall), dimension(:), allocatable, target :: walls
        type(oversetwall), target :: fullWall
        integer(kind=intType), dimension(:), allocatable :: link, indicesToGet
 
        ! Data for the ADT
        integer(kind=intType) :: intInfo(3), intInfo2(3)
        real(kind=realtype) :: coor(4), uvw(5), uvw2(5)
        real(kind=realtype), dimension(3, 2) :: dummy
        real(kind=realtype), parameter :: tol = 1e-12
        integer(kind=intType), dimension(:), pointer :: frontLeaves, frontLeavesNew, &
                                                        BBint, wallFamList
        type(adtbboxtargettype), dimension(:), pointer :: BB
        real(kind=realtype), dimension(3) :: xp
 
        ! The first thing we do is gather all the surface nodes to
        ! each processor such that every processor can make it's own copy of
        ! the complete surface mesh to use to search. Note that this
        ! procedure *DOES NOT SCALE IN MEMORY*...ie eventually the surface
        ! mesh will become too large to store on a single processor,
        ! although this will probably not happen until the sizes get up in
        ! the hundreds of millions of cells.
 
        allocate (walls(nclusters))
        wallfamlist => bcfamgroups(ibcgroupwalls)%famList
        call buildclusterwalls(level, sps, .false., walls, wallfamlist, size(wallfamlist))
 
        if (oversetpresent) then
            ! Finally build up a "full wall" that is made up of all the cluster
            ! walls.
 
            nnodes = 0
            ncells = 0
            do i = 1, nclusters
                nnodes = nnodes + walls(i)%nNodes
                ncells = ncells + walls(i)%nCells
            end do
 
            allocate (fullwall%x(3, nnodes))
            allocate (fullwall%conn(4, ncells))
            allocate (fullwall%ind(nnodes))
 
            nnodes = 0
            ncells = 0
            ii = 0
            do i = 1, nclusters
 
                ! Add in the nodes/elements from this cluster
 
                do j = 1, walls(i)%nNodes
                    nnodes = nnodes + 1
                    fullwall%x(:, nnodes) = walls(i)%x(:, j)
                    fullwall%ind(nnodes) = walls(i)%ind(j)
                end do
 
                do j = 1, walls(i)%nCells
                    ncells = ncells + 1
                    fullwall%conn(:, ncells) = walls(i)%conn(:, j) + ii
                end do
 
                ! Increment the node offset
                ii = ii + walls(i)%nNodes
            end do
 
            ! Finish the setup of the full wall.
            fullwall%nCells = ncells
            fullwall%nNodes = nnodes
            call buildserialquad(ncells, nnodes, fullwall%x, fullwall%conn, fullwall%ADT)
        end if
 
        ! Allocate the (pointer) memory that may be resized as necessary for
        ! the singlePoint search routine.
        allocate (stack(100), bb(20), bbint(20), frontleaves(25), frontleavesnew(25))
 
        ! We need to store the 4 global node indices defining the quad that
        ! each point has the closest point wrt. We also ned to store the uv
        ! values. This allows us to recompute the exact surface point, after
        ! the rquired nodes are fetched from (a possibly) remote proc.
 
        do nn = 1, ndom
            call setpointers(nn, level, sps)
 
            ! Check if elemID and uv are allocated yet.
            if (.not. associated(flowdoms(nn, level, sps)%surfNodeIndices)) then
                allocate (flowdoms(nn, level, sps)%surfNodeIndices(4, 2:il, 2:jl, 2:kl))
                allocate (flowdoms(nn, level, sps)%uv(2, 2:il, 2:jl, 2:kl))
            end if
 
            ! Set the cluster for this block
            c = clusters(cumdomproc(myid) + nn)
 
            do k = 2, kl
                do j = 2, jl
                    do i = 2, il
 
                        ! Compute the coordinates of the cell center
                        coor(1) = eighth * (x(i - 1, j - 1, k - 1, 1) + x(i, j - 1, k - 1, 1) &
                                            + x(i - 1, j, k - 1, 1) + x(i, j, k - 1, 1) &
                                            + x(i - 1, j - 1, k, 1) + x(i, j - 1, k, 1) &
                                            + x(i - 1, j, k, 1) + x(i, j, k, 1))
 
                        coor(2) = eighth * (x(i - 1, j - 1, k - 1, 2) + x(i, j - 1, k - 1, 2) &
                                            + x(i - 1, j, k - 1, 2) + x(i, j, k - 1, 2) &
                                            + x(i - 1, j - 1, k, 2) + x(i, j - 1, k, 2) &
                                            + x(i - 1, j, k, 2) + x(i, j, k, 2))
 
                        coor(3) = eighth * (x(i - 1, j - 1, k - 1, 3) + x(i, j - 1, k - 1, 3) &
                                            + x(i - 1, j, k - 1, 3) + x(i, j, k - 1, 3) &
                                            + x(i - 1, j - 1, k, 3) + x(i, j - 1, k, 3) &
                                            + x(i - 1, j, k, 3) + x(i, j, k, 3))
 
                        if (.not. oversetpresent) then
                            ! No overset present. Simply search our own wall,
                            ! walls(c), (the only one we have) up to the wall
                            ! cutoff.
                            coor(4) = walldistcutoff**2
                            intinfo(3) = 0 ! Must be initialized since the search
                            ! may not find closer point.
                            call mindistancetreesearchsinglepoint(walls(c)%ADT, coor, intinfo, &
                                                                  uvw, dummy, 0, bb, frontleaves, frontleavesnew)
 
                            cellid = intinfo(3)
                            if (cellid > 0) then
                                do kk = 1, 4
                                    flowdoms(nn, level, sps)%surfNodeIndices(kk, i, j, k) = &
                                        walls(c)%ind(walls(c)%conn(kk, cellid))
                                end do
                                flowdoms(nn, level, sps)%uv(:, i, j, k) = uvw(1:2)
                            else
                                ! Just set dummy values. These will never be used.
                                flowdoms(nn, level, sps)%surfNodeIndices(:, i, j, k) = 0
                                flowdoms(nn, level, sps)%uv(:, i, j, k) = 0
                            end if
 
                            ! We are done with this point.
                            cycle
                        end if
 
                        ! This is now the overset (possibly) overlapping surface
                        ! mesh case. It is somewhat more complex since we use
                        ! the same searches to flag cells that are inside the
                        ! body.
 
                        coor(4) = walldistcutoff**2
                        intinfo(3) = 0
                        call mindistancetreesearchsinglepoint(fullwall%ADT, coor, &
                                                              intinfo, uvw, dummy, 0, bb, frontleaves, frontleavesnew)
                        cellid = intinfo(3)
 
                        if (cellid > 0) then
                            ! We found the cell:
 
                            ! If the cell is outside of near-wall distance or our
                            ! cluster doesn't have any owned cells. Just accept it.
                            if (uvw(4) > nearwalldist**2 .or. walls(c)%nCells == 0) then
 
                                do kk = 1, 4
                                    flowdoms(nn, level, sps)%surfNodeIndices(kk, i, j, k) = &
                                        fullwall%ind(fullwall%conn(kk, cellid))
                                end do
                                flowdoms(nn, level, sps)%uv(:, i, j, k) = uvw(1:2)
 
                            else
 
                                ! This point is *closer* than the nearWallDist AND
                                ! it has a wall. Search on our own wall.
 
                                coor(4) = large
                                call mindistancetreesearchsinglepoint(walls(c)%ADT, coor, &
                                                                      intinfo2, uvw2, dummy, 0, bb, &
                                                                      frontleaves, frontleavesnew)
                                cellid2 = intinfo2(3)
 
                                if (uvw2(4) < nearwalldist**2) then
                                    ! Both are close to the wall. Accept the one
                                    ! from our own wall unconditionally.
                                    do kk = 1, 4
                                        flowdoms(nn, level, sps)%surfNodeIndices(kk, i, j, k) = &
                                            walls(c)%ind(walls(c)%conn(kk, cellid2))
                                    end do
                                    flowdoms(nn, level, sps)%uv(:, i, j, k) = uvw2(1:2)
                                else
                                    ! The full wall distance is better. Take that.
 
                                    do kk = 1, 4
                                        flowdoms(nn, level, sps)%surfNodeIndices(kk, i, j, k) = &
                                            fullwall%ind(fullwall%conn(kk, cellid))
                                    end do
                                    flowdoms(nn, level, sps)%uv(:, i, j, k) = uvw(1:2)
 
                                end if
                            end if
                        else
 
                            ! What happend here is a cell is outside the
                            ! wallDistCutoff. We don't care about wall distance
                            ! info here so just set dummy info.
 
                            flowdoms(nn, level, sps)%surfNodeIndices(:, i, j, k) = 0
                            flowdoms(nn, level, sps)%uv(:, i, j, k) = 0
 
                        end if
                    end do
                end do
            end do
        end do
 
        ! Now determine all the node indices this processor needs to get.
        mm = 0
        allocate (indicestoget(ncellslocal(level) * 4), link(ncellslocal(level) * 4))
        do nn = 1, ndom
            call setpointers(nn, level, sps)
            do k = 2, kl
                do j = 2, jl
                    do i = 2, il
                        do kk = 1, 4
                            mm = mm + 1
                            indicestoget(mm) = flowdoms(nn, level, sps)%surfNodeIndices(kk, i, j, k)
                        end do
                    end do
                end do
            end do
        end do
 
        ! This unique-ifies the indices.
        call unique(indicestoget, 4 * ncellslocal(level), nunique, link)
 
        ! we need to update the stored indices to use the ordering of the nodes we will receive.
        mm = 0
        do nn = 1, ndom
            call setpointers(nn, level, sps)
            do k = 2, kl
                do j = 2, jl
                    do i = 2, il
                        do kk = 1, 4
                            mm = mm + 1
                            flowdoms(nn, level, sps)%surfNodeIndices(kk, i, j, k) = link(mm)
                        end do
                    end do
                end do
            end do
        end do
        deallocate (link)
 
        ! Now create the index set for the nodes we need to get. We have to
        ! expand "indices to get" to include the DOF. Use link for this
        ! temporary array operation.
 
        allocate (link(nunique * 3))
        do i = 1, nunique
            link((i - 1) * 3 + 1) = indicestoget(i) * 3
            link((i - 1) * 3 + 2) = indicestoget(i) * 3 + 1
            link((i - 1) * 3 + 3) = indicestoget(i) * 3 + 2
        end do
 
        call iscreategeneral(adflow_comm_world, nunique * 3, link, petsc_copy_values, is1, ierr)
        call echk(ierr, __file__, __line__)
        deallocate (link)
 
        ! Create the volume vector the nodes will be scatter from. Note that
        ! this vector contains all the spectal instances. It is therefore
        ! only allocated on the first call with sps=1
        if (sps == 1) then
            call veccreatempi(adflow_comm_world, 3 * nnodeslocal(level) * ntimeintervalsspectral, &
                              petsc_determine, xvolumevec(level), ierr)
            call echk(ierr, __file__, __line__)
        end if
 
        ! This is the vector we will scatter the nodes into.
        call veccreatempi(adflow_comm_world, 3 * nunique, petsc_determine, &
                          xsurfvec(level, sps), ierr)
        call echk(ierr, __file__, __line__)
 
        call vecgetownershiprange(xsurfvec(level, sps), i, j, ierr)
        call echk(ierr, __file__, __line__)
 
        call iscreatestride(adflow_comm_world, j - i, i, 1, is2, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Create the actual final scatter context.
        call vecscattercreate(xvolumevec(level), is1, xsurfvec(level, sps), is2, &
                              wallscatter(level, sps), ierr)
        call echk(ierr, __file__, __line__)
 
        call isdestroy(is1, ierr)
        call echk(ierr, __file__, __line__)
 
        call isdestroy(is2, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Deallocate all the remaining temporary data
        deallocate (stack, bb, frontleaves, frontleavesnew, bbint)
 
        do i = 1, nclusters
            deallocate (walls(i)%x, walls(i)%conn, walls(i)%ind)
            call destroyserialquad(walls(i)%ADT)
        end do
        deallocate (walls)
 
        if (oversetpresent) then
            deallocate (fullwall%x, fullwall%conn, fullwall%ind)
            call destroyserialquad(fullwall%ADT)
        end if
 
    end subroutine determinewallassociation
 
    subroutine updatexsurf(level)
 
        use blockpointers
        use inputtimespectral
        use utils, only: echk, setpointers
        implicit none
 
        ! Input Parameters
        integer(kind=intType), intent(in) :: level
 
        ! Working Parameters
        integer(kind=intType) :: ii, i, j, k, l, nn, sps, ierr
 
        ! Fill up xVolumeVec
        call vecgetarrayf90(xvolumevec(level), xvolume, ierr)
        call echk(ierr, __file__, __line__)
 
        ii = 0
        do nn = 1, ndom
            do sps = 1, ntimeintervalsspectral
                call setpointers(nn, level, sps)
                do k = 1, kl
                    do j = 1, jl
                        do i = 1, il
                            do l = 1, 3
                                ii = ii + 1
                                xvolume(ii) = x(i, j, k, l)
                            end do
                        end do
                    end do
                end do
            end do
        end do
        call vecrestorearrayf90(xvolumevec(level), xvolume, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Perform the scatter from the global x vector to xSurf. SPS loop since the xSurfVec is done by SPS instance.
        do sps = 1, ntimeintervalsspectral
            call vecscatterbegin(wallscatter(level, sps), xvolumevec(level), &
                                 xsurfvec(level, sps), insert_values, scatter_forward, ierr)
            call echk(ierr, __file__, __line__)
 
            call vecscatterend(wallscatter(level, sps), xvolumevec(level), &
                               xsurfvec(level, sps), insert_values, scatter_forward, ierr)
            call echk(ierr, __file__, __line__)
        end do
 
    end subroutine updatexsurf
 
    subroutine destroywalldistancedata
 
        use constants
        use block, only: flowdoms
        implicit none
 
        ! Working
        integer(kind=intType) :: level, nLevels, l
 
        nlevels = ubound(flowdoms, 2)
        do l = 1, nlevels
            call destroywalldistancedatalevel(l)
        end do
 
        deallocate (xsurfvec, xvolumevec, wallscatter)
 
    end subroutine destroywalldistancedata
 
    subroutine destroywalldistancedatalevel(level)
        use constants
        use inputtimespectral, only: ntimeintervalsspectral
        use utils, onlY: echk
        use block, only: flowdoms
 
        implicit none
 
        ! Input Parameters
        integer(kind=intType), intent(in) :: level
 
        ! Working
        integer(kind=intType) :: ierr, sps
 
        ! Determine if we need to deallocate the PETSc data for
        ! this level
        if (walldistancedataallocated(level)) then
            call vecdestroy(xvolumevec(level), ierr)
            call echk(ierr, __file__, __line__)
 
            do sps = 1, ntimeintervalsspectral
                call vecdestroy(xsurfvec(level, sps), ierr)
                call echk(ierr, __file__, __line__)
 
                call vecscatterdestroy(wallscatter(level, sps), ierr)
                call echk(ierr, __file__, __line__)
            end do
 
            walldistancedataallocated(level) = .false.
        end if
    end subroutine destroywalldistancedatalevel
 
#endif
end module walldistance