c-----------------------------------------------------------------------
      subroutine integrate_surface_vs_z(uz,u,nelgz,nelgxy,w)

c     output uz(lz1,nelgz) = surfaced average value of u as a function of z
c     input:  u(lz1,ly1,lz1,nelt)
c             nelgz  = global number of elements in z direction
c             nelgxy = global number of elements in x-y planes
c
c     w is a work array of size = uz

c
      include 'SIZE'
      include 'TOTAL'

      real z(lz1,nelgz),w(lz1,nelgz)
      real u(lx1,ly1,lz1,nelv)

      integer e,ex,ey,ez,exy,eg,f

      n = nelgz*lz1

      do ipass=1,2  ! ipass=1, get weights, ipass=2, get uz


       call rzero(uz,n)

       if (ipass.eq.1) then
          call rone(w,n)       ! Initialize weight to 1
          w1=1
          w2=0
       else
          w1=0
          w2=1
       endif

       do e=1,nelv
       do f=1,4     ! Only need to check faces 1-4 because 5-6 are z-planes

         if (cbc(e,f,1).eq.'W  ') then ! On surface
            eg = lglel(e)
            ez = 1+(eg-1)/nelxy  ! Which z level

            call facind(i0,i1,j0,j1,k0,k1,lx1,ly1,lz1,f)

            l=0
            do k=k0,k1
            do j=j0,j1
            do i=i0,i1
               l=l+1
               scale = w(k,ez)*area(l,1,f,e)
               uz(k,ez) = uz(k,ez) + scale*(w1+w2*u(i,j,k,e))
            enddo
            enddo
            enddo

         endif
       enddo
       enddo

       if (ipass.eq.1) then ! Sum denominator weights across all processors
          call gop     (uz,w,'+  ',n)
          call invcol1 (uz,n)
          call copy    (w,uz,n)  ! Save weights in w
       else
          call gop     (uz,w,'+  ',n)  ! Sum weighted values across processors
       endif
      enddo

      return
      end
c-----------------------------------------------------------------------