imexQ_oc_integrate Subroutine

  subroutine imexQ_oc_integrate(this, lev, qSDC, fSDC, dt, fintSDC, flags)
    class(pf_imexQ_oc_t), intent(inout) :: this
    class(pf_level_t), intent(in   ) :: lev
    class(pf_encap_t), intent(in   ) :: qSDC(:), fSDC(:, :) !qSDC unused?
    real(pfdp),       intent(in)     :: dt
    class(pf_encap_t), intent(inout) :: fintSDC(:)
    integer, intent(in), optional    :: flags

    integer :: n, m, Nnodes, which
    Nnodes=lev%nnodes
    
    which = 0
    if(present(flags)) then 
      which = flags
    else
      print *, "flags not present in integrate", which
      stop
    end if
!     print *, "IMEXQ_OC INTEGRATE ", which
        
    do n = 1, Nnodes-1
       !  Forward in y
       if( (which .eq. 0) .or. (which .eq. 1) ) then
          call fintSDC(n)%setval(0.0_pfdp, 1)
          do m = 1, Nnodes
!               do p = 1, npieces
!                 call Lev%encap%axpy(fintSDC(n), dt*Lev%sdcmats%qmat(n,m), fSDC(m,p),1)
!               end do
            if (this%explicit) &
               call fintSDC(n)%axpy(dt*lev%sdcmats%qmat(n,m), fSDC(m,1), 1)
            if (this%implicit) &
               call fintSDC(n)%axpy(dt*lev%sdcmats%qmat(n,m), fSDC(m,2), 1)
          end do
       end if
       
       !  Backward in p
       if( (which .eq. 0) .or. (which .eq. 2) ) then
          call fintSDC(Nnodes-n)%setval(0.0_pfdp, 2)
          do m = 1, Nnodes
!               do p = 1, npieces
!                 call Lev%encap%axpy(fintSDC(Nnodes-n), dt*Lev%sdcmats%qmat(n,m), fSDC(Nnodes+1-m,p),2)
!               end do
            if (this%explicit) &
               call fintSDC(Nnodes-n)%axpy(dt*lev%sdcmats%qmat(n,m), fSDC(Nnodes+1-m,1), 2)
            if (this%implicit) &
               call fintSDC(Nnodes-n)%axpy(dt*lev%sdcmats%qmat(n,m), fSDC(Nnodes+1-m,2), 2)
          end do
       end if
    end do
  end subroutine imexQ_oc_integrate