ColumnIntegral_ColumnOperators3D_t Subroutine

  subroutine ColumnIntegral_ColumnOperators3D_t(this,q,qbar)
    !! Column integral of a volume field,
    !!
    !!   qbar_ij = sum_{k=1..nz} sum_{l=0..N} w_l (Z_zeta)^(k)_ijl q^(k)_ijl ,
    !!
    !! for every column and variable. With Z_zeta in metres, qbar has the unit of q
    !! times metres; the column integral of 1 is the local depth. On Gauss-Lobatto nodes
    !! the quadrature is exact when Z_zeta q is a polynomial of degree 2N-1 in zeta in
    !! every layer, on Gauss nodes when it is of degree 2N+1.
    !!
    !!  Input
    !!    - q : volume field on the extruded mesh
    !!  Output
    !!    - qbar : surface field, overwritten
    implicit none
    class(ColumnOperators3D_t),intent(in) :: this
    class(Scalar3D),intent(in) :: q
    class(Scalar2D),intent(inout) :: qbar
    ! Local
    integer :: i,j,l,k,e2,iel,ivar,nz

    call CheckFields_ColumnOperators3D_t(this,q,qbar,"ColumnIntegral")
    nz = this%nLayers

    ! The sum runs over the layers from the bottom up and, inside a layer, over the
    ! nodes from the bottom up, the same order as the device kernel. The innermost loop
    ! is the contiguous index i.
    do ivar = 1,q%nVar
      do e2 = 1,this%nColumns
        do j = 1,this%N+1
          do i = 1,this%N+1
            qbar%interior(i,j,e2,ivar) = 0.0_prec
          enddo
        enddo
        do k = 1,nz
          iel = (e2-1)*nz+k
          do l = 1,this%N+1
            do j = 1,this%N+1
              do i = 1,this%N+1
                qbar%interior(i,j,e2,ivar) = qbar%interior(i,j,e2,ivar)+ &
                                             this%interp%qWeights(l)* &
                                             this%geometry%dxds%interior(i,j,l,iel,1,3,3)* &
                                             q%interior(i,j,l,iel,ivar)
              enddo
            enddo
          enddo
        enddo
      enddo
    enddo

  endsubroutine ColumnIntegral_ColumnOperators3D_t