tecplotIO.F90

Go to the documentation of this file.

module tecplotio
 
    use constants, only: realtype, inttype, maxstringlen, maxcgnsnamelen
    implicit none
    save
 
    character(len=maxStringLen) :: sci6 = "(ES14.6)"
 
    type slice
 
        ! nNodes : Number of nodes for this slice
        ! ind(2, nNodes) : Indices of the nodes in the global node list on either
        !                  side of node
        ! w(2, nNodes) : Weights used to multiply the two global nodes defined in
        !                ind to get compute nodal values (positions, forces etc)
        ! pL, vL, pD, vD, pM, vM : Pressure and viscous components of lift, drag, and moment
        ! CLp, CLv, CDp, CDv, CMp, CMv : Coefficients of pressure and viscous lift, drag, and moment
        ! chord: chord of section
        ! pt, normal: The point and the normal that defines the slicing plane
        ! dir_vec: a direction vector that we use to filter sliced line elements.
        ! use_dir: flag to determine if we use the dir vec or not. if we are doing
        !          a regular slice, e.g. wing section, we dont want to use the dir
        !          and want to use the full slice. if we are doing a cylindrical slice
        !          we then want to pick which direction since we would get 2 slices
        !          on something like a nacelle in this case.
 
        character(len=maxStringLen) :: slicename
        integer(kind=intType) :: sps
        integer(kind=intType), dimension(:, :), allocatable :: ind, conn
        real(kind=realtype), dimension(:, :), allocatable :: w, vars
        integer(kind=intType) :: nnodes
        real(kind=realtype) :: pl, vl, pd, vd, pm, vm, clp, clv, cdp, cdv, cmp, cmv
        real(kind=realtype) :: chord, twist, thickness
        real(kind=realtype), dimension(3) :: le, te
        real(kind=realtype), dimension(3) :: pt, normal, dir_vec
        logical :: use_dir
        integer(kind=intType), allocatable, dimension(:) :: famlist
        ! here we declare that 'slice'-types also now have fx,fy and fz:
        real(kind=realtype) :: fx, fy, fz
    end type slice
 
    type liftdist
        ! nSegments: Number of nodes to use for distribution
        ! normal: Slice direction (normal of the plane)
        ! normal_ind: Index of direction..1 for x, 2 for y, 3 for z
        ! distName: Name of lift distribution
        ! slices: The list of slices this distribution will use
        ! delta: The current delta spacing for the distribution
        ! slicePoints: The list of points where the slices are taken
        character(len=maxStringLen) :: distname
        integer(kind=intType) :: nsegments, normal_ind
        integer(kind=intType), dimension(:), allocatable :: famlist
        real(kind=realtype) :: normal(3)
        real(kind=realtype) :: delta
        real(kind=realtype), dimension(:, :), allocatable :: slicepts
    end type liftdist
 
    logical :: liftdistinitialized = .false.
    integer(kind=intType) :: mscon1(16, 5), mscon2(4, 2)
 
    ! Data for the user supplied slices:
    integer(kind=intType), parameter :: nslicemax = 1000
    integer(kind=intType) :: nparaslices = 0
    integer(kind=intType) :: nabsslices = 0
    type(slice), dimension(:, :), allocatable :: paraslices, absslices
 
    ! Data for the user supplied lift distributions
    integer(kind=intType), parameter :: nliftdistmax = 100
    integer(kind=intType) :: nliftdists = 0
    type(liftdist), dimension(nLiftDistMax), target :: liftdists
 
    ! Tecplot Variable names of the data in the lift distribution data file:
    character(len=maxCGNSNameLen), dimension(:), allocatable :: liftdistname
    integer(kind=intType), parameter :: nliftdistvar = 26
 
