VerticalSweepUp_ColumnOperators3D_t Subroutine

  subroutine VerticalSweepUp_ColumnOperators3D_t(this,r,qBottom,q)
    !! Upward sweep: solves, layer by layer from the bottom of every column,
    !!
    !!   (D + M^{-1} e_0 e_0^T) q^(k) = r^(k) + M^{-1} e_0 q_in^(k) ,
    !!   q_in^(1) = qBottom ,  q_in^(k) = q^(k-1)_N ,
    !!
    !! on every vertical line (i,j) and for every variable, as q^(k) = q_in^(k) 1 + A^{-1} r^(k).
    !! This is the upwind DG discretization of dq/dzeta = r with the inflow at the bottom
    !! of each element. r is in reference scaling: its unit is that of q, per unit zeta,
    !! and no metric factor is applied (see the module header). Requires Gauss-Lobatto
    !! control nodes.
    !!
    !!  Input
    !!    - r : right hand side, volume field
    !!    - qBottom : inflow at the bottom of every column, surface field
    !!  Output
    !!    - q : solution, volume field, overwritten. It must be a different object from r.
    implicit none
    class(ColumnOperators3D_t),intent(in) :: this
    class(Scalar3D),intent(in) :: r
    class(Scalar2D),intent(in) :: qBottom
    class(Scalar3D),intent(inout) :: q
    ! Local
    integer :: i,j,l,m,k,e2,iel,ivar,nz,np

    call CheckSweep_ColumnOperators3D_t(this,r,qBottom,q,"VerticalSweepUp")
    nz = this%nLayers
    np = this%N+1

    ! For each node l the particular solution sum_m upInverse(m,l) r_m is accumulated in
    ! q in the order m = 1..N+1, and the inflow is added last, the same operations as
    ! the device kernel. The innermost loop is the contiguous index i.
    do ivar = 1,r%nVar
      do e2 = 1,this%nColumns
        do k = 1,nz
          iel = (e2-1)*nz+k
          do l = 1,np
            do j = 1,np
              do i = 1,np
                q%interior(i,j,l,iel,ivar) = 0.0_prec
              enddo
            enddo
            do m = 1,np
              do j = 1,np
                do i = 1,np
                  q%interior(i,j,l,iel,ivar) = q%interior(i,j,l,iel,ivar)+ &
                                               this%upInverse(m,l)*r%interior(i,j,m,iel,ivar)
                enddo
              enddo
            enddo
          enddo
          ! Inflow: the bottom trace for the bottom layer, else the top node of the layer below
          if(k == 1) then
            do l = 1,np
              do j = 1,np
                do i = 1,np
                  q%interior(i,j,l,iel,ivar) = qBottom%interior(i,j,e2,ivar)+q%interior(i,j,l,iel,ivar)
                enddo
              enddo
            enddo
          else
            do l = 1,np
              do j = 1,np
                do i = 1,np
                  q%interior(i,j,l,iel,ivar) = q%interior(i,j,np,iel-1,ivar)+q%interior(i,j,l,iel,ivar)
                enddo
              enddo
            enddo
          endif
        enddo
      enddo
    enddo

  endsubroutine VerticalSweepUp_ColumnOperators3D_t