module terms_mod
  use main_data
  use define_state
  use dwr_mod
  use euler_problem
  use lin_solvers
  use inviscid_fluxes
  use stdgm_mod
  use matrix_oper_int
  use mesh_mod
  use eval_jumps
  use modelPorous
  use higher_order_local
  use pm_fluxes

  implicit none

  public :: ComputeST_Terms
  public :: ComputeST_Terms_new
  public :: computeDualResidualPlus
  public :: ComputeSTDiag_Terms_loc

  !  public :: FillBigVector
  !SetF_q_Cw_fast


contains

 !> evaluation of the new big matrix and vector blocks for STDGM method
 !> if implicitly == False -> vector F(u) saved in elem%vec(rhs,:)
 !> if implicitly == True  -> sets the matrix C(*,*) elem%vec(rhs = q(u),
  subroutine ComputeST_Terms( deg_plus )
    logical, intent (in) :: deg_plus
    class(element), pointer :: elem
    class(Time_rule), pointer :: T_rule
    real, dimension(:,:), allocatable :: T_mat
    real, dimension(:), allocatable :: T_vec
    integer :: i, j, k, l, r, kk, rr, ll
    integer :: alpha, Qdeg, dof, Tdof, s_dim, f_dof, f_Tdof, f_dim

    integer :: m, mm, n, nn, wdof, wTdof
    real :: cTime, tt, tt1
    real :: val, val_vec, val_vecD, val0, val1, val2, val3, val_sum
    real :: local_eta, tval
    logical ::iprint, CD_tdt
    integer :: itest
    integer :: nd1, nd2, d1, d2, p_mod_max, q_mod_max

    !print*,' !!!!!!! starting ComputeST_Terms', deg_plus , state%nlSolver%implicitly


    !if(state%nlSolver%implicitly) then
    !   print*,''
    !   print*,'ComputeST_Terms( deg_plus )', deg_plus, state%nlSolver%implicitly
    !   print*,''
    !   print*,''
    !endif
!    print*,'ComputeST_Terms( deg_plus )', deg_plus, "impl = ", state%nlSolver%implicitly
!    print*,'wActual = ', state%wActual

    if ( deg_plus ) then
      p_mod_max = state%p_mod_max
      q_mod_max = state%q_mod_max

      grid%elem(:)%deg_plus = .true.

    else
      p_mod_max = 0
      q_mod_max = 0
      grid%elem(:)%deg_plus = .false.
    end if

    !print*, "degplus in calling UpdateMatrixElementShapeForBigBlock", deg_plus

    !CD_tdt = .false.  ! variant CD_tdt - : normal form with linearization and tjumps
    CD_tdt = .true.   ! variant CD_tdt +  : abnormal form without linearization and tjumps

    if(state%nlSolver%implicitly) then
      state%num_call_C = state%num_call_C + 1
    else
      state%num_call_F = state%num_call_F + 1
    endif

    !print*, 'ComputeST_Terms with implicitly = ' ,state%nlSolver%implicitly, 'deg_plus:', deg_plus
    call cpu_time(tt1)

    cTime = state%time%ctime

    local_eta = 1. / state%time%tau(1)

    ! adding of the term in front of the time derivative
!    if(state%model%varying_time_term) then
      s_dim =  maxval( grid%elem(:)%dof_plus)*ndim
      allocate( T_mat(1: s_dim, 1:s_dim) )
      allocate( T_vec(1: s_dim) )
!    endif

    !NEW
    ! set degree with respect to time
    Qdeg = state%time%Qnum

    ! TODO FR - Qdeg should be connected with q_mod_max !
    T_rule => state%time%T_rule(Qdeg)

    associate ( time => state%time)
    select type ( time )
      class is ( TimeTDG_t )
      class default
        stop 'computeST_Terms only for STDG'
    end select
    end associate

!!! IMPLICITLY = TRUE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

         ! rhsST = q(w)
    if(state%nlSolver%implicitly) then
         !print*, 'ComputeST_Terms called with implicitly true, deg_plus = ', deg_plus

         call ClearMatrixBigBlocks()
         ! TODO - length of STDGM vectors
         call ClearVectorBlocks_STDGM()    ! VD Newton - nenuluje se zbytecne?

         state%linSolver%precond_update = .true.

         do alpha = 1, Qdeg ! temporarily max_Tdof =  max time quadrature nodes

            do i = 1, grid%nelem
              elem => grid%elem(i)

                ! save the wST space-time solution in quadrature index alpha to w
                call Transfer_wST_to_wActual_Elem(elem , alpha, Qdeg)
                ! save the wST space-time solution in quadrature index alpha to w
                call Transfer_wST_to_w_Elem(elem , alpha, Qdeg)

            enddo


         !we have to run ComputeTerms() in the time-quadrature nodes
            !cTime = state%time%ctime
            state%time%ctime = state%time%ttime + state%time%tau(1) * T_rule%lambda(alpha)

            call ComputeTerms(deg_plus )

            !print*, 'First entry of elem1%block =' , size(grid%elem(1)%block(0)%Mb)
            state%time%ctime = cTime


            do i =1, grid%nelem
              elem => grid%elem(i)

              Tdof = elem%Tdof + q_mod_max
              dof = DOFtriang( elem%deg + p_mod_max)
              s_dim = ndim *dof

!              print*, "bigBlock = ", i;
!              print*, "blocck", size(elem%block(0)%Mb,1), size(elem%block(0)%Mb,2), dof
!              print*, "bigBlockss", elem%bigBlock(0)%dimensions(), dof

              ! if(elem%i ==80) then
              !    write(31,*)'  _______ DV__________ ', alpha
              !    do j=1,elem%dof
              !       write(31,'(a10, 2i5, 200es14.6)') 'Mbl w1:', elem%i, j, &
              !            elem%bigBlock(3)%Mb(j,:,1,1,1,1)
              !    enddo
              !    write(31,*)
              !    do j=1,dof
              !       write(31,'(a10, i5, 200es14.6)') 'Mblock:', j,elem%block(3)%Mb(j,:)
              !    enddo
              !    write(31,*)'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
              ! endif

              do l = 1, Tdof
                 !diag blocks of bigBlock
                 do r = 1, Tdof
                    tval = T_rule%phi(l,alpha)*T_rule%phi(r,alpha) * T_rule%weights(alpha)
                    nd1 = 1
                    do d1 = 1,ndim
                       nd2 = 1
                       do d2 = 1, ndim
                          !print*, "cdsace" , dof, d1, d2, l, r, elem%bigBlock(0)%Mb( 1, 1, d1, d2, l, r )
                          !print*, "cdaaw" , nd1,nd1+dof-1,nd2,nd2+dof-1
                          !print*, "block -- ", elem%block(0)%Mb(nd1:nd1+dof-1,nd2:nd2+dof-1)
                          
                          elem%bigBlock(0)%Mb( 1:dof, 1:dof, d1, d2, l, r ) = &
                               elem%bigBlock(0)%Mb( 1:dof, 1:dof, d1, d2, l, r )  + &
                               elem%block(0)%Mb(nd1:nd1+dof-1,nd2:nd2+dof-1) * tval
                          nd2 = nd2 + dof
                       end do
                       nd1 = nd1 + dof
                    end do ! d1
                    
                 enddo !r
                 
                 !offdiagonal blocks of bigbLock
                 do j = 1, elem%flen
                   if(elem%face(neigh,j) > 0) then
                      !f_dof = DOFtriang( elem%face(fdeg,j) + p_mod_max )
                      !FR - alocBig, VD Newton
                      if (deg_plus) then
                        f_dof = elem%face(fdof_plus,j)
                      else
                        f_dof = elem%face(fdof,j)
                      endif


                      f_Tdof = elem%face(fTdof,j) + q_mod_max
                      f_dim = ndim * f_dof


                      do r = 1, f_Tdof
                         tval = T_rule%phi(l,alpha)*T_rule%phi(r,alpha) * T_rule%weights(alpha)
                         nd1 = 1
                         do d1 = 1,ndim
                            nd2 = 1
                            do d2 = 1, ndim
                               !write(*,'(a8, 10i5)') 'hd37u2w;', elem%i, j, elem%deg, &
                               !     nd1, dof, nd2, f_dof, &
                               !     size(elem%block(j)%Mb, 1), size(elem%block(j)%Mb, 2)

!                               write(*,'(a8, 10i5)') 'hd37u2w;', elem%i, j, elem%deg, &
!                                    nd1, dof, nd2, f_dof, &
!                                    size(elem%block(j)%Mb, 1), size(elem%block(j)%Mb, 2)
!                               print*, "---"
!                               print*, "size block = ", size(elem%block(j)%Mb, 1), size(elem%block(j)%Mb, 2)
!                               print*, "size big block = ", elem%bigBlock(j)%dimensions()
!                               print*, "sizes === ", nd1,nd2 !+f_dof-1

                               elem%bigBlock(j)%Mb( 1:dof, 1:f_dof, d1, d2, l, r ) = &
                                    elem%bigBlock(j)%Mb( 1:dof, 1:f_dof, d1, d2, l, r ) + &
                                    elem%block(j)%Mb(nd1:nd1+dof-1,nd2:nd2+f_dof-1) * tval
                               nd2 = nd2 + f_dof
                            end do
                            !print*, "after inner cycle"
                            nd1 = nd1 + dof
                         end do ! d1

                      enddo !r=1,f_Tdof
                   endif
                enddo !j
             enddo !l = 1,Tdof

             ! TODO FR_ST - needs to be tested how varying_time_term works???
             !if (i==1) print*, "FR_ST - needs to be tested how varying_time_term works???"
             ! adding of the term in front of the time derivative - time deriv term itself
             ! for (state%modelName /= 'porous' ) it should give simply the Mass matrix
             if (state%model%varying_time_term) then
               call Compute_Time_deriv(elem, .false., dof*ndim,  T_mat(1:dof*ndim,1:dof*ndim) )

               do l=1,Tdof
                  do r = 1, Tdof
                     ! s_dim = dof*ndim
                     ll = (l -1) * s_dim
                     rr = (r -1) * s_dim

                     tval = T_rule%phi(l,alpha)*T_rule%Dphi(r,alpha) &
                          * T_rule%weights(alpha) * local_eta

                     nd1 = 1
                     do d1 =1,ndim
                        nd2 = 1
                        do d2 = 1,ndim
                           elem%bigBlock(0)%Mb( 1:dof, 1:dof, d1, d2, l, r ) = &
                                elem%bigBlock(0)%Mb( 1:dof, 1:dof, d1, d2, l, r ) &
                                + tval * T_mat(nd1:nd1+dof-1,nd2:nd2+dof-1)
                           nd2 = nd2 + dof
                        end do
                        nd1 = nd1 + dof
                     end do
                  enddo  ! r
               enddo ! l

             end if

!               print*, 'norm of bigBlock 0 =' , norm2(elem%bigBlock(0)%Mb)
!               print*, 'bigBlock  6,6 =' , elem%bigBlock(0)%Mb(6,6,1,1,1,1)
            enddo !i

         end do !alpha =  1,Qdeg


         if(state%model%precomputed_arrays ) then  ! we have to precompute the  values
          stop "precomputed_arrays are not done in ComputeST_Terms!"
         endif


         if(state%model%varying_time_term) then
            do i = 1, grid%nelem
               elem => grid%elem(i)

               Tdof = elem%getActualTDof()
               dof = elem%getActualDof()
               s_dim = ndim *dof

               ! save the wST space-time solution in quadrature index alpha to w
               call Transfer_wST_to_w_Elem(elem , -1, Qdeg)
               call Transfer_wST_to_wActual_Elem(elem , -1, Qdeg)

               ! computing of the matrix in front of the time derivative term at t_{m-1}^+
               call Compute_Time_deriv(elem, .false., s_dim, T_mat(1:s_dim,1:s_dim) )

               do l=1,Tdof
                  do r = 1, Tdof
                     ll = (l -1) * s_dim
                     rr = (r -1) * s_dim

                     tval = T_rule%phi(l,-1)*T_rule%phi(r,-1) * local_eta
                     ! TODO - FR_ST where we use loc_eta: comparison between (w',v) and the jump term

                     nd1 = 1
                     do d1 =1,ndim
                        nd2 = 1
                        do d2 = 1,ndim
                           elem%bigBlock(0)%Mb( 1:dof, 1:dof, d1, d2, l, r ) = &
                                elem%bigBlock(0)%Mb( 1:dof, 1:dof, d1, d2, l, r ) &
                                + tval * T_mat(nd1:nd1+dof-1,nd2:nd2+dof-1)
                           nd2 = nd2 + dof
                        end do
                        nd1 = nd1 + dof
                     end do
                  enddo  ! r
               enddo ! l

            enddo ! i=1,grid%nelem

         end if





!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! IMPLICITLY = FALSE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      else   ! implicitly == .false.

        state%cons = 0. ! auxiliary array for verifying the code

        iprint = .false.

        call ClearVectorBlocks_STDGM()  ! rhsST = 0


        cTime = state%time%ctime  ! storing of the actual value

        itest = -13   ! 9

        do alpha = 1, Qdeg ! temporarily max_Tdof =  max time quadrature nodes

            ! verifying the consistency of the scheme
            state%cons(0,1) = T_rule%weights(alpha)

            do i = 1, grid%nelem
               elem => grid%elem(i)
               call Transfer_wST_to_w_Elem(elem , alpha, Qdeg)
               call Transfer_wST_to_wActual_Elem(elem , alpha, Qdeg)
            enddo

            !we have to run ComputeTerms() in the time-quadrature nodes
            state%time%ctime = state%time%ttime + state%time%tau(1) * state%time%T_rule(Qdeg)%lambda(alpha)

            call ComputeTerms( deg_plus )


            do i =1, grid%nelem
               elem => grid%elem(i)

               !write(*,'(a10, i5, 300es12.4)') '!! ELM:', i, elem%vec(rhs, :)

               !NEW for adaptation
               wdof = elem%dof
               wTdof = elem%Tdof

               !NEW for adaptation
               if(elem%deg_plus) then
                  dof = elem%dof_plus
                  Tdof = elem%Tdof_plus
               else
                  dof = elem%dof
                  Tdof = elem%Tdof
               endif

               s_dim = ndim * dof


               !vec -> rhsST
               do k = 1,ndim
                  do l = 1, Tdof
                     val_vec = T_rule%phi(l,alpha) * T_rule%weights(alpha)

                     elem%rhsST(k, 1:dof, l) = elem%rhsST(k, 1:dof, l)  &
                          +  val_vec *  elem%vec( rhs, (k-1)*dof + 1 : k*dof ) 	! dofA

                     if(k==1 .and. l == 1) state%cons(1,1) = state%cons(1,1)&
                          + val_vec *  elem%vec( rhs, (k-1)*dof + 1 )

                  enddo !l
               enddo !k


               ! adding of the term in front of the time derivative - time deriv term itself
               if(state%model%varying_time_term) then

                  ! computing of the matrix in front of the time derivative term at t_alpha
                  !print*,'Compute_Time_deriv C', s_dim, elem%dof_plus,dofA
                  call Compute_Time_deriv(elem, elem%deg_plus, s_dim,  T_mat(1:s_dim,1:s_dim) )
                  !compute the vector in the time derivative
                  call Compute_Time_deriv_Vector(elem, dof,  T_vec(1: s_dim) )


                  do k = 1,ndim
                     do l = 1, Tdof

                        val_vec = T_rule%phi(l,alpha) * T_rule%weights(alpha) * local_eta

                        val_vecD = T_rule%Dphi(l,alpha) * T_rule%weights(alpha) * local_eta

                        do j=1, dof

                           if(.not. CD_tdt) then
                              ! CD_tdt -
                              elem%rhsST(k, j, l) = elem%rhsST(k, j, l)  &
                                   -  val_vec  &
                                   * dot_product(T_rule%Dphi(1:wTdof,alpha) , &
                                   matmul( T_mat( (k-1)* dof + j, (k-1)*dof+1:(k-1)*dof + wdof ),&
                                   elem%wST(k, 1:wdof, 1:wTdof) ) )
                              ! CD_tdt -
                           else
                              ! CD_tdt +
                              !write(*,'(a8,3i5, 30es12.4)') 'T_v$RE',elem%i, dof, (k-1)* dof +j, T_vec(1: s_dim)

                              elem%rhsST(k, j, l) = elem%rhsST(k, j, l) + val_vecD * T_vec( (k-1)* dof +j)
                              ! CD_tdt +
                           endif


                        enddo
                     enddo  ! l=1,Tdof
                  enddo ! k=1,ndim
               endif ! state%model%varying_time_term


            enddo !i

         end do !alpha =  1,Qdeg


         ! adding the time derivative and the the w_{m-1}^+ part:
         if(state%model%precomputed_arrays ) then  ! we have to precompute the  values
            do i = 1, grid%nelem
               elem => grid%elem(i)

               ! save the wST space-time solution in quadrature index alpha to w
               call Transfer_wST_to_w_Elem(elem , -1, Qdeg)
               call Transfer_wST_to_wActual_Elem(elem , -1, Qdeg)
            enddo
            !if(state%nlSolver%iter <= 5) &
            call ComputeCapacityConductivity(.true.) ! DO NOT comment, used in inv_fluxes !!!!!!!!!!
            !write(*,'(a7,200es12.4)')'CAP Y:',grid%elem(741)%xi(0, 1:3, iRe+4),grid%elem(741)%xi(0, 1:3, iRe+3)

         endif

         !val1 = 0.

         do i =1, grid%nelem
            elem => grid%elem(i)

            !NEW for adaptation
            wdof = elem%dof
            wTdof = elem%Tdof

            !NEW for adaptation
            if(elem%deg_plus) then
               dof = elem%dof_plus
               Tdof = elem%Tdof_plus
            else
               dof = elem%dof
               Tdof = elem%Tdof
            endif

            ! original subroutine for the adding the time derivative and the the w_{m-1}^+ part
            if(.not.  state%model%varying_time_term ) then
               do k = 1, ndim
                  do l = 1,Tdof
                     do j = 1,dof   ! dofA
                        !NEW for adaptation
                        val_vec = local_eta * dot_product( state%time%refTimeMatrix%Mb(l,1:wTdof) , &
                             matmul( elem%Mass%Mb(j, 1:wdof), elem%wST(k,1:wdof,1:wTdof) ))

                        elem%rhsST(k, j, l) = elem%rhsST(k, j, l) - val_vec

                     enddo !j
                  enddo !l

               enddo !k
            endif


            ! adding of the term in front of the time derivative -  jump term,  w_{m-1}^+ part:
            if(state%model%varying_time_term) then

               call Transfer_wST_to_w_Elem(elem , -1, Qdeg)
               call Transfer_wST_to_wActual_Elem(elem , -1, Qdeg)

               ! computing of the matrix in front of the time derivative term at t_{m-1}^+
               s_dim = ndim * dof

               !compute the vector in the time derivative
               call Compute_Time_deriv_Vector(elem, dof,  T_vec(1: s_dim) )


               do k = 1,ndim
                  do l = 1, Tdof

                     val_vec = T_rule%phi(l,-1) * local_eta

                     do j=1, dof
                        ! formula without linearization
                        if(.not. CD_tdt) then
                           ! CD_tdt -
                           elem%rhsST(k, j, l) = elem%rhsST(k, j, l) - val_vec * T_vec( (k-1)*dof +j )
                        endif
                     enddo
                  enddo  ! l=1,Tdof
               enddo ! k=1,ndim

               ! !! VARIANT CD_tdt +
               call Transfer_wST_to_w_Elem(elem , 0, Qdeg)
               call Transfer_wST_to_wActual_Elem(elem , 0, Qdeg)

               ! computing of the matrix in front of the time derivative term at t_{m-1}^+
               s_dim = ndim * dof
               !print*,'Compute_Time_deriv D', s_dim, elem%dof_plus
               call Compute_Time_deriv(elem, elem%deg_plus, s_dim, T_mat(1:s_dim, 1:s_dim) )

               !compute the vector in the time derivative
               call Compute_Time_deriv_Vector(elem, dof,  T_vec(1: s_dim) )

               do k = 1, ndim
                  do l = 1, Tdof

                     val_vec = T_rule%phi(l, 0) * local_eta

                     do j=1, dof

                        ! CD_tdt -  .. is empty
                        if(CD_tdt) then
                           ! CD_tdt +
                           elem%rhsST(k, j, l) = elem%rhsST(k, j, l) - val_vec * T_vec( (k-1)*dof +j )
                        endif
                     enddo
                  enddo  ! l=1,Tdof
               enddo ! k=1,ndim
               ! !! END OF THE TEST


            endif   !state%model%varying_time_term
         enddo !i

    !         write(*,'(a15, 12es12.4)') 'F^m: aft mass' , grid%elem(1)%rhsST(1,:,:)

         ! the w_{m-1}^- part:   "state%model%varying_time_term" given inside

         if(state%model%precomputed_arrays ) then  ! we have to precompute the  values
            do i = 1, grid%nelem
               elem => grid%elem(i)

               ! save the wST space-time solution in quadrature index alpha to w
               call Transfer_wST_to_w_Elem(elem , -1, Qdeg)
               call Transfer_wST_to_wActual_Elem(elem , -1, Qdeg)
            enddo
            !if(state%nlSolver%iter <= 5) &
            call ComputeCapacityConductivity(.true.) ! DO NOT comment, used in inv_fluxes !!!!!!!
            !write(*,'(a7,200es12.4)')'CAP Z:',grid%elem(741)%xi(0, 1:3, iRe+4),grid%elem(741)%xi(0, 1:3, iRe+3)

         endif


         !val = 0.
         do i = 1, grid%nelem
            elem => grid%elem(i)

            ! state%cons(4,1) filled inside
            call Elem_wSTfinToRhsST( grid%elem(i) , state%nlSolver%implicitly)

            ! if(i == 132) then
            ! !if(i == 249) then
            !    print*,'dof = ',dof, 'Tdof =', Tdof

            !    do k= 1, elem%Tdof
            !       !do j=1, dof
            !          write(*,'(a4, 3i5, 300es12.4)') 'wST<',Tdof, i,k, elem%wST(1, :, k) 
            !       !enddo
            !    enddo
            !    do k= 1, Tdof
            !       !do j=1, dof
            !          write(*,'(a4, 3i5, 300es12.4)') 'dex<',Tdof, i,k, elem%rhsST(1, :, k) 
            !       !enddo
            !    enddo
            ! endif
            !write(*,'(a6,3i5,30es12.4)') "T_vec%ED#", state%nlSolver%iter, &
            !           elem%i, dof,  elem%rhsST(:, :, :)


         enddo !i


      endif   !state%nlSolver%implicitly



    ! putting back the original value
    state%time%ctime = cTime

    if(state%model%varying_time_term) deallocate( T_mat, T_vec)

    !deallocate(cons_test)

    call cpu_time(tt)

    ! print*,'###########################',state%nlSolver%implicitly, state%nlSolver%iter
    ! do i=1, grid%nelem
    !    elem => grid%elem(i)
    !    do k=1, ndim
    !       write(*,'(a10, 2i5, 2es12.4, a3,200es12.4)') &
    !            'STterms w:', elem%i, k, elem%xc(:), '|', &
    !            elem%rhsST(k,1:elem%dof, 1:elem%Tdof)
    !    enddo
        
    ! enddo
    

  end subroutine ComputeST_Terms

  !> evaluation of the new big matrix and vector blocks for STDGM method
  !> new version needed for the DWR_P method for computations wST_plus
 !> if implicitly == False -> vector F(u) saved in elem%vec(rhs,:)
 !> if implicitly == True  -> sets the matrix C(*,*) elem%vec(rhs = q(u),
 !> new version:
 !> if wSTplus == TRUE -> wSTplus -> elem%w(0,:)
 !> else wST -> elem%w(0,:)
  subroutine ComputeST_Terms_new( deg_plus, wSTplus )
    logical, intent (in) :: deg_plus
    logical, intent (in) :: wSTplus
    class(element), pointer :: elem
    class(Time_rule), pointer :: T_rule
    real, dimension(:,:), allocatable :: T_mat
    real, dimension(:), allocatable :: T_vec
    integer :: i, j, k, l, r, kk, rr, ll
    integer :: alpha, Qdeg, dof, Tdof, s_dim, f_dof, f_Tdof, f_dim

    integer :: m, mm, n, nn, wdof, wTdof
    real :: cTime, tt, tt1
    real :: val, val_vec, val_vecD, val0, val1, val2, val3, val_sum
    real :: local_eta, tval
    logical ::iprint, CD_tdt
    integer :: itest
    integer :: nd1, nd2, d1, d2 !, p_mod_max, q_mod_max
    integer :: p_mod_loc, q_mod_loc

    p_mod_loc = state%getP_mod()
    q_mod_loc = state%getQ_mod()

    if (.not. deg_plus .and. wSTplus) &
      stop "ComputeST_Terms_new: it make no sense to have wSTplus with deg_plus=FALSE!"

    if (.not. deg_plus .and. state%getP_mod() /= 0) &
      stop "deg_plus=FALSE && p_mod > 0 in ComputeST_Terms_new"

    if ( deg_plus ) then
      !p_mod_max = state%p_mod_max
      !q_mod_max = state%q_mod_max

      call state%setP_mod( state%p_mod_max )
      call state%setQ_mod( state%q_mod_max )
      grid%elem(:)%deg_plus = .true.
    else
      !p_mod_max = 0
      !q_mod_max = 0
      grid%elem(:)%deg_plus = .false.
    end if

    !CD_tdt = .false.  ! variant CD_tdt - : normal form with linearization and tjumps
    CD_tdt = .true.   ! variant CD_tdt +  : abnormal form without linearization and tjumps

    if(state%nlSolver%implicitly) then
      state%num_call_C = state%num_call_C + 1
    else
      state%num_call_F = state%num_call_F + 1
    endif

    !print*, 'ComputeST_Terms with implicitly = ' ,state%nlSolver%implicitly, 'deg_plus:', deg_plus
    call cpu_time(tt1)

    cTime = state%time%ctime

    local_eta = 1. / state%time%tau(1)

    ! adding of the term in front of the time derivative
    if(state%model%varying_time_term) then
      s_dim =  maxval( grid%elem(:)%dof_plus)*ndim
      allocate( T_mat(1: s_dim, 1:s_dim) )
      allocate( T_vec(1: s_dim) )
    endif

    !NEW
    ! set degree with respect to time
    Qdeg = state%time%Qnum

    ! TODO FR - Qdeg should be connected with q_mod_max !
    T_rule => state%time%T_rule(Qdeg)

    associate ( time => state%time)
    select type ( time )
      class is ( TimeTDG_t )
      class default
        stop 'computeST_Terms only for STDG'
    end select
    end associate

!!! IMPLICITLY = TRUE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

         ! rhsST = q(w)
      if(state%nlSolver%implicitly) then
         !print*, 'ComputeST_Terms called with implicitly true, deg_plus = ', deg_plus

         call ClearMatrixBigBlocks()
         ! TODO - length of STDGM vectors
         call ClearVectorBlocks_STDGM()    ! VD Newton - nenuluje se zbytecne?

         state%linSolver%precond_update = .true.

         do alpha = 1, Qdeg ! temporarily max_Tdof =  max time quadrature nodes

            do i = 1, grid%nelem
              elem => grid%elem(i)

               ! save the wST space-time solution in quadrature index alpha to w
              if (wSTplus) then
                call Transfer_wSTplus_to_wActual_Elem(elem , alpha, Qdeg)
              else
                call Transfer_wST_to_wActual_Elem(elem , alpha, Qdeg)
              end if

            enddo

            !we have to run ComputeTerms() in the time-quadrature nodes
            !cTime = state%time%ctime
            state%time%ctime = state%time%ttime + state%time%tau(1) * T_rule%lambda(alpha)

            call ComputeTerms(deg_plus )

            state%time%ctime = cTime

            do i =1, grid%nelem
              elem => grid%elem(i)
              Tdof = elem%getActualTDof()
              dof = elem%getActualDof()
              s_dim = ndim *dof

              do l = 1, Tdof
                  !diag blocks of bigBlock
                do r = 1, Tdof
                  tval = T_rule%phi(l,alpha)*T_rule%phi(r,alpha) * T_rule%weights(alpha)
                  nd1 = 1
                  do d1 = 1,ndim
                    nd2 = 1
                    do d2 = 1, ndim
                      elem%bigBlock(0)%Mb( 1:dof, 1:dof, d1, d2, l, r ) = &
                        elem%bigBlock(0)%Mb( 1:dof, 1:dof, d1, d2, l, r )  + &
                        elem%block(0)%Mb(nd1:nd1+dof-1,nd2:nd2+dof-1) * tval
                      nd2 = nd2 + dof
                    end do
                    nd1 = nd1 + dof
                  end do ! d1

                enddo !r

                !offdiagonal blocks of bigBlock
                do j = 1, elem%flen
                   if(elem%face(neigh,j) > 0) then
                      !f_dof = DOFtriang( elem%face(fdeg,j) + p_mod_max )
                      !FR - alocBig, VD Newton
                      if (deg_plus) then
                        f_dof = elem%face(fdof_plus,j)
                      else
                        f_dof = elem%face(fdof,j)
                      endif
                                            
                      f_dof = grid%elem(elem%face(neigh,j))%getActualDof()
                      f_Tdof = grid%elem(elem%face(neigh,j))%getActualTDof()
                      f_dim = ndim * f_dof

                      do r = 1, f_Tdof
                         tval = T_rule%phi(l,alpha)*T_rule%phi(r,alpha) * T_rule%weights(alpha)
                         nd1 = 1
                         do d1 = 1,ndim
                            nd2 = 1
                            do d2 = 1, ndim
                              elem%bigBlock(j)%Mb( 1:dof, 1:f_dof, d1, d2, l, r ) = &
                                    elem%bigBlock(j)%Mb( 1:dof, 1:f_dof, d1, d2, l, r ) + &
                                    elem%block(j)%Mb(nd1:nd1+dof-1,nd2:nd2+f_dof-1) * tval
                               nd2 = nd2 + f_dof
                            end do
                            !print*, "after inner cycle"
                            nd1 = nd1 + dof
                         end do ! d1

                      enddo !r=1,f_Tdof
                   endif
                enddo !j
             enddo !l = 1,Tdof


              ! TODO FR_ST - needs to be tested how varying_time_term works???
             !if (i==1) print*, "FR_ST - needs to be tested how varying_time_term works???"
             ! adding of the term in front of the time derivative - time deriv term itself
             ! for (state%modelName /= 'porous' ) it should give simply the Mass matrix
             if (state%model%varying_time_term) then
               call Compute_Time_deriv(elem, .false., dof*ndim,  T_mat(1:dof*ndim,1:dof*ndim) )

               do l=1,Tdof
                do r = 1, Tdof
                  ! s_dim = dof*ndim
                  ll = (l -1) * s_dim
                  rr = (r -1) * s_dim

                  tval = T_rule%phi(l,alpha)*T_rule%Dphi(r,alpha) &
                          * T_rule%weights(alpha) * local_eta

                  nd1 = 1
                  do d1 =1,ndim
                    nd2 = 1
                    do d2 = 1,ndim
                      elem%bigBlock(0)%Mb( 1:dof, 1:dof, d1, d2, l, r ) = &
                        elem%bigBlock(0)%Mb( 1:dof, 1:dof, d1, d2, l, r ) &
                        + tval * T_mat(nd1:nd1+dof-1,nd2:nd2+dof-1)
                      nd2 = nd2 + dof
                    end do
                    nd1 = nd1 + dof
                  end do
                enddo  ! r
               enddo ! l

             end if
!               print*, 'norm of bigBlock 0 =' , norm2(elem%bigBlock(0)%Mb)
!               print*, 'bigBlock  6,6 =' , elem%bigBlock(0)%Mb(6,6,1,1,1,1)
            enddo !i

         end do !alpha =  1,Qdeg

         if(state%model%precomputed_arrays ) then  ! we have to precompute the  values
           stop "precomputed_arrays are not done in ComputeST_Terms!"
         endif

         if(state%model%varying_time_term) then
          do i = 1, grid%nelem
            elem => grid%elem(i)

            Tdof = elem%getActualTDof()
            dof = elem%getActualDof()
            s_dim = ndim *dof

            ! save the wST space-time solution in quadrature index alpha to w
            call Transfer_wST_to_w_Elem(elem , -1, Qdeg)
            call Transfer_wST_to_wActual_Elem(elem , -1, Qdeg)

             ! computing of the matrix in front of the time derivative term at t_{m-1}^+
             call Compute_Time_deriv(elem, .false., s_dim, T_mat(1:s_dim,1:s_dim) )

             do l=1,Tdof
                do r = 1, Tdof
                  ll = (l -1) * s_dim
                  rr = (r -1) * s_dim

                  tval = T_rule%phi(l,-1)*T_rule%phi(r,-1) * local_eta

                  nd1 = 1
                  do d1 =1,ndim
                    nd2 = 1
                    do d2 = 1,ndim
                      elem%bigBlock(0)%Mb( 1:dof, 1:dof, d1, d2, l, r ) = &
                        elem%bigBlock(0)%Mb( 1:dof, 1:dof, d1, d2, l, r ) &
                        + tval * T_mat(nd1:nd1+dof-1,nd2:nd2+dof-1)
                      nd2 = nd2 + dof
                    end do
                    nd1 = nd1 + dof
                  end do
                enddo  ! r
             enddo ! l
          enddo ! i=1,grid%nelem
         end if


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! IMPLICITLY = FALSE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      else   ! implicitly == .false.

        state%cons = 0. ! auxiliary array for verifying the code

        iprint = .false.

        call ClearVectorBlocks_STDGM()  ! rhsST = 0


        cTime = state%time%ctime  ! storing of the actual value

        itest = -13   ! 9

        do alpha = 1, Qdeg ! temporarily max_Tdof =  max time quadrature nodes

            ! verifying the consistency of the scheme
            state%cons(0,1) = T_rule%weights(alpha)

            do i = 1, grid%nelem
              elem => grid%elem(i)
              if (wSTplus) then
                call Transfer_wSTplus_to_wActual_Elem(elem , alpha, Qdeg)
              else
                call Transfer_wST_to_wActual_Elem(elem , alpha, Qdeg)
              endif
            enddo

            !we have to run ComputeTerms() in the time-quadrature nodes
            state%time%ctime = state%time%ttime + state%time%tau(1) * state%time%T_rule(Qdeg)%lambda(alpha)

            call ComputeTerms( deg_plus )

            do i =1, grid%nelem
               elem => grid%elem(i)


               !write(*,'(a10, i5, 300es12.4)') '!! ELM:', i, elem%vec(rhs, :)

               !NEW for adaptation
               wdof = elem%dof
               wTdof = elem%Tdof

               !NEW for adaptation
!               if(elem%deg_plus) then
!                  dof = elem%dof_plus
!                  Tdof = elem%Tdof_plus
!               else
!                  dof = elem%dof
!                  Tdof = elem%Tdof
!               endif
               dof = elem%getActualDof()
               Tdof = elem%getActualTDof()

               s_dim = ndim * dof

               !vec -> rhsST
               do k = 1,ndim
                  do l = 1, Tdof
                     val_vec = T_rule%phi(l,alpha) * T_rule%weights(alpha)

                     elem%rhsST(k, 1:dof, l) = elem%rhsST(k, 1:dof, l)  &
                          +  val_vec *  elem%vec( rhs, (k-1)*dof + 1 : k*dof )   ! dofA

                     if(k==1 .and. l == 1) state%cons(1,1) = state%cons(1,1)&
                          + val_vec *  elem%vec( rhs, (k-1)*dof + 1 )

                  enddo !l
               enddo !k

               ! adding of the term in front of the time derivative - time deriv term itself
               if(state%model%varying_time_term) then

                  ! computing of the matrix in front of the time derivative term at t_alpha
                  !print*,'Compute_Time_deriv C', s_dim, elem%dof_plus
                  call Compute_Time_deriv(elem, elem%deg_plus, s_dim,  T_mat(1:s_dim,1:s_dim) )
                  !compute the vector in the time derivative
                  call Compute_Time_deriv_Vector(elem, dof,  T_vec(1: s_dim) )


                  do k = 1,ndim
                     do l = 1, Tdof

                        val_vec = T_rule%phi(l,alpha) * T_rule%weights(alpha) * local_eta

                        val_vecD = T_rule%Dphi(l,alpha) * T_rule%weights(alpha) * local_eta

                        do j=1, dof

                           if(.not. CD_tdt) then
                              ! CD_tdt -
                              elem%rhsST(k, j, l) = elem%rhsST(k, j, l)  &
                                   -  val_vec  &
                                   * dot_product(T_rule%Dphi(1:wTdof,alpha) , &
                                   matmul( T_mat( (k-1)* dof + j, (k-1)*dof+1:(k-1)*dof + wdof ),&
                                   elem%wST(k, 1:wdof, 1:wTdof) ) )
                              ! CD_tdt -
                           else
                              ! CD_tdt +
                              elem%rhsST(k, j, l) = elem%rhsST(k, j, l) + val_vecD * T_vec( (k-1)* dof +j)
                              ! CD_tdt +
                           endif


                        enddo
                     enddo  ! l=1,Tdof
                  enddo ! k=1,ndim
               endif ! state%model%varying_time_term
            enddo !i

         end do !alpha =  1,Qdeg


         ! adding the time derivative and the the w_{m-1}^+ part:
         if(state%model%precomputed_arrays ) then  ! we have to precompute the  values
            do i = 1, grid%nelem
               elem => grid%elem(i)

               ! save the wST space-time solution in quadrature index alpha to w
              if (wSTplus) then
                call Transfer_wSTplus_to_wActual_Elem(elem , -1, Qdeg)
              else
                call Transfer_wST_to_wActual_Elem(elem , -1, Qdeg)
              endif
            enddo
            !if(state%nlSolver%iter <= 5) &
            call ComputeCapacityConductivity(.true.) ! DO NOT comment, used in inv_fluxes !!!!!!!!!!
            !write(*,'(a7,200es12.4)')'CAP Y:',grid%elem(741)%xi(0, 1:3, iRe+4),grid%elem(741)%xi(0, 1:3, iRe+3)

         endif

         !val1 = 0.

         do i =1, grid%nelem
            elem => grid%elem(i)

            !NEW for adaptation
            wdof = elem%dof
            wTdof = elem%Tdof
!
!            !NEW for adaptation
!            if(elem%deg_plus) then
!               dof = elem%dof_plus
!               Tdof = elem%Tdof_plus
!            else
!               dof = elem%dof
!               Tdof = elem%Tdof
!            endif
            dof = elem%getActualDof()
            Tdof = elem%getActualTDof()

            ! original subroutine for the adding the time derivative and the the w_{m-1}^+ part
            if(.not.  state%model%varying_time_term ) then
               do k = 1, ndim
                  do l = 1,Tdof
                     do j = 1,dof   ! dofA
                        !NEW for adaptation
                        print*, "problem here - wST/wSTplus"
                        val_vec = local_eta * dot_product( state%time%refTimeMatrix%Mb(l,1:wTdof) , &
                             matmul( elem%Mass%Mb(j, 1:wdof), elem%wST(k,1:wdof,1:wTdof) ))

                        elem%rhsST(k, j, l) = elem%rhsST(k, j, l) - val_vec

                     enddo !j
                  enddo !l

               enddo !k
            endif

            ! adding of the term in front of the time derivative -  jump term,  w_{m-1}^+ part:
            if(state%model%varying_time_term) then
              if (wSTplus) then
                call Transfer_wSTplus_to_wActual_Elem(elem , -1, Qdeg)
              else
                call Transfer_wST_to_wActual_Elem(elem , -1, Qdeg)
              endif
               ! computing of the matrix in front of the time derivative term at t_{m-1}^+
               s_dim = ndim * dof

               !compute the vector in the time derivative
               call Compute_Time_deriv_Vector(elem, dof,  T_vec(1: s_dim) )


               do k = 1,ndim
                  do l = 1, Tdof

                     val_vec = T_rule%phi(l,-1) * local_eta

                     do j=1, dof
                        ! formula without linearization
                        if(.not. CD_tdt) then
                           ! CD_tdt -
                           elem%rhsST(k, j, l) = elem%rhsST(k, j, l) - val_vec * T_vec( (k-1)*dof +j )
                        endif
                     enddo
                  enddo  ! l=1,Tdof
               enddo ! k=1,ndim


               ! !! VARIANT CD_tdt +
               if (wSTplus) then
                 call Transfer_wSTplus_to_wActual_Elem(elem , 0, Qdeg)
               else
                call Transfer_wST_to_wActual_Elem(elem , 0, Qdeg)
               endif

               ! computing of the matrix in front of the time derivative term at t_{m-1}^+
               s_dim = ndim * dof
               !print*,'Compute_Time_deriv D', s_dim, elem%dof_plus
               call Compute_Time_deriv(elem, elem%deg_plus, s_dim, T_mat(1:s_dim, 1:s_dim) )

               !compute the vector in the time derivative
               call Compute_Time_deriv_Vector(elem, dof,  T_vec(1: s_dim) )

               do k = 1, ndim
                  do l = 1, Tdof

                     val_vec = T_rule%phi(l, 0) * local_eta

                     do j=1, dof

                        ! CD_tdt -  .. is empty
                        if(CD_tdt) then
                           ! CD_tdt +
                           elem%rhsST(k, j, l) = elem%rhsST(k, j, l) - val_vec * T_vec( (k-1)*dof +j )
                        endif
                     enddo
                  enddo  ! l=1,Tdof
               enddo ! k=1,ndim
               ! !! END OF THE TEST


            endif   !state%model%varying_time_term
         enddo !i

    !         write(*,'(a15, 12es12.4)') 'F^m: aft mass' , grid%elem(1)%rhsST(1,:,:)

         ! the w_{m-1}^- part:   "state%model%varying_time_term" given inside

         if(state%model%precomputed_arrays ) then  ! we have to precompute the  values
            do i = 1, grid%nelem
               elem => grid%elem(i)

               ! save the wST space-time solution in quadrature index alpha to w
              if (wSTplus) then
                call Transfer_wSTplus_to_wActual_Elem(elem , -1, Qdeg)
              else
                call Transfer_wST_to_wActual_Elem(elem , -1, Qdeg)
              endif
            enddo
            !if(state%nlSolver%iter <= 5) &
            call ComputeCapacityConductivity(.true.) ! DO NOT comment, used in inv_fluxes !!!!!!!
            !write(*,'(a7,200es12.4)')'CAP Z:',grid%elem(741)%xi(0, 1:3, iRe+4),grid%elem(741)%xi(0, 1:3, iRe+3)

         endif


         !val = 0.
         do i = 1, grid%nelem
            elem => grid%elem(i)

            ! state%cons(4,1) filled inside
            print*, "control here - wST vs wSTplus in wSTfin"
            call Elem_wSTfinToRhsST( grid%elem(i) , state%nlSolver%implicitly)

         enddo !i

      endif   !state%nlSolver%implicitly



    ! putting back the original value
    state%time%ctime = cTime

    if(state%model%varying_time_term) deallocate( T_mat, T_vec)

    !deallocate(cons_test)
    call state%setP_mod(p_mod_loc)
    call state%setQ_mod(q_mod_loc)

    call cpu_time(tt)

    !write(*,'(a10, i5, 200es12.4)') 'STterms w:', grid%elem(1)%i,  grid%elem(1)%rhsST

  end subroutine ComputeST_Terms_new

  !> evaluation of the LOCAL PART (DIAGONAL BLOCK) of thenew big matrix and vector blocks for STDGM method
 !> if implicitly == False -> vector F(u) saved in elem%vec(rhs,:)
 !> if implicitly == True  -> sets the matrix C_K(*,*) elem%vec(rhs = q(u),
 !> new version:
 !> if wSTritz == TRUE -> wSTritz -> elem%wActual
 !> else wST -> elem%wActual
  subroutine ComputeSTDiag_Terms_loc(elemm, deg_plus, wSTRitz )
    class(element), intent(inout) :: elemm
    logical, intent (in) :: deg_plus
    logical, intent (in) :: wSTRitz
    class(element), pointer :: elem1
    class(Time_rule), pointer :: T_rule
    real, dimension(:,:), allocatable :: T_mat
    real, dimension(:), allocatable :: T_vec
    integer :: i, j, k, l, r, kk, rr, ll
    integer :: alpha, Qdeg, dof, Tdof, s_dim, f_dof, f_Tdof, f_dim
    integer :: m, mm, n, nn, wDof, wTdof
    real :: cTime, tt, tt1
    real :: val, val_vec, val_vecD, val0, val1, val2, val3, val_sum
    real :: local_eta, tval
    logical ::iprint, CD_tdt
    integer :: itest
    integer :: nd1, nd2, d1, d2, p_mod_loc, q_mod_loc

!    print*, "ComputeSTDiag_Terms_loc called wit impl: ", state%nlSolver%implicitly

    p_mod_loc = state%getP_mod()
    q_mod_loc = state%getQ_mod()

    if (.not. deg_plus .and. state%getP_mod() /= 0) &
      stop "deg_plus=FALSE && p_mod > 0 in ComputeST_Terms_new"

    if ( deg_plus ) then
      call state%setP_mod( state%p_mod_max )
      call state%setQ_mod( state%q_mod_max )
      grid%elem(:)%deg_plus = .true.
    else
      grid%elem(:)%deg_plus = .false.
    end if

    !print*, 'ComputeST_Terms with implicitly = ' ,state%nlSolver%implicitly, 'deg_plus:', deg_plus
    call cpu_time(tt1)

    cTime = state%time%ctime
    local_eta = 1. / state%time%tau(1)

    ! adding of the term in front of the time derivative
    if(state%model%varying_time_term) then
      stop "local ComputeTerms do not work for time varying term"
    endif

    !NEW
    ! set degree with respect to time
    Qdeg = state%time%Qnum

    ! TODO FR - Qdeg should be connected with q_mod_max !
    T_rule => state%time%T_rule(Qdeg)

    associate ( time => state%time)
    select type ( time )
      class is ( TimeTDG_t )
      class default
        stop 'computeST_Terms only for STDG'
    end select
    end associate




!!! IMPLICITLY = TRUE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         ! rhsST = q(w)
      if(state%nlSolver%implicitly) then

        ! clear the diagonal block
        call elemm%bigBlock(0)%zero()

        do alpha = 1, Qdeg ! temporarily max_Tdof =  max time quadrature nodes
              if (wSTRitz) then
                call Transfer_wSTRitz_to_wActual_Elem(elemm , alpha, Qdeg)
              else
                call Transfer_wST_to_wActual_Elem(elemm , alpha, Qdeg)
              end if



!              ! for Ritz test
!              if (elemm%i==1) then
!                print*, "Transfering on the neighbors COmmented - TEST!!!"
!              end if
!              do i = 1,  elemm%flen
!                  if ( elemm%face(neigh,i) > 0 ) then
!                    elem1 => grid%elem( elemm%face(neigh,i) )
!                    if (wSTRitz) then
!                      call Transfer_wSTplus_to_wActual_Elem(elem1 , alpha, Qdeg)
!                    else
!                      call Transfer_wST_to_wActual_Elem(elem1 , alpha, Qdeg)
!                    end if
!                  end if
!              end do

              ! wST is needed on neighbors for local BC
              do i = 1,  elemm%flen
                  if ( elemm%face(neigh,i) > 0 ) then
                    elem1 => grid%elem( elemm%face(neigh,i) )
                    call Transfer_wST_to_wActual_Elem(elem1 , alpha, Qdeg)
                  end if
              end do

            !we have to run ComputeTerms() in the time-quadrature nodes
            !cTime = state%time%ctime
            state%time%ctime = state%time%ttime + state%time%tau(1) * T_rule%lambda(alpha)
            ! local compute terms
            call Compute_ONLY_ONE_ELEMENT_Terms( elemm )

            state%time%ctime = cTime

            Tdof = elemm%getActualTDof()
            dof = elemm%getActualDof()
            s_dim = ndim *dof

            do l = 1, Tdof
                !diag blocks of bigBlock
              do r = 1, Tdof
                tval = T_rule%phi(l,alpha)*T_rule%phi(r,alpha) * T_rule%weights(alpha)
                nd1 = 1
                do d1 = 1,ndim
                  nd2 = 1
                  do d2 = 1, ndim
                    elemm%bigBlock(0)%Mb( 1:dof, 1:dof, d1, d2, l, r ) = &
                      elemm%bigBlock(0)%Mb( 1:dof, 1:dof, d1, d2, l, r )  + &
                      elemm%block(0)%Mb(nd1:nd1+dof-1,nd2:nd2+dof-1) * tval
                    nd2 = nd2 + dof
                  end do
                  nd1 = nd1 + dof
                end do ! d1
              enddo !r
           enddo !l = 1,Tdof


           ! adding of the term in front of the time derivative - time deriv term itself
           if(state%model%varying_time_term) then
              stop 'FR - time varying time term needs to be done later'
           endif ! state%model%varying_time_term

         end do !alpha =  1,Qdeg

         ! adding of the term in front of the time derivative - jump term,  w_{m-1}^+ part:
         if(state%model%varying_time_term) then
            stop 'FR - time varying time term 2 needs to be done later'
         endif !(state%model%varying_time_term


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! IMPLICITLY = FALSE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      else   ! implicitly == .false.

        iprint = .false.

        dof = elemm%getActualDof()
        Tdof = elemm%getActualTDof()
        s_dim = ndim * dof

        ! control size
        if (size(elemm%rhsST) < ndim*dof*Tdof) then
          print*, "elem%dof, dof, ndim, Tdof:", elemm%dof, dof, ndim, Tdof
          print*, "size of rhsST:", size(elemm%rhsST)
          stop "elem%rhsST is too small in ComputeST_Terms_loc"
        endif

        ! clear vector blocks
        elemm%rhsST(:,:,:) = 0.0

        cTime = state%time%ctime  ! storing of the actual value

        itest = -13   ! 9

        do alpha = 1, Qdeg ! temporarily max_Tdof =  max time quadrature nodes
          if (wSTRitz) then
              call Transfer_wSTRitz_to_wActual_Elem(elemm , alpha, Qdeg)
            else
              call Transfer_wST_to_wActual_Elem(elemm , alpha, Qdeg)
          end if

! test Ritz with updating the neighbors
!          if (elemm%i==1) then
!                print*, "Transfering on the neighbors COmmented - changed to wSTRitz!!!"
!          end if
!
!          ! wST is needed on neighbors for local BC
!          do i = 1,  elemm%flen
!              if ( elemm%face(neigh,i) > 0 ) then
!                elem1 => grid%elem( elemm%face(neigh,i) )
!                if (wSTRitz) then
!                  call Transfer_wSTplus_to_wActual_Elem(elem1 , alpha, Qdeg)
!                else
!                  call Transfer_wST_to_wActual_Elem(elem1 , alpha, Qdeg)
!                end if
!              end if
!          end do

          !  wST is needed on neighbors for local BC
          do i = 1,  elemm%flen
              if ( elemm%face(neigh,i) > 0 ) then
                elem1 => grid%elem( elemm%face(neigh,i) )
                call Transfer_wST_to_wActual_Elem(elem1 , alpha, Qdeg)
              end if
          end do

          !we have to run ComputeTerms() in the time-quadrature nodes
          state%time%ctime = state%time%ttime + state%time%tau(1) * state%time%T_rule(Qdeg)%lambda(alpha)

          call Compute_ONLY_ONE_ELEMENT_Terms( elemm )

          !vec -> rhsST
          do k = 1,ndim
            do l = 1, Tdof
               val_vec = T_rule%phi(l,alpha) * T_rule%weights(alpha)

               elemm%rhsST(k, 1:dof, l) = elemm%rhsST(k, 1:dof, l)  &
                    +  val_vec *  elemm%vec( rhs, (k-1)*dof + 1 : k*dof )   	! dofA
            enddo !l
          enddo !k
        end do !alpha =  1,Qdeg

        if(state%model%precomputed_arrays ) then  ! we have to precompute the  values
          stop "ComputeST_Terms_loc does not work for state%model%precomputed_arrays"
        endif

        ! adding the time derivative and the the w_{m-1}^+ part:
        if (wSTRitz) then
          ! makes sense onl if deg_plus i.e. wSTplus has dof == wdof
          wdof = elemm%getActualDof() !elemm%dof_plus ! the dimension of the function wSTplus
          wTdof = elemm%getActualTDof()
          ! adding the time derivative and the the w_{m-1}^+ part
          do k = 1, ndim
            do l = 1,Tdof
               do j = 1,dof   ! dimension of the residual
                  !NEW for adaptation
                  val_vec = local_eta * dot_product( state%time%refTimeMatrix%Mb(l,1:wTdof) , &
                    matmul( elemm%Mass%Mb(j, 1:wdof), elemm%wST_Ritz(k,1:wdof,1:wTdof) ))

                  elemm%rhsST(k, j, l) = elemm%rhsST(k, j, l) - val_vec

               enddo !j
            enddo !l

          enddo !k

        else
          wdof = elemm%dof
          wTdof = elemm%Tdof
          ! adding the time derivative and the the w_{m-1}^+ part
          do k = 1, ndim
            do l = 1,Tdof
               do j = 1,dof   ! dofA
                  !NEW for adaptation
                  val_vec = local_eta * dot_product( state%time%refTimeMatrix%Mb(l,1:wTdof) , &
                       matmul( elemm%Mass%Mb(j, 1:wdof), elemm%wST(k,1:wdof,1:wTdof) ))
                  elemm%rhsST(k, j, l) = elemm%rhsST(k, j, l) - val_vec
               enddo !j
            enddo !l
          enddo !k
        end if
        ! here is always wST since it is initial condition from previous time step
        call Elem_wSTfinToRhsST( elemm , state%nlSolver%implicitly)

      endif   !state%nlSolver%implicitly

    ! putting back the original value
    state%time%ctime = cTime

    !deallocate(cons_test)
    call state%setP_mod(p_mod_loc)
    call state%setQ_mod(q_mod_loc)

    call cpu_time(tt)

  end subroutine ComputeSTDiag_Terms_loc

  !> the sizes depend on state%p_mod, state%q_mod
  !> since dualResVector is always of size p_mod_max it should be called only when
  !> p_mod == p_mod_max
  !> plus -> zSTplus, !plus -> zST is used
  subroutine computeDualResidualPlus( grid, DWR, plus )
    type( DWR_t ), intent(inout) :: DWR
    class( mesh ), intent(inout) :: grid
    !integer, intent(in) :: p_mod, q_mod
    logical, intent(in) :: plus
    real, allocatable, dimension(:) :: xx ,bb
    integer :: bigNsize
    logical :: loc_implicitly

!    print*, "### computeDualResidualPlus called with plus = ", plus

    if (state%getP_mod() /= state%p_mod_max) &
      stop "computeDualResidualPlus should be called only when p_mod==p_mod_max!"

    bigNsize = state%bigNSize( state%getP_mod(), state%getQ_mod() )
    allocate( xx(1:bigNsize), source = 0.0 )
    allocate( bb(1:bigNsize), source = 0.0 )

    if (plus) then
      call copyZSTplus_toLongVectorMod( grid, state%getP_mod(), state%getQ_mod(), xx )
    else
      call copyZST_toLongVectorMod( grid, state%getP_mod(), state%getQ_mod(), xx )
    end if

    !print*, "control = ", norm2(DWR%x), norm2(xx) ! this seems working OK now

    if ( norm2( grid%elem(1)%bigBlock(0)%Mb ) < 1.E-12 ) &
      stop 'bigBlock is zero in computeDualResidualPlus'

    ! bb = A(plus)^T * zz
    call bMVprodBIG_Dual( bb, xx, bigNsize )

    call DWR%fillDualResVector( grid, bigNsize, bb(1:bigNsize) )
    deallocate( xx, bb )

  end subroutine computeDualResidualPlus





end module terms_mod