contains
    subroutine addparaslice(sliceName, pt, normal, dir_vec, use_dir, famList, n)
        !
        !       This subroutine is intended to be called from python.
        !       This routine will add a parametric slice to the list of user
        !       supplied slices.
        use constants
        use communication
        use surfacefamilies
        use surfaceutils
        use inputtimespectral
        implicit none
 
        ! Input parameters
        character(len=*), intent(in) :: sliceName
        real(kind=realtype), dimension(3), intent(in) :: pt, normal, dir_vec
        logical, intent(in) :: use_dir
        integer(kind=intType), intent(in) :: n, famList(n)
 
        ! Working
        integer(kind=intType) :: sps, sizeNode, sizeCell
        integer(kind=intType), dimension(:), pointer :: wallList
        real(kind=realtype), dimension(:, :), allocatable :: pts
        integer(kind=intType), dimension(:, :), allocatable :: conn
        integer(kind=intType), dimension(:), allocatable :: elemFam, cgnsBlockID
 
        if (.not. allocated(paraslices)) then
            allocate (paraslices(nslicemax, ntimeintervalsspectral))
        end if
 
        ! We have to add a slice for each spectral instance.
        do sps = 1, ntimeintervalsspectral
            nparaslices = nparaslices + 1
 
            if (nparaslices > nslicemax) then
                print *, 'Error: Exceeded the maximum number of slices. Increase nSliceMax'
                stop
            end if
 
            ! Slices are created on walls and walls only. Retrieve the
            ! points, connectivity and familyID of all the walls.
            walllist => bcfamgroups(ibcgroupwalls)%famList
            call getsurfacesize(sizenode, sizecell, walllist, size(walllist), .true.)
            allocate (pts(3, sizenode), conn(4, sizecell), elemfam(sizecell), cgnsblockid(sizecell))
            call getsurfaceconnectivity(conn, cgnsblockid, sizecell, walllist, size(walllist), .true.)
            call getsurfacepoints(pts, sizenode, sps, walllist, size(walllist), .true.)
            call getsurfacefamily(elemfam, sizecell, walllist, size(walllist), .true.)
 
            ! Create actual slice
            call createslice(pts, conn, elemfam, paraslices(nparaslices, sps), pt, normal, dir_vec, &
                             use_dir, slicename, famlist)
 
            ! Clean up memory.
            deallocate (pts, conn, elemfam)
        end do
 
    end subroutine addparaslice
 
    subroutine addabsslice(sliceName, pt, normal, dir_vec, use_dir, famList, n)
        !
        !       This subroutine is intended to be called from python.
        !       This routine will add an absolute slice to the list of user
        !       supplied slices.
        use constants
        use communication
        use surfacefamilies
        use surfaceutils
        use inputtimespectral
        implicit none
 
        ! Input parameters
        character(len=*), intent(in) :: sliceName
        real(kind=realtype), dimension(3), intent(in) :: pt, normal, dir_vec
        logical, intent(in) :: use_dir
        integer(kind=intType), intent(in) :: n, famList(n)
 
        ! Working
        integer(kind=intType) :: sps, sizeNode, sizeCell
        integer(kind=intType), dimension(:), pointer :: wallList
        real(kind=realtype), dimension(:, :), allocatable :: pts
        integer(kind=intType), dimension(:, :), allocatable :: conn
        integer(kind=intType), dimension(:), allocatable :: elemFam, cgnsBlockID
 
        if (.not. allocated(absslices)) then
            allocate (absslices(nslicemax, ntimeintervalsspectral))
        end if
 
        do sps = 1, ntimeintervalsspectral
            nabsslices = nabsslices + 1
 
            if (nabsslices > nslicemax) then
                print *, 'Error: Exceeded the maximum number of slices. Increase nSliceMax'
                stop
            end if
 
            walllist => bcfamgroups(ibcgroupwalls)%famList
            call getsurfacesize(sizenode, sizecell, walllist, size(walllist), .true.)
            allocate (pts(3, sizenode), conn(4, sizecell), elemfam(sizecell), cgnsblockid(sizecell))
            call getsurfaceconnectivity(conn, cgnsblockid, sizecell, walllist, size(walllist), .true.)
            call getsurfacepoints(pts, sizenode, sps, walllist, size(walllist), .true.)
            call getsurfacefamily(elemfam, sizecell, walllist, size(walllist), .true.)
            call createslice(pts, conn, elemfam, absslices(nabsslices, sps), pt, normal, dir_vec, &
                             use_dir, slicename, famlist)
 
            ! Clean up memory.
            deallocate (pts, conn, elemfam)
        end do
 
    end subroutine addabsslice
 
    subroutine addliftdistribution(nSegments, normal, normal_ind, distName, famList, n)
        !
        !       This subroutine is intended to be called from python.
        !       This routine will add the description of a lift distribution
 
        use constants
        use communication
        use surfacefamilies
        implicit none
 
        ! Input parameters
        character(len=*), intent(in) :: distName
        integer(kind=intType), intent(in) :: nSegments
        real(kind=realtype), dimension(3) :: normal
        integer(kind=intType), intent(in) :: normal_ind
        integer(kind=intType), intent(in) :: n, famList(n)
 
        nliftdists = nliftdists + 1
        if (nliftdists > nliftdistmax) then
            print *, 'Error: Exceeded the maximum number of lift distributions. &
                 &Increase nLiftDistMax'
            stop
        end if
 
        liftdists(nliftdists)%nSegments = nsegments
        liftdists(nliftdists)%normal = normal
        liftdists(nliftdists)%normal_ind = normal_ind
        liftdists(nliftdists)%distName = distname
 
        allocate (liftdists(nliftdists)%famList(n))
        liftdists(nliftdists)%famList(:) = famlist
 
    end subroutine addliftdistribution
 
    subroutine writetecplot(sliceFile, writeSlices, liftFile, writeLift, &
                            surfFile, writeSurf, famList, nFamList)
        !
        !       This is the master routine for writing tecplot data from adflow.
        !       This routine will write the slice, lift and surface files
        !       depending on the flags writeSlics, writeLift and writeSurface.
        !       The reason for the combined routine is that we can safely only
        !       perform the nodal averaging once which is required for all
        !       three output files.
        use constants
        use inputtimespectral, only: ntimeintervalsspectral
        use surfacefamilies, only: bcfamgroups, bcfamexchange
        use surfaceutils, only: getsurfacesize
        use oversetdata, only: zippermeshes
        use outputmod, only: numberofsurfsolvariables
        implicit none
 
        ! Input Params
        character(len=*), intent(in) :: sliceFile, liftFile, surfFile
        logical, intent(in) :: writeSlices, writeLift, writeSurf
        integer(kind=intType), intent(in) :: nFamList
        integer(kind=intType), intent(in), dimension(nFamList) :: famList
        real(kind=realtype), dimension(:, :, :), allocatable :: nodalvalues
        ! Working
        integer(kind=intType) :: sps, nSolVar, sizeNOde, sizeCell
        integer(kind=intType), dimension(:), pointer :: wallLIST
 
        ! Determine the number of surface variables we have
        call numberofsurfsolvariables(nsolvar)
        walllist => bcfamgroups(ibcgroupwalls)%famList
        call getsurfacesize(sizenode, sizecell, walllist, size(walllist), .true.)
 
        ! Allocate and compute the wall-based surface data for hte slices
        ! and lift distributions.
        allocate (nodalvalues(max(sizenode, 1), nsolvar + 6 + 3, ntimeintervalsspectral))
 
        do sps = 1, ntimeintervalsspectral
            call computesurfaceoutputnodaldata(bcfamexchange(ibcgroupwalls, sps), &
                                               zippermeshes(ibcgroupwalls), .true., nodalvalues(:, :, sps))
        end do
 
        if (writeslices) then
            call writeslicesfile(slicefile, nodalvalues)
        end if
 
        if (writelift) then
            call writeliftdistributionfile(liftfile, nodalvalues)
        end if
 
        deallocate (nodalvalues)
 
        if (writesurf) then
            call writetecplotsurfacefile(surffile, famlist)
        end if
 
    end subroutine writetecplot
 
    subroutine writeslicesfile(fileName, nodalValues)
        !
        ! This subroutine is intended to be called from python.  This
        ! routine will write the user defined slics to an to the (ascii)
        ! tecplot file fileName. ASCII files are used for simplicity since
        ! very little information is actually written.
        use constants
        use communication
        use outputmod
        use inputtimespectral
        use inputphysics
        use inputiteration
        use inputio
        use surfacefamilies
        use surfaceutils
        use utils, only: echk
        implicit none
 
        ! Input Params
        character(len=*), intent(in) :: fileName
        real(kind=realtype), intent(inout), dimension(:, :, :) :: nodalvalues
 
        ! Working parameters
        integer(kind=intType) :: file, i, sps, nSolVar, ierr
        character(len=maxStringLen) :: fname
        character(len=7) :: intString
        character(len=maxCGNSNameLen), dimension(:), allocatable :: solNames
        integer(kind=intType), allocatable, dimension(:) :: famList
        type(slice) :: globalSlice
        integer(kind=intType) :: sizeNode, sizeCell
        integer(kind=intType), dimension(:), pointer :: wallList
        real(kind=realtype), dimension(:, :), allocatable :: pts
        integer(kind=intType), dimension(:, :), allocatable :: conn
        integer(kind=intType), dimension(:), allocatable :: elemFam, cgnsBlockID
 
        ! Only write if we actually have lift distributions
        testwriteslices: if (nparaslices + nabsslices > 0) then
 
            if (myid == 0 .and. printiterations) then
                print "(a)", "#"
                print "(a)", "# Writing slices file(s): "
            end if
 
            do sps = 1, ntimeintervalsspectral
 
                ! If it is time spectral we need to agument the filename
                if (equationmode == timespectral) then
                    write (intstring, "(i7)") sps
                    intstring = adjustl(intstring)
                    fname = trim(filename)//"Spectral"//trim(intstring)
                else
                    fname = filename
                end if
 
                file = 11
                ! Open file on root proc:
                if (myid == 0) then
                    ! Print the filename to stdout
                    print "(a,4x,a)", "#", trim(fname)
                    open (unit=file, file=trim(fname))
 
                    ! Write Header Information
                    write (file, *) "Title = ""ADflow Slice Data"""
                    write (file, "(a)", advance="no") "Variables = "
                    write (file, "(a)", advance="no") " ""CoordinateX"" "
                    write (file, "(a)", advance="no") " ""CoordinateY"" "
                    write (file, "(a)", advance="no") " ""CoordinateZ"" "
                    write (file, "(a)", advance="no") " ""XoC"" "
                    write (file, "(a)", advance="no") " ""YoC"" "
                    write (file, "(a)", advance="no") " ""ZoC"" "
 
                    ! Number of additional variables on slices
                    call numberofsurfsolvariables(nsolvar)
                    allocate (solnames(nsolvar))
                    call surfsolnames(solnames)
 
                    ! Write the rest of the variables
                    do i = 1, nsolvar
                        write (file, "(a,a,a)", advance="no") """", trim(solnames(i)), """ "
                    end do
 
                    write (file, "(1x)")
                    deallocate (solnames)
                end if
                call mpi_bcast(nsolvar, 1, adflow_integer, 0, adflow_comm_world, ierr)
                call echk(ierr, __file__, __line__)
 
                ! Slices are created on walls and walls only. Retrieve the
                ! points, connectivity and familyID of all the walls.
                walllist => bcfamgroups(ibcgroupwalls)%famList
                call getsurfacesize(sizenode, sizecell, walllist, size(walllist), .true.)
                allocate (pts(3, sizenode), conn(4, sizecell), elemfam(sizecell), cgnsblockid(sizecell))
                call getsurfaceconnectivity(conn, cgnsblockid, sizecell, walllist, size(walllist), .true.)
                call getsurfacepoints(pts, sizenode, sps, walllist, size(walllist), .true.)
                call getsurfacefamily(elemfam, sizecell, walllist, size(walllist), .true.)
 
                ! Integration is performed in parallel
                do i = 1, nparaslices
                    call integrateslice(paraslices(i, sps), globalslice, &
                                        nodalvalues(:, :, sps), nsolvar, .true.)
                    if (myid == 0) then
                        call writeslice(globalslice, file, nsolvar)
                    end if
                    call destroyslice(globalslice)
                end do
 
                do i = 1, nabsslices
                    ! 'Destroy' the slice...just dealloc the allocated data.
                    ! before we do, save the family list
                    allocate (famlist(size(absslices(i, sps)%famList)))
                    famlist = absslices(i, sps)%famList
                    call destroyslice(absslices(i, sps))
 
                    ! Make new one in the same location
                    call createslice(pts, conn, elemfam, absslices(i, sps), &
                                     absslices(i, sps)%pt, absslices(i, sps)%normal, absslices(i, sps)%dir_vec, &
                                     absslices(i, sps)%use_dir, absslices(i, sps)%sliceName, famlist)
 
                    call integrateslice(absslices(i, sps), globalslice, &
                                        nodalvalues(:, :, sps), nsolvar, .true.)
                    if (myid == 0) then
                        call writeslice(globalslice, file, nsolvar)
                    end if
                    call destroyslice(globalslice)
                    deallocate (famlist)
                end do
 
                !Close file on root proc
                if (myid == 0) then
                    close (file)
                end if
            end do
 
            if (myid == 0 .and. printiterations) then
                print "(a)", "# Slices file(s) written"
                print "(a)", "#"
            end if
        end if testwriteslices
    end subroutine writeslicesfile
 
    subroutine writeliftdistributionfile(fileName, nodalValues)
        !
        !
        !       This subroutine is intended to be called from python.
        !       This routine will write the added lift distributions
        !       to the (ascii) tecplot file fileName. ASCII files are
        !       used for siplicity since very little informatin is actually
        !       written.
        use constants
        use communication
        use outputmod
        use inputphysics
        use inputtimespectral
        use inputiteration
        use surfacefamilies
        implicit none
 
        ! Input Params
        character(len=*), intent(in) :: fileName
        real(kind=realtype), dimension(:, :, :), allocatable :: nodalvalues
 
        ! Working parameters
        integer(kind=intType) :: file, sps
        character(len=maxStringLen) :: fname
        character(len=7) :: intString
 
        ! Only write if we actually have lift distributions
        testwriteliftdists: if (nliftdists > 0) then
 
            if (myid == 0 .and. printiterations) then
                print "(a)", "#"
                print "(a)", "# Writing lift distribution file(s):"
            end if
 
            do sps = 1, ntimeintervalsspectral
 
                ! If it is time spectral we need to agument the filename
                if (equationmode == timespectral) then
                    write (intstring, "(i7)") sps
                    intstring = adjustl(intstring)
                    fname = trim(filename)//"Spectral"//trim(intstring)
                else
                    fname = filename
                end if
 
                file = 11
                ! Open file on root proc:
                if (myid == 0) then
                    ! Print the filename to stdout
                    print "(a,4x,a)", "#", trim(fname)
 
                    open (unit=file, file=trim(fname))
                end if
 
                call writeliftdistributions(sps, file, nodalvalues(:, :, sps))
 
                ! Close file on root proc
                if (myid == 0) then
                    close (file)
                end if
            end do
 
            if (myid == 0 .and. printiterations) then
                print "(a)", "# Lift distribution file(s) written"
                print "(a)", "#"
            end if
 
        end if testwriteliftdists
    end subroutine writeliftdistributionfile
 
    subroutine writeliftdistributions(sps, fileID, nodalValues)
        !
        !       This subroutine writes the liftdistribution for the specified
        !       spectral instance. It is assumed that the required file handles
        !       are already open and can be written to
        use constants
        use communication
        use outputmod
        use su_cgns
        use cgnsnames
        use surfacefamilies, only: bcfamgroups, bcfamexchange
        use surfaceutils
        use utils, only: echk
        use sorting, only: faminlist
        implicit none
 
        ! Input parameters
        integer(kind=intType), intent(in) :: sps, fileID
        real(kind=realtype), dimension(:, :), intent(in) :: nodalvalues
 
        real(kind=realtype), dimension(3) :: xmin, xmax, xmin_local, xmax_local
        real(kind=realtype), parameter :: tol = 1e-8
        type(liftdist), pointer :: d
        integer(kind=intType) :: i, j, ii, jj, idist, ierr
        real(kind=realtype), dimension(:, :), allocatable :: values
        character(len=maxCGNSNameLen), dimension(:), allocatable :: liftdistnames
        real(kind=realtype) :: dmin, dmax, suml, sumd, summ, span, delta, xcur(3)
        type(slice) :: localslice, globalslice
        integer(kind=intType) :: sizenode, sizecell
        integer(kind=intType), dimension(:), pointer :: walllist
        real(kind=realtype), dimension(:, :), allocatable :: pts
        integer(kind=intType), dimension(:, :), allocatable :: conn
        integer(kind=intType), dimension(:), allocatable :: elemfam, cgnsblockid
 
        ! Slices are created on walls and walls only. Retrieve the
        ! points, connectivity and familyID of all the walls.
        walllist => bcfamgroups(ibcgroupwalls)%famList
        call getsurfacesize(sizenode, sizecell, walllist, size(walllist), .true.)
        allocate (pts(3, sizenode), conn(4, sizecell), elemfam(sizecell), cgnsblockid(sizecell))
        call getsurfaceconnectivity(conn, cgnsblockid, sizecell, walllist, size(walllist), .true.)
        call getsurfacepoints(pts, sizenode, sps, walllist, size(walllist), .true.)
        call getsurfacefamily(elemfam, sizecell, walllist, size(walllist), .true.)
 
        do idist = 1, nliftdists
 
            d => liftdists(idist)
            xmin_local = huge(real(zero))
            xmax_local = -huge(real(zero))
 
            ! Get the bounding box for the entire geometry we have been slicing.
            elemloop: do i = 1, size(conn, 2)
                if (faminlist(elemfam(i), d%FamList)) then
 
                    ! Extract each of the 4 nodes on this quad:
                    do jj = 1, 4
                        xcur = pts(:, conn(jj, i))
                        ! Check the max/min on each index
                        do ii = 1, 3
                            xmin_local(ii) = min(xmin_local(ii), xcur(ii))
                            xmax_local(ii) = max(xmax_local(ii), xcur(ii))
                        end do
                    end do
                end if
            end do elemloop
 
            ! Globalize all min/max values.
            call mpi_allreduce(xmin_local, xmin, 3, adflow_real, mpi_min, &
                               adflow_comm_world, ierr)
            call echk(ierr, __file__, __line__)
 
            call mpi_allreduce(xmax_local, xmax, 3, adflow_real, mpi_max, &
                               adflow_comm_world, ierr)
            call echk(ierr, __file__, __line__)
 
            d%delta = (xmax(d%normal_ind) - xmin(d%normal_ind)) / dble((d%nSegments - 1))
            allocate (d%slicePts(3, d%nSegments))
 
            ! Zero out all segments
            d%slicePts = zero
 
            ! These are the variable names for the lift distribution:
            allocate (liftdistnames(nliftdistvar))
            ! Set the names here
            liftdistnames(1) = "eta"
            liftdistnames(2) = cgnscoorx
            liftdistnames(3) = cgnscoory
            liftdistnames(4) = cgnscoorz
            liftdistnames(5) = "Lift"
            liftdistnames(6) = "Drag"
            liftdistnames(7) = "Moment"
            liftdistnames(8) = "Normalized Lift"
            liftdistnames(9) = "Normalized Drag"
            liftdistnames(10) = "Normalized Moment"
            liftdistnames(11) = "CL"
            liftdistnames(12) = "CD"
            liftdistnames(13) = "CM"
            liftdistnames(14) = "CLp"
            liftdistnames(15) = "CDp"
            liftdistnames(16) = "CMp"
            liftdistnames(17) = "CLv"
            liftdistnames(18) = "CDv"
            liftdistnames(19) = "CMv"
            liftdistnames(20) = "Elliptical"
            liftdistnames(21) = "thickness"
            liftdistnames(22) = "twist"
            liftdistnames(23) = "chord"
            liftdistnames(24) = "fx"
            liftdistnames(25) = "fy"
            liftdistnames(26) = "fz"
 
            ! Only write header info for first distribution only
            if (myid == 0) then
                if (idist == 1) then
                    write (fileid, *) "Title= ""ADflow Lift Distribution Data"""
                    write (fileid, "(a)", advance="no") "Variables = "
                    do i = 1, nliftdistvar ! here we just write var-names in the header
                        write (fileid, "(a,a,a)", advance="no") """", trim(liftdistnames(i)), """ "
                    end do
                    write (fileid, "(1x)")
                end if
 
                write (fileid, "(a,a,a)") "Zone T= """, trim(d%distName), """"
                write (fileid, *) "I= ", d%nSegments
                write (fileid, *) "DATAPACKING=BLOCK"
            end if
 
            allocate (values(d%nSegments, nliftdistvar))
            values = zero
 
            do i = 1, d%nSegments
                if (i == 1) then
                    d%slicePts(d%normal_ind, i) = xmin(d%normal_ind) + tol
                else if (i == d%nSegments) then
                    d%slicePts(d%normal_ind, i) = xmin(d%normal_ind) + (i - 1) * d%delta - tol
                else
                    d%slicePts(d%normal_ind, i) = xmin(d%normal_ind) + (i - 1) * d%delta
                end if
            end do
 
            ! Scaled Eta values
            dmin = minval(d%slicePts(d%normal_ind, :))
            dmax = maxval(d%slicePts(d%normal_ind, :))
 
            ! next line we save ETA as the first column. This corresponds to
            ! when we open a lift file, e.g. fc_000_lift.dat where we find
            ! ETA as the first entry... (also visible through tec360)
            values(:, 1) = (d%slicePts(d%normal_ind, :) - dmin) / (dmax - dmin)
            ! Coordinate Varaibles
            if (d%normal_ind == 1) then! X slices
                values(:, 2) = d%slicePts(1, :)
            else if (d%normal_ind == 2) then ! Y slices
                values(:, 3) = d%slicePts(2, :)
            else if (d%normal_ind == 3) then ! Z slices
                values(:, 4) = d%slicePts(3, :)
            end if
 
            !      as mentioned above nSegments are number of nodes in the
            !      distribution, e.g. maybe you have asked for 150 points out along
            !      your wing where you want the lift. Ok, we now for each of these
            !      150 points have to create a slice and integrate our metrics around
            !      said slice, save the integrated metrics (e.g. lift, drag) for the
            !      current point. We then move on to the next point, make a new slice
            !      get the new point's metrics and save them ... we do this 150 times
            do i = 1, d%nSegments
                ! Make new one in the same location
                ! we just pass the normal again as the direction vector because its not used
                call createslice(pts, conn, elemfam, localslice, d%slicePts(:, i), &
                                 d%normal, d%normal, .false., "does_not_matter", d%famList)
                call integrateslice(localslice, globalslice, nodalvalues, 0, .false.)
 
                ! Total lift, drag, and moment
                values(i, 5) = globalslice%pL + globalslice%vL
                values(i, 6) = globalslice%pD + globalslice%vD
                values(i, 7) = globalslice%pM + globalslice%vM
 
                ! Total CL, CD, and CM
                values(i, 11) = globalslice%CLp + globalslice%CLv
                values(i, 12) = globalslice%CDp + globalslice%CDv
                values(i, 13) = globalslice%CMp + globalslice%CMv
 
                ! Pressure lift, drag, and moment coefficients
                values(i, 14) = globalslice%CLp
                values(i, 15) = globalslice%CDp
                values(i, 16) = globalslice%CMp
 
                ! Viscous lift, drag, and moment coefficients
                values(i, 17) = globalslice%CLv
                values(i, 18) = globalslice%CDv
                values(i, 19) = globalslice%CMv
 
                ! t/c, twist, chord
                values(i, 21) = globalslice%thickness
                values(i, 22) = globalslice%twist
                values(i, 23) = globalslice%chord
 
                ! here we now save our new, added values that we have desired to use
                values(i, 24) = globalslice%fx
                values(i, 25) = globalslice%fy
                values(i, 26) = globalslice%fz
 
                call destroyslice(localslice)
                call destroyslice(globalslice)
 
            end do
 
            ! Sum up the lift, drag, and moment values from the slices:
            suml = zero
            sumd = zero
            summ = zero
            do i = 1, d%nSegments - 1
                suml = suml + half * (values(i, 5) + values(i + 1, 5))
                sumd = sumd + half * (values(i, 6) + values(i + 1, 6))
                summ = summ + half * (values(i, 7) + values(i + 1, 7))
            end do
 
            ! Now compute the normalized lift, drag and elliptical since
            ! we know the sum. Note delta is non-dimensional!
            delta = values(2, 1) - values(1, 1)
 
            ! This is the "nondimensional" span...it basically takes into account if you have
            ! a wing not at the sym plane
            span = maxval(values(:, 1)) - minval(values(:, 1))
            dmin = minval(values(:, 1))
            suml = suml * delta
            sumd = sumd * delta
            summ = summ * delta
 
            do i = 1, d%nSegments
 
                ! Normalized Lift, Drag, and Moment
                values(i, 8) = values(i, 5) / suml
                values(i, 9) = values(i, 6) / sumd
                values(i, 10) = values(i, 7) / abs(summ)
 
                ! elliptical value
                values(i, 20) = four / pi / span * sqrt(one - (values(i, 1) - dmin)**2 / span**2)
            end do
 
            ! Write all variables in block format
            if (myid == 0) then
                do j = 1, nliftdistvar
                    do i = 1, d%nSegments
                        write (fileid, sci6) values(i, j)
                    end do
                end do
            end if
 
            ! Deallocate slice list and point list
            deallocate (d%slicePts)
 
            ! Destroy temp variables
            deallocate (liftdistnames, values)
        end do
 
    end subroutine writeliftdistributions
 
    subroutine writetecplotsurfacefile(fileName, famList)
        use constants
        use communication, only: myid, adflow_comm_world, nproc
        use inputtimespectral, only: ntimeintervalsspectral
        use inputphysics, only: equationmode
        use inputiteration, only: printiterations
        use inputio, only: precisionsurfgrid, precisionsurfsol
        use outputmod, only: surfsolnames, numberofsurfsolvariables
        use surfacefamilies, onlY: bcfamexchange, famnames, familyexchange
        use utils, only: echk, setpointers, setbcpointers
        use bcpointers, only: xx
        use sorting, only: faminlist
        use extraoutput, only: surfwriteblank
        use oversetdata, only: zippermesh, zippermeshes
        use surfaceutils
#include <petsc/finclude/petsc.h>
        use petsc
        implicit none
 
        ! Input Params
        character(len=*), intent(in) :: filename
        integer(kind=intType), intent(in), dimension(:) :: famlist
 
        ! Working parameters
        integer(kind=intType) :: i, j, nn, mm, fileid, ivar, ii, ierr, isize
        integer(Kind=intType) :: nsolvar, ibeg, iend, jbeg, jend, sps, sizenode, sizecell
        integer(kind=intType) :: ibcgroup, ifam, iproc, ncells, nnodes, ncellstowrite, izone, lastzonesharing
        character(len=maxStringLen) :: fname
        character(len=7) :: intstring
        integer(kind=intType), dimension(:), allocatable :: nodesizes, nodedisps
        integer(kind=intType), dimension(:), allocatable :: cellsizes, celldisps
        character(len=maxCGNSNameLen), dimension(:), allocatable :: solnames
        real(kind=realtype), dimension(:, :), allocatable :: nodalvalues
        integer(kind=intType), dimension(:, :), allocatable :: conn, localconn
        real(kind=realtype), dimension(:, :), allocatable :: vars
        integer(kind=intType), dimension(:), allocatable :: mask, elemfam, localelemfam, cgnsblockid
        logical :: blanksave, bcgroupneeded, datawritten
        type(zippermesh), pointer :: zipper
        type(familyexchange), pointer :: exch
        if (myid == 0 .and. printiterations) then
            print "(a)", "#"
            print "(a)", "# Writing tecplot surface file(s):"
        end if
 
        ! Number of surface variables. Note that we *explictly*
        ! remove the potential for writing the surface blanks as
        ! these are not necessary for the tecplot IO as we write the
        ! zipper mesh directly. We must save and restore the
        ! variable in case the CGNS otuput still wants to write it.
        blanksave = surfwriteblank
        surfwriteblank = .false.
        call numberofsurfsolvariables(nsolvar)
        allocate (solnames(nsolvar))
        call surfsolnames(solnames)
 
        spectralloop: do sps = 1, ntimeintervalsspectral
 
            ! If it is time spectral we need to agument the filename
            if (equationmode == timespectral) then
                write (intstring, "(i7)") sps
                intstring = adjustl(intstring)
                fname = trim(filename)//"Spectral"//trim(intstring)
            else
                fname = filename
            end if
 
            fileid = 11
            ! Open file on root proc:
 
            if (myid == 0) then
                ! Print the filename to stdout
                print "(a,4x,a)", "#", trim(fname)
 
                open (unit=fileid, file=trim(fname), form='UNFORMATTED', access='stream', status='replace')
 
                ! Tecplot magic number
                write (fileid) "#!TDV112"
 
                ! Integer value of 1 (4 bytes)
                call writeinteger(1)
 
                ! Integer for FileType: 0 = Full, 1= Grid, 2 = Solution
                call writeinteger(0)
 
                ! Write the title of the file
                call writestring("ADflow Surface Solution Data")
 
                ! Write the number of variable names
                call writeinteger(3 + nsolvar)
 
                ! Write the variable names
                call writestring("CoordinateX")
                call writestring("CoordinateY")
                call writestring("CoordinateZ")
 
                ! Write the rest of the variables
                do i = 1, nsolvar
                    call writestring(trim(solnames(i)))
                end do
 
                deallocate (solnames)
 
            end if
 
            ! First pass through to generate and write header information
            masterbcloop1: do ibcgroup = 1, nfamexchange
 
                ! Pointers for easier reading
                exch => bcfamexchange(ibcgroup, sps)
                zipper => zippermeshes(ibcgroup)
 
                ! First thing we do is figure out if we actually need to do
                ! anything with this BCgroup at all. If none the requested
                ! families are in this BCExcahnge we don't have to do
                ! anything.
 
                bcgroupneeded = .false.
                do i = 1, size(famlist)
                    if (faminlist(famlist(i), exch%famList)) then
                        bcgroupneeded = .true.
                    end if
                end do
 
                ! Keep going if we don't need this.
                if (.not. bcgroupneeded) then
                    cycle
                end if
 
                ! Get the sizes of this BCGroup
                call getsurfacesize(sizenode, sizecell, exch%famList, size(exch%famList), .true.)
                call mpi_reduce(sizenode, nnodes, 1, adflow_integer, mpi_sum, 0, adflow_comm_world, ierr)
                call echk(ierr, __file__, __line__)
 
                ! Now gather up the cell family info.
                allocate (celldisps(0:nproc), cellsizes(nproc))
 
                call mpi_gather(sizecell, 1, adflow_integer, &
                                cellsizes, 1, adflow_integer, 0, adflow_comm_world, ierr)
                call echk(ierr, __file__, __line__)
 
                allocate (localelemfam(sizecell))
                call getsurfacefamily(localelemfam, sizecell, exch%famList, size(exch%famList), .true.)
 
                if (myid == 0) then
                    celldisps(0) = 0
                    do iproc = 1, nproc
                        celldisps(iproc) = celldisps(iproc - 1) + cellsizes(iproc)
                    end do
                    ncells = sum(cellsizes)
                    allocate (elemfam(ncells))
                end if
 
                call mpi_gatherv(localelemfam, &
                                 size(localelemfam), adflow_integer, elemfam, &
                                 cellsizes, celldisps, adflow_integer, 0, adflow_comm_world, ierr)
                call echk(ierr, __file__, __line__)
 
                ! Deallocate temp memory
                deallocate (localelemfam, cellsizes, celldisps)
 
                rootproc: if (myid == 0 .and. ncells > 0) then
                    do ifam = 1, size(exch%famList)
 
                        ! Check if we have to write this one:
                        faminclude: if (faminlist(exch%famList(ifam), famlist)) then
                            ncellstowrite = 0
                            do i = 1, ncells
                                ! Check if this elem is to be included
                                if (elemfam(i) == exch%famList(ifam)) then
                                    ncellstowrite = ncellstowrite + 1
                                end if
                            end do
 
                            if (ncellstowrite > 0) then
                                call writefloat(zonemarker) ! Zone Marker
                                call writestring(trim(famnames(exch%famList(ifam)))) ! Zone Name
                                call writeinteger(-1) ! Parent Zone (-1 for None)
                                call writeinteger(-1) ! Strand ID (-2 for tecplot assignment)
                                call writedouble(zero) ! Solution Time
                                call writeinteger(-1)! Zone Color (Not used anymore) (-1)
                                call writeinteger(3) ! Zone Type (3 for FEQuadrilateral)
                                call writeinteger(0)! Data Packing (0 for block)
                                call writeinteger(0)! Specify Var Location (0=don't specify, all at nodes)
                                call writeinteger(0) ! Are raw 1-to-1 face neighbours supplied (0 for false)
                                call writeinteger(nnodes) ! Number of nodes in FE Zone
                                call writeinteger(ncellstowrite) ! Number of elements in FE Zone
                                call writeinteger(0) ! ICellDim, jCellDim, kCellDim (for future use, set to 0)
                                call writeinteger(0)
                                call writeinteger(0)
                                call writeinteger(0) ! Aux data specified (0 for no)
                            end if
                        end if faminclude
                    end do
                end if rootproc
                if (myid == 0) then
                    deallocate (elemfam)
                end if
            end do masterbcloop1
 
            if (myid == 0) then
                call writefloat(datasectionmarker) ! Eohmarker to mark difference between header and data section
            end if
 
            ! Now do everything again but for real.
            masterbcloop: do ibcgroup = 1, nfamexchange
 
                ! Pointers for easier reading
                exch => bcfamexchange(ibcgroup, sps)
                zipper => zippermeshes(ibcgroup)
 
                ! First thing we do is figure out if we actually need to do
                ! anything with this BCgroup at all. If none the requested
                ! families are in this BCExcahnge we don't have to do
                ! anything.
 
                bcgroupneeded = .false.
                do i = 1, size(famlist)
                    if (faminlist(famlist(i), exch%famList)) then
                        bcgroupneeded = .true.
                    end if
                end do
 
                ! Keep going if we don't need this.
                if (.not. bcgroupneeded) then
                    cycle
                end if
 
                ! Get the sizes of this BCGroup
                call getsurfacesize(sizenode, sizecell, exch%famList, size(exch%famList), .true.)
                allocate (nodalvalues(max(sizenode, 1), nsolvar + 3 + 6))
                ! Compute the nodal data
 
                call computesurfaceoutputnodaldata(bcfamexchange(ibcgroup, sps), &
                                                   zippermeshes(ibcgroup), .false., nodalvalues(:, :))
 
                ! Gather up the number of nodes to be set to the root proc:
                allocate (nodesizes(nproc), nodedisps(0:nproc))
                nodesizes = 0
                nodedisps = 0
 
                call mpi_allgather(sizenode, 1, adflow_integer, nodesizes, 1, adflow_integer, &
                                   adflow_comm_world, ierr)
                call echk(ierr, __file__, __line__)
                nodedisps(0) = 0
                do iproc = 1, nproc
                    nodedisps(iproc) = nodedisps(iproc - 1) + nodesizes(iproc)
                end do
 
                isize = 3 + 6 + nsolvar
                if (myid == 0) then
                    nnodes = sum(nodesizes)
                else
                    nnodes = 1
                end if
 
                ! Only root proc actually has any space allocated
                allocate (vars(nnodes, isize))
 
                ! Gather values to the root proc.
                do i = 1, isize
                    call mpi_gatherv(nodalvalues(:, i), sizenode, &
                                     adflow_real, vars(:, i), nodesizes, nodedisps, &
                                     adflow_real, 0, adflow_comm_world, ierr)
                    call echk(ierr, __file__, __line__)
                end do
                deallocate (nodalvalues)
 
                ! Now gather up the connectivity
                allocate (celldisps(0:nproc), cellsizes(nproc))
 
                call mpi_gather(sizecell, 1, adflow_integer, &
                                cellsizes, 1, adflow_integer, 0, adflow_comm_world, ierr)
                call echk(ierr, __file__, __line__)
 
                if (allocated(cgnsblockid)) then
                    deallocate (cgnsblockid)
                end if
 
                allocate (localconn(4, sizecell), localelemfam(sizecell), cgnsblockid(sizecell))
                call getsurfaceconnectivity(localconn, cgnsblockid, sizecell, exch%famList, size(exch%famList), .true.)
                call getsurfacefamily(localelemfam, sizecell, exch%famList, size(exch%famList), .true.)
 
                if (myid == 0) then
                    celldisps(0) = 0
                    do iproc = 1, nproc
                        celldisps(iproc) = celldisps(iproc - 1) + cellsizes(iproc)
                    end do
                    ncells = sum(cellsizes)
                    allocate (conn(4, ncells))
                    allocate (elemfam(ncells))
                end if
 
                ! We offset the conn array by nodeDisps(iProc) which
                ! automagically adjusts the connectivity to account for the
                ! number of nodes from different processors
 
                call mpi_gatherv(localconn + nodedisps(myid), &
                                 4 * size(localconn, 2), adflow_integer, conn, &
                                 cellsizes * 4, celldisps * 4, adflow_integer, 0, adflow_comm_world, ierr)
                call echk(ierr, __file__, __line__)
 
                call mpi_gatherv(localelemfam, &
                                 size(localelemfam), adflow_integer, elemfam, &
                                 cellsizes, celldisps, adflow_integer, 0, adflow_comm_world, ierr)
                call echk(ierr, __file__, __line__)
 
                ! Local values are finished
                deallocate (localconn, localelemfam)
                izone = 0
                rootproc2: if (myid == 0 .and. ncells > 0) then
 
                    ! Need zero based for binary output.
                    conn = conn - 1
 
                    allocate (mask(ncells))
                    datawritten = .false.
                    do ifam = 1, size(exch%famList)
 
                        ! Check if we have to write this one:
                        faminclude2: if (faminlist(exch%famList(ifam), famlist)) then
 
                            ! Create a temporary mask
                            mask = 0
                            ncellstowrite = 0
                            do i = 1, ncells
                                ! Check if this elem is to be included
                                if (elemfam(i) == exch%famList(ifam)) then
                                    mask(i) = 1
                                    ncellstowrite = ncellstowrite + 1
                                end if
                            end do
 
                            actualwrite2: if (ncellstowrite > 0) then
 
                                ! Write Zone Data
                                call writefloat(zonemarker)
 
                                ! Data Type for each variable (1 for float, 2 for double)
                                do i = 1, 3
                                    if (precisionsurfgrid == precisionsingle) then
                                        call writeinteger(1)
                                    else if (precisionsurfgrid == precisiondouble) then
                                        call writeinteger(2)
                                    end if
                                end do
 
                                do i = 1, nsolvar
                                    if (precisionsurfsol == precisionsingle) then
                                        call writeinteger(1)
                                    else if (precisionsurfsol == precisiondouble) then
                                        call writeinteger(2)
                                    end if
                                end do
 
                                call writeinteger(0) ! Has passive variables (0 for no)
 
                                if (.not. datawritten) then
 
                                    ! Save the current 0-based zone so we know which to share with.
                                    lastzonesharing = izone
 
                                    call writeinteger(0) ! Has variable sharing (0 for no)
                                    call writeinteger(-1)  ! Zone based zone number to share connectivity (-1 for no sharing)
 
                                    ! Min/Max Value for coordinates
                                    do j = 1, 3
                                        call writedouble(minval(vars(:, j)))
                                        call writedouble(maxval(vars(:, j)))
                                    end do
 
                                    ! Min/Max Value for solution variables
                                    do j = 1, nsolvar
                                        call writedouble(minval(vars(:, j + 9)))
                                        call writedouble(maxval(vars(:, j + 9)))
                                    end do
 
                                    ! Dump the coordinates
                                    do j = 1, 3
                                        if (precisionsurfgrid == precisionsingle) then
                                            call writefloats(vars(1:nnodes, j))
                                        else if (precisionsurfsol == precisiondouble) then
                                            call writedoubles(vars(1:nnodes, j))
                                        end if
                                    end do
 
                                    ! Dump the solution variables
                                    do j = 1, nsolvar
                                        if (precisionsurfsol == precisionsingle) then
                                            call writefloats(vars(1:nnodes, j + 9))
                                        else if (precisionsurfsol == precisiondouble) then
                                            call writedoubles(vars(1:nnodes, j + 9))
                                        end if
                                    end do
 
                                else
                                    ! This will be sharing data from another zone.
                                    call writeinteger(1) ! Has variable sharing (0 for no)
 
                                    ! Write out the zone number for sharing the data.
                                    do j = 1, 3 + nsolvar
                                        call writeinteger(lastzonesharing)
                                    end do
 
                                    call writeinteger(-1)  ! Zone based zone number to share connectivity (-1 for no sharing
                                end if
 
                                ! Dump the connectivity
                                j = 0
                                do i = 1, ncells
                                    ! Check if this elem is to be included
                                    if (mask(i) == 1) then
                                        call writeintegers(conn(:, i))
                                        j = j + 1
                                    end if
                                end do
 
                                izone = izone + 1
                            end if actualwrite2
                        end if faminclude2
                    end do
                    deallocate (mask)
                end if rootproc2
                if (myid == 0) then
                    deallocate (conn, elemfam)
                end if
                deallocate (cellsizes, celldisps, nodesizes, nodedisps, vars)
            end do masterbcloop
 
            if (myid == 0) then
                close (fileid)
            end if
 
        end do spectralloop
 
        ! Restore the modified option
        surfwriteblank = blanksave
        if (myid == 0 .and. printiterations) then
            print "(a)", "# Tecplot surface file(s) written"
            print "(a)", "#"
        end if
 
    contains
 
        subroutine writefloat(adflowRealVal)
            use iso_fortran_env, only: real32
            implicit none
            real(kind=realtype) :: adflowrealval
            real(kind=real32) :: float
            float = adflowrealval
            write (fileid) float
        end subroutine writefloat
 
        subroutine writedouble(adflowRealVal)
            use iso_fortran_env, only: real64
            implicit none
            real(kind=realtype) :: adflowrealval
            real(kind=real64) :: dble
            dble = adflowrealval
            write (fileid) dble
        end subroutine writedouble
 
        subroutine writefloats(adflowRealVals)
            use iso_fortran_env, only: real32
            implicit none
            real(kind=realtype) :: adflowrealvals(:)
            real(kind=real32) :: floats(size(adflowrealvals))
            integer :: i
            floats = adflowrealvals
            write (fileid) floats
 
        end subroutine writefloats
 
        subroutine writedoubles(adflowRealVals)
            use iso_fortran_env, only: real64
            implicit none
            real(kind=realtype) :: adflowrealvals(:)
            real(kind=real64) :: dbles(size(adflowrealvals))
            integer :: i
            dbles = adflowrealvals
            write (fileid) dbles
 
        end subroutine writedoubles
 
        subroutine writeinteger(adflowIntegerVal)
            use iso_fortran_env, only: int32
            implicit none
            integer(kind=intType) :: adflowIntegerVal
            integer(kind=int32) :: int
 
            int = adflowintegerval
            write (fileid) int
        end subroutine writeinteger
 
        subroutine writeintegers(adflowIntegerVals)
            use iso_fortran_env, only: int32
            implicit none
            integer(kind=intType) :: adflowIntegerVals(:), i
            integer(kind=int32) :: ints(size(adflowintegervals))
            ints = adflowintegervals
            write (fileid) ints
 
        end subroutine writeintegers
 
        subroutine writestring(str)
 
            implicit none
 
            character(len=*) :: str
            integer(kind=intType) :: i
 
            do i = 1, len(str)
                write (fileid) iachar(str(i:i))
            end do
            write (fileid) 0
 
        end subroutine writestring
 
    end subroutine writetecplotsurfacefile
 
    subroutine initializeliftdistributiondata
        use constants
        use communication
        use blockpointers
        use inputphysics
        use outputmod
        use inputtimespectral
        use surfacefamilies
        implicit none
 
        ! Working Variables
        integer(kind=intType) :: nPts, nCells, sps
 
        if (liftdistinitialized) then
            return
        else
 
            ! Data for the marching squares method: Which edges are cut by
            ! the contour. We haven't dealt with the case of an ambiguous
            ! contour which is case 6 and 11. Since most of the slices we are
            ! doing are planes this won't matter.
            mscon1(1, :) = (/0, 0, 0, 0, 0/)
            mscon1(2, :) = (/1, 4, 0, 0, 0/)
            mscon1(3, :) = (/1, 2, 0, 0, 0/)
            mscon1(4, :) = (/4, 2, 0, 0, 0/)
            mscon1(5, :) = (/2, 3, 0, 0, 0/)
            mscon1(6, :) = (/2, 3, 0, 0, 0/) ! Should be 2, 3, 1, 4
            mscon1(7, :) = (/1, 3, 0, 0, 0/)
            mscon1(8, :) = (/4, 3, 0, 0, 0/)
            mscon1(9, :) = (/4, 3, 0, 0, 0/)
            mscon1(10, :) = (/1, 3, 0, 0, 0/)
            mscon1(11, :) = (/2, 3, 0, 0, 0/) ! Should be 2, 3, 1, 4
            mscon1(12, :) = (/2, 3, 0, 0, 0/)
            mscon1(13, :) = (/4, 2, 0, 0, 0/)
            mscon1(14, :) = (/1, 2, 0, 0, 0/)
            mscon1(15, :) = (/1, 4, 0, 0, 0/)
            mscon1(16, :) = (/0, 0, 0, 0, 0/)
 
            mscon2(1, :) = (/1, 2/)
            mscon2(2, :) = (/2, 3/)
            mscon2(3, :) = (/3, 4/)
            mscon2(4, :) = (/4, 1/)
            liftdistinitialized = .true.
        end if
 
    end subroutine initializeliftdistributiondata
 
    subroutine computesurfaceoutputnodaldata(exch, zipper, includeTractions, nodalValues)
        !
        !       This purpose of this subroutine is to compute all nodal values
        !
        use constants
        use communication
        use inputphysics
        use blockpointers
        use surfacefamilies, only: bcfamgroups, familyexchange
        use outputmod, only: storesurfsolinbuffer, numberofsurfsolvariables, &
                             surfsolnames
        use surfaceutils
        use utils, only: setpointers, echk
        use sorting, only: faminlist
        use oversetdata, only: zippermesh
#include <petsc/finclude/petsc.h>
        use petsc
        implicit none
        ! Input Param
        type(familyexchange) :: exch
        logical :: includeTractions
        real(kind=realtype), dimension(:, :), intent(inout) :: nodalvalues
        type(zippermesh) :: zipper
        ! Working params
        integer(kind=intType) :: i, j, ii, jj, kk, nn, mm, isol, ierr, npts, ncells
        integer(kind=intType) :: nfields, nsolvar, ibeg, iend, jbeg, jend, ind(4), ni, nj
        integer(kind=intType) :: sizeNode, sizeCell, iDim
        integer(kind=intType), dimension(3, 2) :: cellRangeCGNS
        character(len=maxCGNSNameLen), dimension(:), allocatable :: solNames
        real(kind=realtype), dimension(:), allocatable :: buffer
        real(kind=realtype), dimension(:), pointer :: weightptr, localptr
        real(kind=realtype), dimension(:, :), allocatable :: tmp
        logical :: viscousSubFace
 
        nodalvalues = zero
 
        call numberofsurfsolvariables(nsolvar)
        allocate (solnames(nsolvar))
        call surfsolnames(solnames)
 
        ! The tractions have a sepcial routine so call that first before we
        ! mess with the family information.
 
        if (includetractions) then
 
            ! And put the traction where it needs to go if necessary.  Note that
            ! this computation will only work corectly when the wallExchange
            ! which include all wallgroups is used.
 
            call computenodaltractions(exch%sps)
            ii = 0
            do nn = 1, ndom
                call setpointers(nn, 1_inttype, exch%sps)
                do mm = 1, nbocos
                    ibeg = bcdata(mm)%inBeg; iend = bcdata(mm)%inEnd
                    jbeg = bcdata(mm)%jnBeg; jend = bcdata(mm)%jnEnd
 
                    if (faminlist(bcdata(mm)%famID, exch%famList)) then
                        do j = jbeg, jend
                            do i = ibeg, iend
                                ii = ii + 1
                                nodalvalues(ii, 4) = bcdata(mm)%Tp(i, j, 1)
                                nodalvalues(ii, 5) = bcdata(mm)%Tp(i, j, 2)
                                nodalvalues(ii, 6) = bcdata(mm)%Tp(i, j, 3)
 
                                nodalvalues(ii, 7) = bcdata(mm)%Tv(i, j, 1)
                                nodalvalues(ii, 8) = bcdata(mm)%Tv(i, j, 2)
                                nodalvalues(ii, 9) = bcdata(mm)%Tv(i, j, 3)
                            end do
                        end do
                    end if
                end do
            end do
 
            ! Not quite dont yet with the nodal tractions; we need to send
            ! the nodal tractions that the duplices the *zipper* mesh needs
            ! to the root proc. We have a special scatter for this.
 
            if (zipper%allocated) then
 
                ! Loop over the 6 tractions
                do idim = 1, 6
 
                    ! Copy the values into localPtr
                    call vecgetarrayf90(exch%nodeValLocal, localptr, ierr)
                    call echk(ierr, __file__, __line__)
 
                    do i = 1, exch%nNodes
                        localptr(i) = nodalvalues(i, idim + 3)
                    end do
 
                    call vecrestorearrayf90(exch%nodeValLocal, localptr, ierr)
                    call echk(ierr, __file__, __line__)
 
                    ! Now use the *zipper* scatter
                    call vecscatterbegin(zipper%scatter, exch%nodeValLocal, &
                                         zipper%localVal, insert_values, scatter_forward, ierr)
                    call echk(ierr, __file__, __line__)
 
                    call vecscatterend(zipper%scatter, exch%nodeValLocal, &
                                       zipper%localVal, insert_values, scatter_forward, ierr)
                    call echk(ierr, __file__, __line__)
 
                    ! Copy the zipper values out on the root proc
                    if (myid == 0) then
                        call vecgetarrayf90(zipper%localVal, localptr, ierr)
                        call echk(ierr, __file__, __line__)
 
                        do i = 1, size(localptr)
                            nodalvalues(exch%nNodes + i, idim + 3) = localptr(i)
                        end do
 
                        call vecrestorearrayf90(zipper%localVal, localptr, ierr)
                        call echk(ierr, __file__, __line__)
                    end if
                end do
            end if
        end if
 
        ! Get the current set of surface points for the family we just set.
        call getsurfacesize(sizenode, sizecell, exch%famList, size(exch%famlist), .true.)
        allocate (tmp(3, sizenode))
        call getsurfacepoints(tmp, sizenode, exch%sps, exch%famList, size(exch%famList), .true.)
 
        do i = 1, sizenode
            nodalvalues(i, 1:3) = tmp(1:3, i)
        end do
        deallocate (tmp)
        ! For the remainder of the variables, use arithematic averaging.
        call vecgetarrayf90(exch%nodeValLocal, localptr, ierr)
        call echk(ierr, __file__, __line__)
 
        localptr = zero
        ! ii is the index into the pointer array
        ii = 0
        do nn = 1, ndom
            call setpointers(nn, 1_inttype, exch%sps)
            do mm = 1, nbocos
                ibeg = bcdata(mm)%inBeg; iend = bcdata(mm)%inEnd
                jbeg = bcdata(mm)%jnBeg; jend = bcdata(mm)%jnEnd
                ni = iend - ibeg + 1
                nj = jend - jbeg + 1
                if (faminlist(bcdata(mm)%famID, exch%FamList)) then
                    do j = 0, nj - 2
                        do i = 0, ni - 2
 
                            ! Scatter 1 to each node.
                            ind(1) = ii + (j) * ni + i + 1
                            ind(2) = ii + (j) * ni + i + 2
                            ind(3) = ii + (j + 1) * ni + i + 2
                            ind(4) = ii + (j + 1) * ni + i + 1
                            do jj = 1, 4
                                localptr(ind(jj)) = localptr(ind(jj)) + one
                            end do
                        end do
                    end do
                    ii = ii + ni * nj
                end if
            end do
        end do
 
        call vecrestorearrayf90(exch%nodeValLocal, localptr, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Globalize the area
        call vecset(exch%sumGlobal, zero, ierr)
        call echk(ierr, __file__, __line__)
 
        call vecscatterbegin(exch%scatter, exch%nodeValLocal, &
                             exch%sumGlobal, add_values, scatter_forward, ierr)
        call echk(ierr, __file__, __line__)
 
        call vecscatterend(exch%scatter, exch%nodeValLocal, &
                           exch%sumGlobal, add_values, scatter_forward, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Now compute the inverse of the weighting so that we can multiply
        ! instead of dividing.
 
        call vecgetarrayf90(exch%sumGlobal, localptr, ierr)
        call echk(ierr, __file__, __line__)
 
        localptr = one / localptr
 
        call vecrestorearrayf90(exch%sumGlobal, localptr, ierr)
        call echk(ierr, __file__, __line__)
 
        varloop: do isol = 1, nsolvar
 
            ! Extract the poitner to the local array
            call vecgetarrayf90(exch%nodeValLocal, localptr, ierr)
            call echk(ierr, __file__, __line__)
 
            ! ii is the continous running element pointer
            ii = 0
            localptr = zero
            ! Do each variable separately.
            domainloop: do nn = 1, ndom
                call setpointers(nn, 1, exch%sps)
 
                bocoloop: do mm = 1, nbocos
                    if (faminlist(bcdata(mm)%famID, exch%famList)) then
 
                        ! storeSurfSolInBuffer needs to know if the subface is
                        ! viscous or not.
                        viscoussubface = .false.
                        if (bctype(mm) == nswalladiabatic .or. &
                            bctype(mm) == nswallisothermal) then
                            viscoussubface = .true.
                        end if
 
                        ! Determine the cell range *in the original cgns
                        ! ordering*. The reason for this is the storeSurfSolInBufer
                        ! is normally used for writing CGNS and thus it is that
                        ! ording that is important. However, that isn't too hard to
                        ! deal with.
 
                        jbeg = bcdata(mm)%jnBeg + 1
                        jend = bcdata(mm)%jnEnd
                        ibeg = bcdata(mm)%inBeg + 1
                        iend = bcdata(mm)%inEnd
 
                        ! Dummy value of 1 for the face values not set.
                        cellrangecgns = 1
                        select case (bcfaceid(mm))
                        case (imin, imax)
                            ! I range meaningless
                            cellrangecgns(2, 1) = ibeg + jbegor - 1
                            cellrangecgns(2, 2) = iend + jbegor - 1
 
                            cellrangecgns(3, 1) = jbeg + kbegor - 1
                            cellrangecgns(3, 2) = jend + kbegor - 1
 
                        case (jmin, jmax)
                            ! J range meaningless
                            cellrangecgns(1, 1) = ibeg + ibegor - 1
                            cellrangecgns(1, 2) = iend + ibegor - 1
 
                            cellrangecgns(3, 1) = jbeg + kbegor - 1
                            cellrangecgns(3, 2) = jend + kbegor - 1
 
                        case (kmin, kmax)
                            ! J range meaningless
                            cellrangecgns(1, 1) = ibeg + ibegor - 1
                            cellrangecgns(1, 2) = iend + ibegor - 1
 
                            cellrangecgns(2, 1) = jbeg + jbegor - 1
                            cellrangecgns(2, 2) = jend + jbegor - 1
                        end select
 
                        ! Allocate enough space for the 1D buffer
                        allocate (buffer((iend - ibeg + 1) * (jend - jbeg + 1)))
                        kk = 0
                        call storesurfsolinbuffer(exch%sps, buffer, kk, nn, bcfaceid(mm), &
                                                  cellrangecgns, solnames(isol), viscoussubface, .false.)
 
                        ! Now since the storeSurfSol just put things in a flat
                        ! array and are face based, here we take the 1D face data
                        ! and scatter to the nodes.
 
                        ibeg = bcdata(mm)%inBeg; iend = bcdata(mm)%inEnd
                        jbeg = bcdata(mm)%jnBeg; jend = bcdata(mm)%jnEnd
                        ni = iend - ibeg + 1
                        nj = jend - jbeg + 1
                        jj = 0
                        do j = 0, nj - 2
                            do i = 0, ni - 2
                                jj = jj + 1
                                ! Scatter value to each node
                                ind(1) = ii + (j) * ni + i + 1
                                ind(2) = ii + (j) * ni + i + 2
                                ind(3) = ii + (j + 1) * ni + i + 2
                                ind(4) = ii + (j + 1) * ni + i + 1
                                do kk = 1, 4
                                    localptr(ind(kk)) = localptr(ind(kk)) + buffer(jj)
                                end do
                            end do
                        end do
 
                        ii = ii + ni * nj
                        deallocate (buffer)
                    end if
                end do bocoloop
            end do domainloop
 
            ! Return our pointer
            call vecrestorearrayf90(exch%nodeValLocal, localptr, ierr)
            call echk(ierr, __file__, __line__)
 
            ! Globalize the function value
            call vecset(exch%nodeValGlobal, zero, ierr)
            call echk(ierr, __file__, __line__)
 
            call vecscatterbegin(exch%scatter, exch%nodeValLocal, &
                                 exch%nodeValGlobal, add_values, scatter_forward, ierr)
            call echk(ierr, __file__, __line__)
 
            call vecscatterend(exch%scatter, exch%nodeValLocal, &
                               exch%nodeValGlobal, add_values, scatter_forward, ierr)
            call echk(ierr, __file__, __line__)
 
            ! Now divide by the weighting sum. We can do this with a
            ! vecpointwisemult since we already divided by the weight.
            call vecpointwisemult(exch%nodeValGlobal, exch%nodeValGlobal, &
                                  exch%sumGlobal, ierr)
            call echk(ierr, __file__, __line__)
 
            ! Push back to the local values
            call vecscatterbegin(exch%scatter, exch%nodeValGlobal, &
                                 exch%nodeValLocal, insert_values, scatter_reverse, ierr)
            call echk(ierr, __file__, __line__)
 
            call vecscatterend(exch%scatter, exch%nodeValGlobal, &
                               exch%nodeValLocal, insert_values, scatter_reverse, ierr)
            call echk(ierr, __file__, __line__)
 
            ! Copy the values into nodalValues
            call vecgetarrayf90(exch%nodeValLocal, localptr, ierr)
            call echk(ierr, __file__, __line__)
 
            do i = 1, size(localptr)
                nodalvalues(i, isol + 9) = localptr(i)
            end do
 
            call vecrestorearrayf90(exch%nodeValLocal, localptr, ierr)
            call echk(ierr, __file__, __line__)
 
            if (zipper%allocated) then
 
                ! Now use the *zipper* scatter
                call vecscatterbegin(zipper%scatter, exch%nodeValLocal, &
                                     zipper%localVal, insert_values, scatter_forward, ierr)
                call echk(ierr, __file__, __line__)
 
                call vecscatterend(zipper%scatter, exch%nodeValLocal, &
                                   zipper%localVal, insert_values, scatter_forward, ierr)
                call echk(ierr, __file__, __line__)
 
                ! Copy the zipper values out on the root proc
                if (myid == 0) then
                    call vecgetarrayf90(zipper%localVal, localptr, ierr)
                    call echk(ierr, __file__, __line__)
 
                    do i = 1, size(localptr)
                        nodalvalues(exch%nNodes + i, 9 + isol) = localptr(i)
                    end do
 
                    call vecrestorearrayf90(zipper%localVal, localptr, ierr)
                    call echk(ierr, __file__, __line__)
                end if
            end if
        end do varloop
        deallocate (solnames)
 
    end subroutine computesurfaceoutputnodaldata
 
    subroutine createslice(pts, conn, elemFam, slc, pt, normal, dir_vec, use_dir, sliceName, famList)
        !
        !       This subroutine creates a slice on a plane defined by pt and
        !       and dir. It only uses the families specified in the famList.
        !       sps define which specral instance to use.
        !
        use constants
        use utils, only: reallocatereal2, reallocateinteger2, pointreduce
        use sorting, only: faminlist
        implicit none
 
        ! Input param
        real(kind=realtype), dimension(:, :), intent(in) :: pts
        integer(kind=intType), dimension(:, :), intent(in) :: conn
        integer(kind=intType), dimension(:), intent(in) :: elemFam
        type(slice), intent(inout) :: slc
        real(kind=realtype), dimension(3), intent(in) :: pt, dir_vec, normal
        logical, intent(in) :: use_dir
        character(len=*), intent(in) :: sliceName
        integer(kind=intType), dimension(:), intent(in) :: famList
 
        ! Working param
        integer(kind=intType) :: i, nMax, nUnique, oldInd, newInd
        integer(kind=intType) :: patchIndices(4), indexSquare, jj, kk, icon, iCoor, num1, num2
        real(kind=realtype) :: f(4), d, ovrdnom, tol, len_vec
        logical :: logic1, foundFam
        real(kind=realtype), dimension(:, :), pointer :: tmpweight, dummy, tmpnodes
        integer(kind=intType), dimension(:, :), pointer :: tmpInd
        integer(kind=intType), dimension(:), allocatable :: link
        real(kind=realtype), dimension(:), allocatable :: fc
        real(kind=realtype), dimension(3) :: elemc, vec
 
        ! Allocate the family list this slice is to use:
        slc%sliceName = slicename
        ! Set the info for the slice:
        slc%pt = pt
        slc%normal = normal
        slc%dir_vec = dir_vec
        slc%use_dir = use_dir
        slc%nNodes = 0
        allocate (slc%famList(size(famlist)))
        slc%famList = famlist
        ! First step is to compute the 'function' value that will be used
        ! for the contour.
 
        ! Equation of plane: ax + by + cz + d = 0
        d = -pt(1) * normal(1) - pt(2) * normal(2) - pt(3) * normal(3)
        ovrdnom = one / sqrt(normal(1)**2 + normal(2)**2 + normal(3)**2)
 
        ! Compute the distance function on all possible surfaces on this
        ! processor.
        allocate (fc(size(pts, 2)))
        do i = 1, size(pts, 2)
            ! Now compute the signed distance
            fc(i) = (normal(1) * pts(1, i) + normal(2) * pts(2, i) + normal(3) * pts(3, i) + d) * ovrdnom
        end do
 
        ! Estimate size of slice by the sqrt of the number of nodes in the
        ! mesh. Exact size doesn't matter as we realloc if necessary.
        nmax = int(sqrt(dble(size(pts, 2))))
        allocate (tmpweight(2, nmax), tmpind(2, nmax), tmpnodes(3, nmax))
 
        icoor = 0
        oldind = 1
 
        ! Loop over all elements
        elemloop: do i = 1, size(conn, 2)
 
            faminclude: if (faminlist(elemfam(i), famlist)) then
 
                ! zero out the centroid
                elemc = zero
 
                ! Extract the indices and function values at each corner
                do jj = 1, 4
                    patchindices(jj) = conn(jj, i)
                    f(jj) = fc(patchindices(jj))
                    elemc = elemc + 0.25_realtype * pts(:, patchindices(jj))
                end do
 
                ! check if we are using the direction to pick sliced elements
                if (use_dir) then
                    ! check if the centroid of this element is in the correct direction
                    vec = elemc - pt
                    ! normalize vector
                    len_vec = sqrt(vec(1) * vec(1) + vec(2) * vec(2) + vec(3) * vec(3))
                    vec = vec / len_vec
                    if ((vec(1) * dir_vec(1) + vec(2) * dir_vec(2) + vec(3) * dir_vec(3)) .lt. zero) then
                        ! we reject this element; just set all signed distances to 1.0
                        do jj = 1, 4
                            f(jj) = one
                        end do
                    end if
                end if
 
                ! Based on the values at each corner, determine which
                ! type contour we have
                indexsquare = 1
 
                if (f(1) .lt. zero) indexsquare = indexsquare + 1
                if (f(2) .lt. zero) indexsquare = indexsquare + 2
                if (f(3) .lt. zero) indexsquare = indexsquare + 4
                if (f(4) .lt. zero) indexsquare = indexsquare + 8
 
                logic1 = .true.
 
                kk = 1
                do while (logic1)
                    ! This is the edge
                    icon = mscon1(indexsquare, kk)
 
                    if (icon == 0) then
                        logic1 = .false.
                    else
 
                        ! num1, num2 are node indices
                        num1 = mscon2(icon, 1)
                        num2 = mscon2(icon, 2)
 
                        icoor = icoor + 1
                        if (icoor > nmax) then
                            ! Need to reallocate the arrays. Make it double the size
                            call reallocatereal2(tmpweight, 2, 2 * nmax, 2, nmax, .true.)
                            call reallocatereal2(tmpnodes, 3, 2 * nmax, 3, nmax, .true.)
                            call reallocateinteger2(tmpind, 2, 2 * nmax, 2, nmax, .true.)
                            nmax = nmax * 2
                        end if
 
                        ! Weight factors
                        tmpweight(2, icoor) = (zero - f(num1)) / (f(num2) - f(num1))
                        tmpweight(1, icoor) = one - tmpweight(2, icoor)
 
                        ! Store the weight factors
                        tmpind(:, icoor) = (/patchindices(num1), patchindices(num2)/)
 
                        ! Store the physical nodes so we know how to reduce
                        tmpnodes(:, icoor) = &
                            tmpweight(1, icoor) * pts(:, tmpind(1, icoor)) + &
                            tmpweight(2, icoor) * pts(:, tmpind(2, icoor))
                        kk = kk + 1
                    end if
                end do
            end if faminclude
        end do elemloop
 
        ! To save disk space, we can compact out the doubly defined nodes that
        ! were created during the slicing process. Then we can allocate the
        ! final weight array and index array to be the exact correct
        ! length
 
        allocate (dummy(3, icoor), link(icoor))
        tol = 1e-12
        call pointreduce(tmpnodes, icoor, tol, dummy, link, nunique)
        allocate (slc%w(2, nunique), slc%ind(2, nunique), slc%conn(2, icoor / 2))
        slc%nNodes = nunique
 
        ! Modify the data accordingly
        do i = 1, icoor
            slc%w(:, link(i)) = tmpweight(:, i)
            slc%ind(:, link(i)) = tmpind(:, i)
        end do
 
        ! The connectivity is actually link reshaped since the original conn
        ! would have been (1,2), (3,4), (5,6) etc.
        do i = 1, icoor / 2
            slc%conn(1, i) = link(2 * i - 1)
            slc%conn(2, i) = link(2 * i)
        end do
 
        deallocate (tmpnodes, tmpweight, tmpind, dummy, link, fc)
    end subroutine createslice
 
    subroutine destroyslice(slc)
        !
        !       This subroutine destroys a slice created by the createSlice
        !       routine
        !
        use constants
        implicit none
 
        ! Input param
        type(slice), intent(inout) :: slc
 
        ! Deallocate weights and indices if they are already allocated
        if (allocated(slc%w)) then
            deallocate (slc%w)
        end if
 
        if (allocated(slc%ind)) then
            deallocate (slc%ind)
        end if
 
        if (allocated(slc%conn)) then
            deallocate (slc%conn)
        end if
 
        if (allocated(slc%famList)) then
            deallocate (slc%famList)
        end if
 
        if (allocated(slc%vars)) then
            deallocate (slc%vars)
        end if
 
    end subroutine destroyslice
 
    subroutine integrateslice(lSlc, gSlc, nodalValues, nFields, doConnectivity)
        !
        !       This subroutine integrates the forces on slice slc and computes
        !       the integrated quantities for lift, drag, cl and cd with
        !       contributions from both the pressure and visoucs forces.
        !       It optionally interpolates solution variables as well.
        !
        use constants
        use inputphysics
        use flowvarrefstate
        use communication
        use utils, only: echk
        implicit none
 
        ! Input Variables
        type(slice) :: lSlc, gSlc
        integer(kind=intType), intent(in) :: nFields
        logical, intent(in) :: doConnectivity
        real(kind=realtype), dimension(:, :), intent(in) :: nodalvalues
        ! Working variables
        integer(kind=intType) :: i, j, i1, i2
        real(kind=realtype), dimension(3) :: x1, x2, pt1, pt2, vt1, vt2, pf, vf, pf_elem, vf_elem
        real(kind=realtype) :: len, dmax, dmin, dist, fact, m(3, 3), tmp(6)
        real(kind=realtype) :: r(3), r_new(3), hyp, te(3), le(3), theta, w1, w2
        integer(kind=intType) :: bestPair(2), normal_ind, iProc, ierr, iSize
        real(kind=realtype), dimension(:, :), allocatable :: tempcoords
        real(kind=realtype), dimension(:, :), allocatable :: localvals
        integer(kind=intType), dimension(:), allocatable :: sliceNodeSizes, sliceCellSizes
        integer(kind=intType), dimension(:), allocatable :: nodeDisps, cellDisps
        real(kind=realtype), dimension(3) :: refpoint, pm, vm
        real(kind=realtype) :: xc, yc, zc
        ! Need another tmp to reduce the fx, fy and fz values
        real(kind=realtype) :: tmp_(3)
 
        ! Copy the info related to slice
        gslc%sliceName = trim(lslc%sliceName)
        gslc%pt = lslc%pt
        gslc%normal = lslc%normal
 
        ! Back out what is the main index of the slice, x, y or z based on
        ! the direction. Not the best approach, but that's ok
        ! TODO this can be improved since we are now doing arbitrary slice directions
        normal_ind = maxloc(abs(gslc%normal), 1)
 
        pf = zero
        vf = zero
        pm = zero  ! pressure moment
        vm = zero  ! viscous moment
        isize = 3 + 6 + nfields
        allocate (localvals(isize, lslc%nNodes))
 
        ! Interpolate the required values
        do i = 1, lslc%nNodes
            i1 = lslc%ind(1, i)
            i2 = lslc%ind(2, i)
            w1 = lslc%w(1, i)
            w2 = lslc%w(2, i)
            localvals(1:isize, i) = w1 * nodalvalues(i1, 1:isize) + w2 * nodalvalues(i2, 1:isize)
        end do
 
        ! first communicate the coordinates and connectivities to get the chord and the LE/TE
        ! we need this information to get the quarter-chord location,
        ! which is needed for the moment computations
 
        ! Gather up the number of nodes to be set to the root proc:
        allocate (slicenodesizes(nproc), nodedisps(0:nproc))
        slicenodesizes = 0
        nodedisps = 0
        call mpi_allgather(lslc%nNodes, 1, adflow_integer, slicenodesizes, 1, adflow_integer, &
                           adflow_comm_world, ierr)
        call echk(ierr, __file__, __line__)
        nodedisps(0) = 0
        do iproc = 1, nproc
            nodedisps(iproc) = nodedisps(iproc - 1) + slicenodesizes(iproc) * isize
        end do
 
        if (myid == 0) then
            gslc%nNodes = sum(slicenodesizes)
            allocate (gslc%vars(isize, gslc%nNodes))
        end if
 
        call mpi_gatherv(localvals, isize * lslc%nNodes, adflow_real, gslc%vars, slicenodesizes * isize, &
                         nodedisps, adflow_real, 0, adflow_comm_world, ierr)
        call echk(ierr, __file__, __line__)
 
        ! We may also need to gather the connectivity if the slice will have
        ! to be written to a file.
        if (doconnectivity) then
            i = size(lslc%conn, 2)
            allocate (celldisps(0:nproc), slicecellsizes(nproc))
 
            call mpi_gather(i, 1, adflow_integer, slicecellsizes, 1, adflow_integer, &
                            0, adflow_comm_world, ierr)
            call echk(ierr, __file__, __line__)
 
            if (myid == 0) then
                celldisps(0) = 0
                do iproc = 1, nproc
                    celldisps(iproc) = celldisps(iproc - 1) + slicecellsizes(iproc) * 2
                end do
                allocate (gslc%conn(2, sum(slicecellsizes)))
            end if
 
            ! We offset the conn array by nodeDisps(iProc) which
            ! automagically adjust the connectivity to account for the
            ! number of nodes from different processors
 
            call mpi_gatherv(lslc%conn + nodedisps(myid) / isize, 2 * size(lslc%conn, 2), adflow_integer, gslc%conn, &
                             slicecellsizes * 2, celldisps, adflow_integer, 0, adflow_comm_world, ierr)
            call echk(ierr, __file__, __line__)
 
            ! Not quite finished yet since we will have gathered nodes from
            ! multiple procs we have to adjust the connectivity
 
            deallocate (slicecellsizes, celldisps)
        end if
 
        deallocate (slicenodesizes, nodedisps)
        ! done communicating the geometry info
 
        ! process the geometry on the root processor to figure out the LE, TE, and the chord
        if (myid == 0) then
            ! TODO the comments below can be used for a _slightly_ more accurate LE/TE detection.
            ! this does not really show any difference other than the twist computation,
            ! but depending on the application, getting the twist distribution from ADflow data might be critical.
            ! so we leave the comments below for a better implementation.
            ! For now, we just work with the 2 max-distance nodes
 
            ! begin template:
            ! loop over the elements and find the TE bends
            ! if we have 2 TE bends, then take the mid-TE point as the TE
            ! find the largest distance node
            ! if we dont have 2 bends, just get the max distance between any 2 nodes, thats our chord line
            ! end template.
 
            ! Compute the chord as the max length between any two nodes...this
            ! is n^2, so should be changed in the future
            dmax = zero
            bestpair = (/1, 1/)
            do i = 1, size(gslc%vars, 2)
                ! extract node:
                x1 = gslc%vars(1:3, i)
 
                do j = i + 1, size(gslc%vars, 2)
                    ! extract node:
                    x2 = gslc%vars(1:3, j)
 
                    dist = sqrt((x1(1) - x2(1))**2 + (x1(2) - x2(2))**2 + (x1(3) - x2(3))**2)
 
                    if (dist > dmax) then
                        dmax = dist
                        bestpair = (/i, j/)
                    end if
                end do
            end do
 
            ! Set chord, protected from zero
            gslc%chord = max(dmax, 1e-12)
 
            ! figure out which node is the LE and the TE
            if (gslc%vars(1, bestpair(1)) < gslc%vars(1, bestpair(2))) then
                ! first entry in the "bestPair" is the LE node
                x1 = gslc%vars(1:3, bestpair(1))
                x2 = gslc%vars(1:3, bestpair(2))
            else
                ! second entry in the "bestPair" is the LE node
                x1 = gslc%vars(1:3, bestpair(2))
                x2 = gslc%vars(1:3, bestpair(1))
            end if
 
            ! using the LE and TE coordinates, compute the quarter chord point
            ! we will use this as reference for the moment computation
            refpoint(1) = 0.75_realtype * x1(1) + 0.25_realtype * x2(1)
            refpoint(2) = 0.75_realtype * x1(2) + 0.25_realtype * x2(2)
            refpoint(3) = 0.75_realtype * x1(3) + 0.25_realtype * x2(3)
        end if
 
        ! communicate the reference point coordinates across processors.
        call mpi_bcast(refpoint, 3, adflow_real, 0, adflow_comm_world, ierr)
        call echk(ierr, __file__, __line__)
 
        ! now we are done with the global geometric operations. go back to integrating the slices locally on each proc
 
        ! loop over elements
        do i = 1, size(lslc%conn, 2)
            ! Compute all the local variables we need.
 
            ! extract nodes:
            i1 = lslc%conn(1, i)
            i2 = lslc%conn(2, i)
 
            x1 = localvals(1:3, i1)
            x2 = localvals(1:3, i2)
 
            ! extract pressure tractions
            pt1 = localvals(4:6, i1)
            pt2 = localvals(4:6, i2)
 
            ! extract viscous tractions
            vt1 = localvals(7:9, i1)
            vt2 = localvals(7:9, i2)
 
            ! Length of this segment
            len = sqrt((x1(1) - x2(1))**2 + (x1(2) - x2(2))**2 + (x1(3) - x2(3))**2)
 
            ! compute the pressure and viscous forces on this element
            pf_elem = half * (pt1 + pt2) * len
            vf_elem = half * (vt1 + vt2) * len
 
            ! Integrate the pressure and viscous forces separately
            pf = pf + pf_elem
            vf = vf + vf_elem
 
            ! compute moment about the global reference locations
            xc = half * (x1(1) + x2(1)) - refpoint(1)
            yc = half * (x1(2) + x2(2)) - refpoint(2)
            zc = half * (x1(3) + x2(3)) - refpoint(3)
 
            ! pressure components
            pm(1) = pm(1) + yc * pf_elem(3) - zc * pf_elem(2)
            pm(2) = pm(2) + zc * pf_elem(1) - xc * pf_elem(3)
            pm(3) = pm(3) + xc * pf_elem(2) - yc * pf_elem(1)
 
            ! viscous components
            vm(1) = vm(1) + yc * vf_elem(3) - zc * vf_elem(2)
            vm(2) = vm(2) + zc * vf_elem(1) - xc * vf_elem(3)
            vm(3) = vm(3) + xc * vf_elem(2) - yc * vf_elem(1)
 
        end do
 
        ! That is as far as we can go in parallel. We now have to gather up
        ! pL, pD, pM, vL, vD, vM as well as the nodes to the root proc.
 
        ! Set the local values we can in the slice
        lslc%pL = liftdirection(1) * pf(1) + liftdirection(2) * pf(2) + liftdirection(3) * pf(3)
        lslc%pD = dragdirection(1) * pf(1) + dragdirection(2) * pf(2) + dragdirection(3) * pf(3)
        lslc%vL = liftdirection(1) * vf(1) + liftdirection(2) * vf(2) + liftdirection(3) * vf(3)
        lslc%vD = dragdirection(1) * vf(1) + dragdirection(2) * vf(2) + dragdirection(3) * vf(3)
 
        ! the moments are a bit different than lift and drag. we keep the 3 components of the moment
        ! in pM in this routine but then the slc%pM variable only has the component of the moment
        ! we are interested in. we use the direction index to get this value out and set it in the slice
        lslc%pM = pm(normal_ind)
        lslc%vM = vm(normal_ind)
 
        ! save the x,y,z-forces into the appropriate real-container
        ! from their type(slice) definition (see the top of this file)
        lslc%fx = pf(1) + vf(1)
        lslc%fy = pf(2) + vf(2)
        lslc%fz = pf(3) + vf(3)
 
        ! Reduce the lift/drag values
        call mpi_reduce((/lslc%pL, lslc%pD, lslc%pM, lslc%vL, lslc%vD, lslc%vM/), tmp, 6, adflow_real, mpi_sum, &
                        0, adflow_comm_world, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Reduce the fx, fy and fz using THE OTHER tmp, i.e. tmp_
        call mpi_reduce((/lslc%fx, lslc%fy, lslc%fz/), tmp_, 3, adflow_real, mpi_sum, &
                        0, adflow_comm_world, ierr)
        call echk(ierr, __file__, __line__)
 
        if (myid == 0) then
            gslc%pL = tmp(1)
            gslc%pD = tmp(2)
            gslc%pM = tmp(3)
            gslc%vL = tmp(4)
            gslc%vD = tmp(5)
            gslc%vM = tmp(6)
 
            ! we now also must remember to store the MPIreduced fx, fy and fz in the
            ! global slice type:
            gslc%fx = tmp_(1)
            gslc%fy = tmp_(2)
            gslc%fz = tmp_(3)
 
            ! Compute factor to get coefficient
            fact = two / (gammainf * pinf * machcoef * machcoef * pref)
 
            ! Take dmax as chord and compute coefficients
            gslc%CLp = gslc%pL / gslc%chord * fact
            gslc%CDp = gslc%pD / gslc%chord * fact
            gslc%CLv = gslc%vL / gslc%chord * fact
            gslc%CDv = gslc%vD / gslc%chord * fact
 
            ! Moment factor has an extra lengthRef
            fact = fact / (lengthref * lref)
 
            ! moments
            gslc%CMp = gslc%pM / gslc%chord * fact
            gslc%CMv = gslc%vM / gslc%chord * fact
 
            ! Default values
            gslc%twist = zero
            gslc%thickness = zero
 
            if (gslc%nNodes == 0) then
                return
            end if
 
            ! Lastly we need the twist and the twist and the thickness
            i1 = bestpair(1)
            i2 = bestpair(2)
            x1 = gslc%vars(1:3, i1)
            x2 = gslc%vars(1:3, i2)
 
            if (x1(1) > x2(1)) then
                te = x1
                le = x2
            else
                te = x2
                le = x1
            end if
 
            ! Save the leading and trailing edges so we can do scaled output
            ! later
            gslc%le = le
            gslc%te = te
 
            ! Finally we need to get the thickness. For this, compute temporary
            ! section nodes and rotate them by the twist values we just computed
            ! and take the max and min
 
            ! Length of hyptoneuse is the same
            hyp = sqrt((x1(1) - x2(1))**2 + (x1(2) - x2(2))**2 + (x1(3) - x2(3))**2)
 
            if (normal_ind == 1) then
                ! Xslice...we don't how what to do here..could be y or z. Don't
                ! do anything.
                gslc%twist = zero
            else if (normal_ind == 2) then
                ! Yslice
                theta = asin((le(3) - te(3)) / hyp)
                gslc%twist = theta * 180.0 / pi
            else
                ! Zslice
                theta = asin((le(2) - te(2)) / hyp)
                gslc%twist = theta * 180.0 / pi
            end if
 
            if (normal_ind == 1) then
                m(1, 1) = one; m(1, 2) = zero; m(1, 3) = zero;
                m(2, 1) = zero; m(2, 2) = one; m(2, 3) = zero;
                m(3, 1) = zero; m(3, 2) = zero; m(3, 3) = one;
            else if (normal_ind == 2) then
                ! Y-rotation matrix
                m(1, 1) = cos(-theta); m(1, 2) = zero; m(1, 3) = sin(-theta);
                m(2, 1) = zero; m(2, 2) = one; m(2, 3) = zero;
                m(3, 1) = -sin(-theta); m(3, 2) = zero; m(3, 3) = cos(-theta);
            else
                ! Z rotation Matrix
                m(1, 1) = cos(theta); m(1, 2) = -sin(theta); m(1, 3) = zero;
                m(2, 1) = sin(theta); m(2, 2) = cos(theta); m(2, 3) = zero;
                m(3, 1) = zero; m(3, 2) = zero; m(3, 3) = one;
            end if
 
            allocate (tempcoords(3, size(gslc%vars, 2)))
            do i = 1, size(gslc%vars, 2)
                ! extract node:
                r = gslc%vars(1:3, i) - te
                r_new = matmul(m, r)
                tempcoords(:, i) = r_new + te
            end do
 
            ! Now get the max and the min and divide by the chord for t/c
            if (normal_ind == 1) then
                gslc%thickness = 0 ! Again, don't know what to do here
            else if (normal_ind == 2) then
                dmax = maxval(tempcoords(3, :))
                dmin = minval(tempcoords(3, :))
                gslc%thickness = (dmax - dmin) / hyp
            else if (normal_ind == 3) then
                dmax = maxval(tempcoords(2, :))
                dmin = minval(tempcoords(2, :))
                gslc%thickness = (dmax - dmin) / hyp
            end if
            deallocate (tempcoords)
        end if
    end subroutine integrateslice
 
    subroutine writeslice(slc, fileID, nFields)
        ! Write the data in slice 'slc' to openfile ID fileID
        use constants
        use inputio
        use commonformats, only: int5
        implicit none
 
        ! Input Parameters
        type(slice), intent(in) :: slc
        integer(kind=intType), intent(in) :: fileid, nfields
 
        ! Working Variables
        integer(kind=intType) :: i, j
        real(kind=realtype) :: tmp, tx, ty, tz
 
        write (fileid, "(a,a,a)") "Zone T= """, trim(slc%sliceName), """"
 
        ! IF we have nodes actually write:
        if (slc%nNodes > 0) then
            write (fileid, *) "Nodes = ", slc%nNodes, " Elements= ", size(slc%conn, 2), " ZONETYPE=FELINESEG"
            write (fileid, *) "DATAPACKING=POINT"
 
            do i = 1, slc%nNodes
                ! Write the coordinates
                do j = 1, 3
                    write (fileid, sci6, advance='no') slc%vars(j, i)
                end do
 
                ! Write the scaled coordiantes with the LE at (0,0,0)
                do j = 1, 3
                    tmp = slc%vars(j, i)
                    write (fileid, sci6, advance='no') (tmp - slc%le(j)) / slc%chord
                end do
 
                ! Write field data. Starts at 9 (after 3 coordindates and the 6 tractions)
                do j = 1, nfields
                    write (fileid, sci6, advance='no') slc%vars(9 + j, i)
                end do
 
                write (fileid, "(1x)")
            end do
 
            do i = 1, size(slc%conn, 2)
                write (fileid, int5) slc%conn(1, i), slc%conn(2, i)
            end do
        else ! Write dummy data so the number of zones are the same
 
            write (fileid, *) "Nodes = ", 2, " Elements= ", 1, " ZONETYPE=FELINESEG"
            write (fileid, *) "DATAPACKING=POINT"
            do i = 1, 2
                do j = 1, 6
                    write (fileid, sci6, advance='no') zero
                end do
 
                do j = 1, nfields
                    write (fileid, sci6, advance='no') zero
                end do
 
                write (fileid, "(1x)")
            end do
            write (fileid, int5) 1, 2
        end if
    end subroutine writeslice
 
end module tecplotio