!> subroutine of evaluation of some quantities, fluxes for porous media flow
!> the main iterative loop
module pm_fluxes
  use main_data
!  use problem_oper
  use eval_sol
  use set_solution
  use stdgm_mod
  use mesh_mod
  use inviscid_fluxes
  use emiss_bc
  
  implicit none

    ! porous media flow
  public:: ComputeCapacityConductivity
  public:: ComputeCapacityConductivity_2

  public:: Eval_Characteristic_paramets
  public:: Elem_water_content

contains



    !> POROUS MEDIA FLOW, compute the HO continuous reconstruction of the capacity and conductivity
  subroutine ComputeCapacityConductivity( onlyTimeTerms )
    logical, intent (in) :: onlyTimeTerms ! is true, do not compute conductivity
    class(element), pointer :: elem
    real, dimension(:,:), allocatable :: wi
    real, dimension(:,:,:), allocatable :: Dwi
    real, dimension(:,:,:), allocatable :: TA, Ksk
    type(volume_rule), pointer :: V_rule
    type(Gauss_rule), pointer :: G_rule
    real, dimension(:,:), allocatable :: Re_1
    real, dimension(:,:), pointer :: phi
    integer :: itype, deg1, Qdof, Qnum, i,j, j1, ifile, l, dof, Ldof, ist, k, k1, ie
    logical :: conforming_P1

    return

    conforming_P1 = .false.
    !conforming_P1 = .true.

    !print*,'#  ComputeCapacityConductivity called',onlyTimeTerms
    !return

    !if(state%nlSolver%iter == 6 .and. onlyTimeTerms .and. state%nlSolver%implicitly) &
    !     print*,"  if(state%nlSolver%iter > 5) return"
    !if(state%nlSolver%iter > 5) return

    ! continuous reconstruction of the actual solution
    !call ComputeHO_LocalProblems( 1 )

    ! conforming P1 projected solution given by the vertex interpolation
    if(conforming_P1)  &
         call Compute_conforming_P1_projection( grid%npoin, porous_data%wR(1:grid%npoin, 0:ndim) )

    !stop "u493uj43i"

    !open(23, file ='conductivity', status='unknown')

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

       ! volume integ nodes
       Qdof = elem%Qdof
       allocate(wi( 1:Qdof, 1:ndim), Dwi(1:Qdof,1:ndim, 1:nbDim))

       call Eval_w_Elem(elem, wi(1:Qdof,1:ndim) )

       call Eval_Dw_Elem(elem, Dwi(1:Qdof, 1:ndim, 1:nbDim) )


       if(conforming_P1) then
          ! the use of the values from the Compute_conforming_P1_projection
          V_rule => state%space%V_rule(elem%Qnum)

          wi = 0.
          do j=1, 3
             j1 = mod(j, 3) + 1
             k = elem%face(idx, j1)
             do l=1, Qdof
                wi(l, 1:ndim) =  wi(l, 1:ndim) + V_rule%lambda(l, j) * porous_data%wR(k, 1:ndim)
             enddo
          enddo
          !!wi(1:Qdof,1:ndim) = wi(1:Qdof,1:ndim) / 2.
       endif

       allocate(  TA(1:Qdof, 1:ndim, 1:ndim ) )

       ! capacity
       call Set_Time_Matrix_porous(elem, ndim,  Qdof, wi( 1:Qdof, 1:ndim), &
            elem%xi(0,1:Qdof, 1:2+iRe),  TA(1:Qdof, 1:ndim, 1:ndim) )

       !if(state%nlSolver%iter <= 5) &
       elem%xi(0, 1:Qdof, nbdim+iRe + 2 ) =  TA(1:Qdof, 1, 1)

       ! if(elem%i == 1724) then
       !      write(100*state%time%iter + 191,'(a5,i5,200es16.8)')'CAP_W:', state%nlSolver%iter, &
       !      elem%xi(0, :, iRe+4)
       !      write(100*state%time%iter + 192,'(a5,i5,200es16.8)')'CAP_W:', state%nlSolver%iter, &
       !           wi(:, 1)
       !      write(100*state%time%iter + 193,'(a5,i5,200es16.8)')'CAP_W:', state%nlSolver%iter, &
       !           wi(:, 1) - elem%xi(0,:, 2)
       !   endif


       ! if( abs( elem%xc(1) ) < 1. .and. abs(elem%xc(2) ) < 1. ) then

       !    ifile = 500+state%nlSolver%iter
       !    do k=1,Qdof
       !       write(ifile, *)  elem%xi(0, k, 1:2 ), wi(k, 1), elem%xi(0, k, 2+iRe+2 )
       !    enddo
       !    write(ifile, *) '  '
       !    write(ifile, *) ' ## #  #   #   # '
       !    write(ifile, *) '  '
       ! endif

       if( .not.  onlyTimeTerms ) then
          ! conductivity
          allocate(Re_1(1:iRe, 1:Qdof) )
          Re_1(1, 1:Qdof) = 1./state%model%Re
          Re_1(2:iRe, 1:Qdof) = transpose( elem%xi(0, 1:Qdof, 2+1:2+iRe-1) )

          allocate( Ksk(1:Qdof, 1:nbdim, 1:nbDim) )
          !if(elem%i == 22) write(31,'(a6, 200es12.4)') 'elem%w:',elem%w(0,:)
          do k=1,Qdof
             call Eval_Diff_Porous_Coeffs(wi(k,1), Dwi(k, 1, 1:nbDim), Ksk(k, 1:nbDim, 1:nbDim), &
                  Re_1(1:iRe, k), 0, elem%xi(0, k, 1:nbDim) )

             !if(elem%i == 22) then
             !   write(31,'(a6,2i5, 200es12.4)') 'xii:',elem%i, k, &
             !        wi(k, :), Dwi(k, 1, :),  Ksk(k,1, 1)
             !endif

          enddo

          !if(state%nlSolver%iter <= 5) &
          elem%xi(0, 1:Qdof, nbdim+iRe + 1 ) =  Ksk(1:Qdof, 1, 1)

          deallocate( Ksk, Re_1)
       endif

       deallocate(wi, Dwi, Ta)

       if( .not.  onlyTimeTerms ) then
          ! edge values
          do ie = 1, elem%flen

             Qdof = elem%face(fGdof,ie)

             allocate( wi(1:Qdof,1:ndim), Dwi(1:Qdof,1:ndim,1:nbDim) )
             call Eval_w_Edge(elem, ie,  wi, .false.)
             call Eval_Dw_Edge(elem,  ie,  Dwi(1:Qdof,1:ndim,1:nbDim),  .false.)

             if(conforming_P1) then
                G_rule => state%space%G_rule(elem%face(fGnum, ie) )

                j =  ie
                j1 = mod( j, 3) + 1
                k = elem%face(idx, j)
                k1 = elem%face(idx, j1)
                do l=1, Qdof
                   wi(l, 1:ndim) =  (1 - G_rule%lambda(l)) * porous_data%wR(k, 1:ndim) &
                        +  G_rule%lambda(l)  * porous_data%wR(k1, 1:ndim)
                   !   !write(71, *)  elem%xi(ie, l, 1:nbDim), wi(l, 1:ndim)
                enddo
             endif

             ! conductivity
             allocate(Re_1(1:iRe, 1:Qdof) )

             Re_1(1, 1:Qdof) = 1./state%model%Re
             Re_1(2:iRe, 1:Qdof) = transpose( elem%xi(ie, 1:Qdof, 2+1:2+iRe-1) )


             allocate( Ksk(1:Qdof, 1:nbdim, 1:nbDim) )

             !if(elem%i == 22) write(31,'(a6, 200es12.4)') 'elem%w:',elem%w(0,:)
             !if(elem%i == 357) then
             !   write(*,'(a6, i5, 200es12.4)') 'Wii:', ie, wi(:, 1) !, Dwi(k, 1, :),  Ksk(k,1, 1)
             !endif

             do k=1,Qdof


                call Eval_Diff_Porous_Coeffs(wi(k,1), Dwi(k, 1, 1:nbDim), Ksk(k, 1:nbDim, 1:nbDim), &
                     Re_1(1:iRe, k), 0, elem%xi(ie, k, 1:nbDim) )

                !if( abs( elem%xi(ie, k, 1) ) < 0.1 .and. abs(elem%xi(ie, k, 2) ) < 0.1 ) then
                !   write(*,'(a6, 2i5, 200es12.4)') 'Ksk:', elem%i, ie, Re_1(1:iRe, k), elem%xi(ie, k, 1:nbDim),&
                !        Ksk(k, 1, 1), Ksk(k,2,2)
                !endif


             enddo

             !if(state%nlSolver%iter <= 5) &
             elem%xi(ie, 1:Qdof, nbdim+iRe + 1 ) =  Ksk(1:Qdof, 1, 1)

             ! if( abs( elem%xc(1) ) < 1. .and. abs(elem%xc(2) ) < 1. ) then
             !    do k=1,Qdof
             !       write(100+state%nlSolver%iter, *)  elem%xi(ie, k, 1:2 ), wi(k,1), &
             !            elem%xi(ie, k, 2+iRe+1 )
             !    enddo
             !    write(100+state%nlSolver%iter, *) '  '
             !    write(100+state%nlSolver%iter, *) '### de3yd83hd  ', &
             !         elem%i, state%nlSolver%implicitly, state%nlSolver%lambda

             !    write(100+state%nlSolver%iter, *) '  '
             ! endif


             deallocate(wi, Dwi, Ksk, Re_1)

          enddo  ! ie = 1, elem%flen

       endif  ! onlyTimeTerms

    enddo ! do i=1,grid%nelem

    ! FR chybelo end do ?
    !print*, 'ComputeCapacityConductivity: '
    !print*, 'FR: I added missing end do - maybe also something more is missing???'

  end subroutine ComputeCapacityConductivity

  !> POROUS MEDIA, detection of element intersection zero level of the pressure head
  subroutine Elem_detect_intersect(elem)
    class(element), intent(inout) :: elem
    real, dimension(:,:), allocatable :: wi
    integer :: je, Qdof

    elem%intersect_zero = .false.

    do je=0, elem%flen  ! volume & edge integ nodes

       if(je ==0) then      ! volume integ nodes
          Qdof  = elem%Qdof
          allocate (wi(1:Qdof, 1:2*ndim) )

          ! H
          call Eval_w_Elem(elem, wi(1:Qdof, 1:ndim) )

          ! h = H - z
          wi(1:Qdof, 2) = wi(1:Qdof, 1) - elem%xi(0, 1:Qdof, 2)

       else                 ! edge integ nodes
          Qdof = elem%face(fGdof,je)
          allocate (wi(1:Qdof, 1:2*ndim) )

          ! H
          call Eval_w_Edge(elem, je, wi(1:Qdof, 1:ndim), .false.)

          ! h = H - z
          wi(1:Qdof, 2) = wi(1:Qdof, 1) - elem%xi(je, 1:Qdof, 2)

       endif

       ! detection of the intersection
       if(maxval( wi(1:Qdof, 2) ) > 0. .and. minval( wi(1:Qdof, 2) )< 0.) then
             elem%intersect_zero = .true.

          !do l=1,Qdof
          !   write(20, *) elem%xi(je, l, 1:2), wi(l,1:2), &
          !        maxval( wi(1:Qdof,2) ), minval( wi(1:Qdof,2) ), '  73y73h'
          !enddo
       endif

       deallocate(wi)

    enddo  ! je = 0, elem%flen
  end subroutine Elem_detect_intersect


  !> POROUS MEDIA FLOW, compute the capacity and conductivity
  !> limitation for the intersection
  subroutine ComputeCapacityConductivity_2( )
    class(element), pointer :: elem
    real, dimension(:,:), allocatable :: wi, ww
    real, dimension(:,:), pointer :: phi
    type(Lagrang_rule), pointer :: L_rule
    type(volume_rule), pointer :: V_rule
    type(Gauss_rule), pointer :: G_rule
    real, dimension(:,:), allocatable :: Re_1
    real :: weight
    integer :: Qdof, i, dof, l, je, Fdof, Qnum,k

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

       call Elem_detect_intersect(elem)

       ! setting of the solution used later for the computation of the  capacity and conductivity
       ! for elements detected by elem%intersect_zero
       if(elem%intersect_zero ) then

          ! P_1 projection of the solution w in element vertexes
          Qnum = 1

          L_rule => state%space%L_rule(Qnum)
          Qdof = L_rule%Qdof

          phi => L_rule%phi(1:dof, 1:Qdof)

          allocate( wi(  1:ndim, 1: Qdof))
          allocate( ww(  1:ndim, 1: elem%dof))

          ! values in Lagrangian nodes (element vertexes)
          do k=1,ndim
             wi(k, 1:Qdof) = matmul( elem%w(0, (k-1)*dof+1 : k*dof), phi(1:dof, 1:Qdof) )
          enddo

          !do l=1,Qdof
          !   write(21, *) grid%x(elem%face(idx, l), 1:2), wi(1,l), '  7de333h'
          !enddo

          !!!write(*,'(a8, 2i5, 30es12.4)') 'wi:',Qnum,Qdof, wi(1, :)
          ! recomputation in DG basis functions
          call Lagr2BasisDiff(elem, Qnum,  wi(1:ndim, 1:Qdof), elem%Qnum, dof, ww(1:ndim, 1:dof) )

          !!!write(*,'(a8, 2i5, 30es12.4)') 'ww:',Qnum,dof, ww(1, :)

          !call PlotElemFunction3D(22, elem, dof, ww(1, 1: dof) )

          ! averaging of the original function and its projection
          weight = 1.0  !0.75

          do k=1, ndim
             ww(k, 1: dof) = weight * ww(k, 1: dof) + (1.-weight) * elem%w(0, (k-1)*dof+1 : k*dof)
          enddo

          !call PlotElemFunction3D(23, elem, dof, ww(1, 1: dof) )

          ! projected solution in integ nodes
          Fdof = max(elem%Qdof, maxval( elem%face(fGdof, :)  ) )

          if( allocated( elem%w_limit ) ) deallocate( elem%w_limit )
          allocate( elem%w_limit(0:elem%flen, 1:Fdof, 1:ndim) )

          ! recomputation of the projection in integ nodes
          do je=0, elem%flen  ! volume & edge integ nodes

             if(je ==0) then      ! volume integ nodes
                Qdof  = elem%Qdof
                phi => state%space%V_rule(elem%Qnum)%phi(1:dof, 1:Qdof)

             else ! boundary edge
                Qnum = elem%face(fGnum,je)
                Qdof = state%space%G_rule(Qnum)%Qdof

                ! the reference test functions
                if(elem%HGnode) then
                   phi => state%space%G_rule(Qnum)%phi(elem%type, elem%HGface(1,je), elem%HGface(2, je),&
                        1:dof, 1:Qdof)
                else
                   phi => state%space%G_rule(Qnum)%phi(elem%type, je, 1, 1:dof, 1:Qdof)
                endif
             endif

             do k=1,ndim
                elem%w_limit(je, 1:Qdof, k) = matmul( ww(k, 1: elem%dof) , phi(1:dof, 1:Qdof) )
             enddo


             !do l=1,Qdof
             !   write(24, *) elem%xi(je, l, 1:2),  elem%w_limit(je, l, 1),  '  222y73h'
             !enddo

          enddo

          deallocate(wi, ww)
          !stop "83deu3d3h"
       end if

    enddo  ! i=1,grid%nelem


    !print*,'#  ComputeCapacityConductivity_2 called',state%nlSolver%implicitly
    !stop
  end subroutine ComputeCapacityConductivity_2


  !>  evaluation of the characteristic parameters of the porous media flow
  subroutine Eval_Characteristic_paramets ( gridA )
    use  porous_data_module
    type(mesh), intent(inout) :: gridA
    class(element), pointer :: elem
    class(Time_rule), pointer :: T_rule
    real, dimension(:), allocatable :: vals
    integer :: i, l, ie, k, alpha, Qdeg, Qdof, dof
    real :: val, flux, rel_diff, r0, r1, r2, r3, DBC, val0
    real :: wce_R, wce_L, wce_P, WC1, WC2, WC_lost
    real :: bilance_loc, bilance_glob, xii(2), t2
    integer :: ifile
    character*10 :: name

    real :: test

    test = 0.

    ifile = 10
    open(ifile, file = 'PM_flow', status='UNKNOWN', position='APPEND')

    !open(93, file = 'pm_data_IN', status='UNKNOWN', position='APPEND')
    !write(93,'(x)')
    
    Qdeg = state%time%Qnum

    T_rule => state%time%T_rule(Qdeg)

    ! verification of the emission BC, output in 'PM_emiss' file
    !call SetEmissionBC( .true. )

    
    associate ( time => state%time)
      select type ( time )
      class is ( TimeTDG_t )

         ! k = 10000
         ! xii(1:2) = 0.
         ! do i=0,k
         !    val = -50 + 40. *i /k
         !    call Eval_water_content(val, val0, grid%elem(1)%xi(0, 1, 2+1:2+iRe-1), xii(1:2))
         !    !write(34, *) val, val0
         ! enddo
         ! !stop "d83d3hj3"


         ! water content at t=0

         !print*,'porous_evals%water_content_init  = ', porous_evals%water_content_init  ,state%time%iter

         ! initialization at t= 0
         if( (state%time%ttime - state%time%tau(1)) /state%time%tau(1)  <= 1E-3) then
            porous_evals%flow_inlet_Time = 0.
            porous_evals%flow_bound_Time = 0.
            porous_evals%DBC_inaccuracy_Time = 0.

            porous_evals%water_content_init = 0.

            porous_evals%water_content_losses = 0.
            porous_evals%algeb_resid_tot = 0.

            do i = 1, gridA%nelem
               elem => gridA%elem(i)
               dof =elem%dof
               do k=1,ndim
                  elem%w(0,(k-1)*dof+1:dof) =  elem%wSTfin(k,1:dof)
                  elem%wActual(k, 1:dof) =  elem%wSTfin(k,1:dof)
               enddo

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

               call Elem_water_content(elem, val)
               porous_evals%water_content_init  = porous_evals%water_content_init  + val

               !write(*,'(a8, i5, 30es14.6)') 'WC in:',i, val, elem%w(0,1:dof)

               !if(i <= -5 .or. i== grid%nelem)  write(*,'(a10, i5, 30es14.6)') &
               !     ' WC0 : ', i, val, porous_evals%water_content_init
               !write(*,'(a8, i5, 30es14.6)') 'WC in:',i, val
               !if(i == 1) val0 = val
            enddo
            !print*,'porous_evals%water_content_init  = ', porous_evals%water_content_init , &
            !     state%time%ttime , state%time%tau(1), state%time%iter

            ! storing of the value at t_{m-1} from the left
            porous_evals%water_content_old = porous_evals%water_content_init

            write(ifile, 998) 0, 0.,  porous_evals%water_content_init,  porous_evals%water_content_init,&
                 0., 0., 0., 0.,0., 0., 0., 0.,0., 0., 0., 0., 0., 0., 0, 0, 0., 0., 0., 0., 0.

            !write(*, *) 'PM ___ INIT',  porous_evals%water_content_init,  porous_evals%water_content_init

         endif

         ! water content at actual time t_m^- and t_{m-1}^+
         porous_evals%water_content  = 0.
         porous_evals%water_content_back  = 0.

         WC1 = 0.
         WC2 = 0.
         WC_lost = 0.

         do i = 1, gridA%nelem
            elem => gridA%elem(i)
            dof =elem%dof

            if (Qdeg /= elem%TQnum) then
               !F@R Verify if it is OK, some nodes could be in wrong position
               stop 'Verify if it is OK, some nodes could be in wrong position'
            endif

            ! water contents at t_{m-1}^-
            do k=1,ndim
               elem%w(0,(k-1)*dof+1:dof) =  elem%wSTfin(k,1:dof)
               elem%wActual(k, 1:dof) =  elem%wSTfin(k,1:dof)
            enddo

            call Elem_water_content(elem, wce_p)


            ! solution w at t_{m-1}^+
            call Transfer_wST_to_w_Elem(elem , -1, Qdeg)
            call Transfer_wST_to_wActual_Elem(elem , -1, Qdeg)

            call Elem_water_content(elem, wce_L)
            porous_evals%water_content_back  = porous_evals%water_content_back  + wce_L


            ! solution w at t_m^-
            call Transfer_wST_to_w_Elem(elem , 0, Qdeg)
            call Transfer_wST_to_wActual_Elem(elem , 0, Qdeg)

            call Elem_water_content(elem, wce_R)
            porous_evals%water_content  = porous_evals%water_content  + wcE_R
            !write(*,'(a8, i5, 30es14.6)') 'WC m-:',i, wcE_R,porous_evals%water_content 

            !write(*,'(a8, i5, 30es14.6)') 'WC m-:',i, val, elem%w(0,1:dof)
            WC1 = WC1 + (wce_R - wce_p)
            WC2 = WC2 + (wce_R - wce_L)
            WC_lost = WC_lost + (wce_L - wce_p)

            test = test + elem%rhsST(1, 1, 1)
         enddo  ! i=1,gridA%nelem

         !!write(*, *) 'PM &&& INIT',  porous_evals%water_content,  porous_evals%water_content_back

         porous_evals%algeb_resid = test
         porous_evals%algeb_resid_tot = porous_evals%algeb_resid_tot + test
         !!print*,'test:;:', test,  elem%rhsST(1, 1, 1)

         ! difference between water contents at t_{m-1}^- and t_{m-1}^+
         r3 = porous_evals%water_content_back - porous_evals%water_content_old

         ! accumulation of the losses
         porous_evals%water_content_losses = porous_evals%water_content_losses + r3

         !write(*,'(a30,2i8, 2es12.4, a6, es12.4)') &
         !     '  porous_evals%water_content  = ', gridA%nelem, state%nsize,&
         !     porous_evals%water_content, porous_evals%water_content - porous_evals%water_content_init,&
         !     ',  t = ', state%time%ttime

         ! write(*,'(a30, 20es18.10)') '%water_content_init  = ', porous_evals%water_content_init
         ! write(*,'(a30, 20es18.10)') '%water_content_old  = ', porous_evals%water_content_old
         ! write(*,'(a30, 20es18.10)') '%water_content_back  = ', porous_evals%water_content_back, &
         !      porous_evals%water_content_back -  porous_evals%water_content_init, wC2, wC_lost
         ! write(*,'(a30, 20es18.10)') '%water_content  = ', porous_evals%water_content, &
         !      porous_evals%water_content - porous_evals%water_content_init, &
         !      porous_evals%water_content - porous_evals%water_content_old, WC1
         ! print*

         !flow throught the boundary
         porous_evals%flow_inlet = 0.
         porous_evals%flow_bound = 0.

         ! violation of the Dirichlet BC
         porous_evals%DBC_inaccuracy = 0

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

            do ie=1, elem%flen
               k = elem%face(neigh, ie)

               if(k <= 0) then  ! boundary edge

                  flux = 0.
                  DBC = 0.

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

                     ! save the wST space-time solution in quadrature index alpha to w
                     state%time%ctime = state%time%ttime - state%time%tau(1)*  T_rule%weights(alpha)
                     call Transfer_wST_to_w_Elem(elem , alpha, Qdeg)
                     call Transfer_wST_to_wActual_Elem(elem , alpha, Qdeg)

                     ! flux though edge at time t
                     call Elem_Flux_edge(elem, ie, val, alpha)

                     ! ST flux
                     flux = flux + val * T_rule%weights(alpha) * state%time%tau(1)
                     !if(i == 391) then
                     !   write(*,'(a10, 2i5, 5es12.4)') &
                     !     'times:',i,alpha, state%time%ctime, state%time%tau(1), val, flux
                     !call PlotElemFunction3D(10, elem, elem%dof, elem%w(0, 1:elem%dof) )
                     !endif

                     !write(*,'(a8, 3i5, 30es12.4)') 'Flux$$.:',i, ie,  alpha, val, flux

                     ! Dirichlet BC
                     if( elem%iBC(ie) /= 0 .and. elem%tBC(ie) /= 3) then
                        call Elem_DBC_jump(elem, ie, val)
                        DBC = DBC + val * T_rule%weights(alpha) * state%time%tau(1)

                        !write(*,'(a8, i5, 30es12.4)') ' $$.:', alpha, T_rule%weights(alpha), val, DBC

                     endif
                     !write(*,'(a10,6i5, 30es12.4)') 'flux::;',alpha, i,ie, k, &
                     !     elem%ibc(ie), elem%tbc(ie), elem%xc(:), &
                     !     val, flux
                  enddo  ! alpha

                  !if(elem%i >= 1 .and. elem%i <= 8) &
                  !     write(*,'(a10,6i5, 30es18.10)') 'flux::;', i,ie, k, elem%face(nei_i, ie), &
                  !     elem%ibc(ie), elem%tbc(ie), flux
                  !     val, flux


                  ! all boundary
                  porous_evals%flow_bound =  porous_evals%flow_bound + flux

                  ! only inlet
                  if( elem%iBC(ie) /= 0 .and. elem%tBC(ie) /= -3) then
                     porous_evals%flow_inlet = porous_evals%flow_inlet + flux

                     porous_evals%DBC_inaccuracy =  porous_evals%DBC_inaccuracy + DBC

                     !write(93,'(a8, es12.4, 4i5, 30es12.4)') ' DBC:', state%time%ttime, &
                     !     state%time%iter, &
                     !     elem%i, elem%iBC(ie),elem%tBC(ie), elem%xc(:), flux,porous_evals%flow_inlet
                     !print*
                     !write(93, *) elem%xc
                  endif

               endif ! if( k <= 0)

            enddo ! ie = 1, elem%flen

         enddo  ! i=1,gridA%nelem
         !close(93)

         ! total flux over (0, t_actual)
         porous_evals%flow_inlet_Time = porous_evals%flow_inlet_Time + porous_evals%flow_inlet

         porous_evals%flow_bound_Time = porous_evals%flow_bound_Time + porous_evals%flow_bound

         porous_evals%DBC_inaccuracy_Time = porous_evals%DBC_inaccuracy_Time + porous_evals%DBC_inaccuracy

         ! write(*,'(a30, 20es18.10)') '%flow_bound  = ', porous_evals%flow_bound, porous_evals%flow_bound_Time
         ! write(*,'(a30, 20es18.10)') '%flow_inlet  = ', porous_evals%flow_inlet, porous_evals%flow_inlet_Time
         ! write(*,'(a30, 20es18.10)') '%flow_NC  = ',porous_evals%flow_bound - porous_evals%flow_inlet, &
         !      porous_evals%flow_bound_Time - porous_evals%flow_inlet_Time
         ! write(*,'(a30, 20es18.10)') '%DBC  = ', porous_evals%DBC_inaccuracy , &
         !       porous_evals%DBC_inaccuracy_Time
         ! write(*,'(a30, 20es18.10)') 'rezid = %DBC-%flow_inlet  = ', &
         !      porous_evals%flow_inlet - porous_evals%DBC_inaccuracy, &
         !      porous_evals%flow_inlet_Time - porous_evals%DBC_inaccuracy_Time
         ! print*,'_________________________________________________'


         ! difference between water content and the flux
         !print*,'WWCDD:', porous_evals%water_content, porous_evals%water_content_init
         r1 = porous_evals%water_content - porous_evals%water_content_init
         r2 = porous_evals%flow_bound_Time
         r0 = porous_evals%flow_inlet_Time
         rel_diff = abs ( r1 - r2) / max(1E-15, max(abs(r1), abs(r2) ) )

         ! write(*,*)
         ! write(*,'(a6, 30es12.4)') 'PM_flo:', &
         !      porous_evals%water_content_init, &
         !      porous_evals%water_content_old, &
         !      porous_evals%water_content_back, &
         !      porous_evals%water_content, &
         !      porous_evals%flow_inlet, porous_evals%flow_bound
         ! write(*,*)


         !bilance among water contents, flow and violation of DBC
         bilance_loc = porous_evals%water_content -porous_evals%water_content_old &
              - porous_evals%flow_bound + porous_evals%DBC_inaccuracy
         bilance_glob = r1 - r2 + porous_evals%DBC_inaccuracy_Time

         if(mod(state%time%iter, 10) == 0) &
              write(*,'(a10, i5, i7, a4, es8.2, a2,2es10.2,a1,es10.2,a1,es9.2,a2,2es10.2)') &
              'PM chars:', gridA%nelem, state%nsize ,&
              ', t=', state%time%ttime, '|', r1, r2, '(',r0,')',rel_diff,'|', &
              bilance_loc, bilance_glob
         !r3,  porous_evals%water_content_losses/ r1

         ! print*,'__________________________________________________________________'
         ! print*,' porous_evals%water_content_ini =', porous_evals%water_content_init
         ! print*,' porous_evals%water_content_old =', porous_evals%water_content_old
         ! print*,' porous_evals%water_content_bak =', porous_evals%water_content_back, r3
         ! print*,' porous_evals%water_content_act =', porous_evals%water_content
         ! print*,' flow difference                =', r1
         ! print*
         ! print*,' porous_evals%flow_bound = ', porous_evals%flow_bound, r2
         ! print*,' BC _ violance           = ',  porous_evals%DBC_inaccuracy
         ! print*,' bilance_loc             = ', porous_evals%water_content -porous_evals%water_content_old&
         !      - porous_evals%flow_bound + porous_evals%DBC_inaccuracy
         ! print*,' bilance                 = ', r1 - r2 + porous_evals%DBC_inaccuracy_Time
         ! print*,'__________________________________________________________________'
         ! write(87,'(a12, 2es12.4)') 'bilances:', &
         !      porous_evals%water_content -porous_evals%water_content_old&
         !      - porous_evals%flow_bound + porous_evals%DBC_inaccuracy, &
         !      r1 - r2 + porous_evals%DBC_inaccuracy_Time

         !write(*,'(a30,20es12.4)') '  porous_evals%flows:  = ', &
         !     porous_evals%flow_inlet, porous_evals%flow_bound, &
         !     porous_evals%flow_inlet_Time, porous_evals%flow_bound_Time

         ! file 'PM_flow' outut
         write(ifile, 998) state%time%iter, state%time%ttime, &                        ! 1..2
              porous_evals%water_content_init,  porous_evals%water_content, &          ! 3..4
              r1, r2, rel_diff, abs(r3), abs(porous_evals%water_content_losses), &     ! 5..9
              porous_evals%water_content_losses / max(1E-15, r1), &                    ! 10
              porous_evals%flow_inlet, porous_evals%flow_bound, &                      ! 11..12
              porous_evals%flow_inlet_Time, porous_evals%flow_bound_Time, &            ! 13..14
              abs(porous_evals%DBC_inaccuracy),  abs(porous_evals%DBC_inaccuracy_Time), & ! 15..16
              bilance_loc, bilance_glob,  gridA%nelem, state%nsize, &                  ! 17..20
              abs(porous_evals%algeb_resid),  abs(porous_evals%algeb_resid_tot), &        ! 21..22
              porous_evals%water_content - porous_evals%water_content_old, &   ! 23
              abs(  porous_evals%water_content - porous_evals%water_content_old - porous_evals%flow_bound)/ &
              max(1E-25, abs(porous_evals%flow_bound) ), &                        ! 24
              porous_evals%flow_inlet/ state%time%tau(1)    ! 25
         close(ifile)

         !open(99, file = 'fort.99', status = 'unknown', position='APPEND')
         !call cpu_time(t2)
         !write(99, 999) grid%nelem, state%nsize , state%space%deg, state%time%deg, &
         !     state%time%iter, state%NlSolver%Aiter, state%linSolver%iter_tot, state%time%ttime, &
         !     state%time%ttime/ state%time%iter, &
         !     r1, r2, rel_diff, &
         !     sqrt(state%estim(resA:resST, 1)),   t2- state%start_time
         !close(99)
         !stop "83y38ih3"

         ! storing the value for the next time step
         porous_evals%water_content_old  = porous_evals%water_content

      class default
         stop 'Eval_Characteristic_paramets, For TimeTDG_t only!!!'
      end select
    end associate ! time
!998 format( i5, 17es12.4, i6, i8, 4es12.4 )
998 format( i5, 3es12.4, 2es16.8, 12es12.4, i6, i8, 5es12.4 )
999 format( 6i7, i10, 2es12.4, 2es16.8, 12es12.4)

  end subroutine Eval_Characteristic_paramets

  !> compute the water content on one element
  subroutine Elem_water_content(elem, val)
    class(element), intent(inout) :: elem
    real, intent(inout) :: val
    real, dimension(:,:), allocatable :: wi
    real, dimension(:), allocatable :: theta
    integer :: Qdof, l

    Qdof = elem%Qdof
    allocate(wi(1:Qdof, 1:ndim), theta(1:Qdof) )

    ! w in integ nodes
    call Eval_w_Elem(elem, wi(1:Qdof, 1:ndim) )

    ! evaluation of the water content in integ nodes
    do l=1,Qdof
       call Eval_water_content(wi(l,1), theta(l), elem%xi(0, l, 2+1:2+iRe-1), &
            elem%xi(0, l,  1:nbDim) )
       !write(33, *) elem%xi(0, l,  1:nbDim), wi(l,1), wi(l,1)-elem%xi(0,l,2), theta(l), '  3y4333'
       !write(*,'(a8, 2i5, 30es12.4)') 'wc_i:',elem%i, l, wi(l, 1), theta(l)
       !if(elem%i == 1 .and. l<=3) write(*,'(a8, i5, 30e12.4)') 'wcc:',elem%i,  wi(l,1), theta(l)

    enddo

    ! integration of the water content over elem
    call IntegrateFunction(elem, theta(1:Qdof),  val)
    !if(elem%i == 1 ) write(*,'(a8, i5, 300es12.4)') 'integ;',elem%i, val

    deallocate(wi, theta)
  end subroutine Elem_water_content



  !> compute the violation of the Dirichlet boundary condition
  subroutine Elem_DBC_jump(elem, ie, val)
    class(element), intent(inout) :: elem
    integer, intent(in) :: ie  ! index of the edge
    real, intent(inout) :: val
    real, dimension(:,:), allocatable :: wi, wB, Re_1
    real, dimension(:), allocatable :: penal, penal2, weights
    integer :: Qdof, Qnum, l

    Qdof = elem%face(fGdof,ie)
    Qnum = elem%face(fGnum,ie)

    allocate(wi(1:Qdof, 1:ndim), wB(1:Qdof, 1:nDim) )
    allocate(penal(1:Qdof), weights(1:Qdof))

    ! w in integ nodes
    call Eval_w_Edge(elem, ie, wi(1:Qdof, 1:ndim), .false.) !!elem%n(ie,:) !

    ! w_B in integ nodes
    call Exact_Porous(Qdof, grid%b_edge(-elem%face(neigh,ie))%x_div(1:Qdof, 1:nbDim), &
         wB(1:Qdof,1:ndim), state%time%ctime )


    weights(1:Qdof) = state%space%G_rule(Qnum)%weights(1:Qdof) !* elem%dn(ie)
    !penal(1:Qdof) = 1. / state%model%Re  ! elem%d_gamma   ! hiden in the weights


    ! FR comment the following, setPenalty_new is used instead


    allocate(Re_1(1:iRe,1:Qdof) )
    if(state%model%Re > 0.) then
       Re_1(1,1:Qdof) = 1./state%model%Re
    endif
    ! ! set the penalty OLD  -- DOES NOT WORK for discontinuous difusion !!!!!!
    ! allocate(penal2(1:Qdof) , source = 0.0)
    ! call setPenaltySigma( elem, ie, Qdof, Re_1(1:iRe, 1:Qdof), elem%d_gamma, penal2(1:Qdof) )
    ! penal(1:Qdof) =  penal2(1:Qdof)

    ! testing new subroutine
    call setPenaltySigma_New( elem, ie, Qdof, penal(1:Qdof) )
    !if ( norm2( penal - penal2 ) > 1.E-14 ) then
    !!  print*, 'New version of setPenalty differs from the old one2!', norm2( penal - penal2 )
    !!  stop
    !endif
    !deallocate(penal2)
    deallocate(Re_1)


    val = 0.
    ! integration of the penalty over the edge
    do l=1,Qdof

       val = val +  weights(l) *  penal(l) *  ( wi(l,1) - wB(l,1) )

    enddo

    !write(*,'(a8, i5, 30es12.4)') ' ed wi:', elem%i, wi(:, 1)
    !write(*,'(a8, i5, 30es12.4)') ' ed wB:', Qdof, wB(:, 1)
    !write(*,'(a8, i5, 30es12.4)') ' diff:', -99, wi(:, 1)-wB(:, 1)
    !print*,'----------------------', val

    !print*,'val =' ,val
    !
    !call  IntegrateFunctionNormalEdge(elem, ie, transpose(Dwi(1:Qdof, 1, 1:nbdim)), K_sk(1:Qdof, 1, 1),&
    !     val)
    !print*,'val =' ,val
    !print*

    deallocate(wi, wB, penal, weights)

  end subroutine Elem_DBC_jump


  !> setting of the quantity "water content" in order to use  for the adaptation
  !> solution is stored in elem%wS(1:ndimL, 1:elem%dof), this array is replaced
  subroutine SetWaterContent_for_Metric(ndimL, elem)
    integer, intent(in) :: ndimL  ! index of the edge
    class(element), intent(inout) :: elem
    real, dimension(:), allocatable :: wi
    real, dimension(:), allocatable :: theta
    real, dimension(:,:), allocatable :: ww
    type(volume_rule), pointer :: V_rule
    integer :: Qdof, dof, l, n

    dof = elem%dof
    Qdof = elem%Qdof
    allocate(wi(1:Qdof), theta(1:Qdof) )

    V_rule => state%space%V_rule(elem%Qnum)

    ! array for the DG reconstruction
    allocate(ww(1:dof, 1:ndimL), source = 0.0 )

    do n = 1, ndimL
       ! w in integ nodes
       call Eval_wS_Elem(elem, elem%wS(n, 1:dof), V_rule, wi(1:Qdof) )

       ! evaluation of the water content in integ nodes
       do l=1,Qdof
          call Eval_water_content(wi(l), theta(l), elem%xi(0, l, 2+1:2+iRe-1), &
               elem%xi(0, l,  1:nbDim) )
          !write(30+n, *) elem%xi(0, l,  1:nbDim), wi(l), theta(l), '  3y4333'
       enddo

       ! computation of the RHS of the projection
       call IntegrateVectorB(elem, dof, theta(1:Qdof), ww(1:dof, n) )

    enddo

    ! solution of the projection
    call SolveLocalMatrixProblem(dof, elem%mass%Mb(1:dof, 1:dof), ndimL, ww(1:dof, 1:ndimL))

    ! assembling the array back
    elem%wS(1:ndimL,  1:dof) = transpose( ww(1:dof, 1:ndimL) )

    ! graphical output
    !do n = 1, ndimL
    !   call PlotElemFunction3D(30+n, elem, dof, elem%wS(n, 1: dof) )
    !end do
    !if(elem%i == grid%nelem) stop "73y33gg"

  end subroutine SetWaterContent_for_Metric



  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! FORCHHEIMER  MODEL
  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
  
    !>  evaluation of the characteristic parameters of the porous media flow
  subroutine Eval_Forchheimer_paramets ( gridA )
    use  porous_data_module
    type(mesh), intent(inout) :: gridA
    class(element), pointer :: elem
    class(Time_rule), pointer :: T_rule
    real, dimension(:), allocatable :: vals
    integer :: i, l, ie, k, alpha, Qdeg, Qdof, dof
    real :: val, flux, rel_diff, r0, r1, r2, r3, DBC, val0
    real :: wce_R, wce_L, wce_P, WC1, WC2, WC_lost
    real :: bilance_loc, bilance_glob, xii(2)
    integer :: ifile

    real :: test

    test = 0.

    ifile = 10
    open(ifile, file = 'PM_flow', status='UNKNOWN', position='APPEND')

    Qdeg = state%time%Qnum

    T_rule => state%time%T_rule(Qdeg)

    ! verification of the emission BC, output in 'PM_emiss' file
    !call SetEmissionBC( .true. )

    
    associate ( time => state%time)
      select type ( time )
      class is ( TimeTDG_t )

         ! k = 10000
         ! xii(1:2) = 0.
         ! do i=0,k
         !    val = -50 + 40. *i /k
         !    call Eval_water_content(val, val0, grid%elem(1)%xi(0, 1, 2+1:2+iRe-1), xii(1:2))
         !    !write(34, *) val, val0
         ! enddo
         ! !stop "d83d3hj3"


         ! water content at t=0

         !print*,'porous_evals%water_content_init  = ', porous_evals%water_content_init  ,state%time%iter

         ! initialization at t= 0
         if( (state%time%ttime - state%time%tau(1)) /state%time%tau(1)  <= 1E-3) then
            porous_evals%flow_inlet_Time = 0.
            porous_evals%flow_bound_Time = 0.
            porous_evals%DBC_inaccuracy_Time = 0.

            porous_evals%water_content_init = 0.

            porous_evals%water_content_losses = 0.
            porous_evals%algeb_resid_tot = 0.

            do i = 1, gridA%nelem
               elem => gridA%elem(i)
               dof =elem%dof
               do k=1,ndim
                  elem%w(0,(k-1)*dof+1:dof) =  elem%wSTfin(k,1:dof)
                  elem%wActual(k, 1:dof) =  elem%wSTfin(k,1:dof)
               enddo

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

               call Elem_water_content_Forch(elem, val)
               porous_evals%water_content_init  = porous_evals%water_content_init  + val

               !write(*,'(a8, i5, 30es14.6)') 'WC in:',i, val, elem%w(0,1:dof)

               !if(i <= -5 .or. i== grid%nelem)  write(*,'(a10, i5, 30es14.6)') &
               !     ' WC0 : ', i, val, porous_evals%water_content_init
               !write(*,'(a8, i5, 30es14.6)') 'WC in:',i, val
               !if(i == 1) val0 = val
            enddo
            !print*,'porous_evals%water_content_init  = ', porous_evals%water_content_init , &
            !     state%time%ttime , state%time%tau(1), state%time%iter

            ! storing of the value at t_{m-1} from the left
            porous_evals%water_content_old = porous_evals%water_content_init

            write(ifile, 998) 0, 0.,  porous_evals%water_content_init,  porous_evals%water_content_init,&
                 0., 0., 0., 0.

            !write(*, *) 'PM ___ INIT',  porous_evals%water_content_init,  porous_evals%water_content_init

         endif

         ! water content at actual time t_m^- and t_{m-1}^+
         porous_evals%water_content  = 0.
         porous_evals%water_content_back  = 0.

         WC1 = 0.
         WC2 = 0.
         WC_lost = 0.

         do i = 1, gridA%nelem
            elem => gridA%elem(i)
            dof =elem%dof

            if (Qdeg /= elem%TQnum) then
               !F@R Verify if it is OK, some nodes could be in wrong position
               stop 'Verify if it is OK, some nodes could be in wrong position'
            endif

            ! water contents at t_{m-1}^-
            do k=1,ndim
               elem%w(0,(k-1)*dof+1:dof) =  elem%wSTfin(k,1:dof)
               elem%wActual(k, 1:dof) =  elem%wSTfin(k,1:dof)
            enddo

            call Elem_water_content_Forch(elem, wce_p)


            ! solution w at t_{m-1}^+
            call Transfer_wST_to_w_Elem(elem , -1, Qdeg)
            call Transfer_wST_to_wActual_Elem(elem , -1, Qdeg)

            call Elem_water_content_Forch(elem, wce_L)
            porous_evals%water_content_back  = porous_evals%water_content_back  + wce_L


            ! solution w at t_m^-
            call Transfer_wST_to_w_Elem(elem , 0, Qdeg)
            call Transfer_wST_to_wActual_Elem(elem , 0, Qdeg)

            call Elem_water_content_Forch(elem, wce_R)
            porous_evals%water_content  = porous_evals%water_content  + wcE_R
            !write(*,'(a8, i5, 30es14.6)') 'WC m-:',i, wcE_R,porous_evals%water_content 

            !write(*,'(a8, i5, 30es14.6)') 'WC m-:',i, val, elem%w(0,1:dof)
            WC1 = WC1 + (wce_R - wce_p)
            WC2 = WC2 + (wce_R - wce_L)
            WC_lost = WC_lost + (wce_L - wce_p)

            test = test + elem%rhsST(1, 1, 1)
         enddo  ! i=1,gridA%nelem

         !!write(*, *) 'PM &&& INIT',  porous_evals%water_content,  porous_evals%water_content_back

         porous_evals%algeb_resid = test
         porous_evals%algeb_resid_tot = porous_evals%algeb_resid_tot + test
         !!print*,'test:;:', test,  elem%rhsST(1, 1, 1)

         ! difference between water contents at t_{m-1}^- and t_{m-1}^+
         r3 = porous_evals%water_content_back - porous_evals%water_content_old

         ! accumulation of the losses
         porous_evals%water_content_losses = porous_evals%water_content_losses + r3

         !write(*,'(a30,2i8, 2es12.4, a6, es12.4)') &
         !     '  porous_evals%water_content  = ', gridA%nelem, state%nsize,&
         !     porous_evals%water_content, porous_evals%water_content - porous_evals%water_content_init,&
         !     ',  t = ', state%time%ttime

         ! write(*,'(a30, 20es18.10)') '%water_content_init  = ', porous_evals%water_content_init
         ! write(*,'(a30, 20es18.10)') '%water_content_old  = ', porous_evals%water_content_old
         ! write(*,'(a30, 20es18.10)') '%water_content_back  = ', porous_evals%water_content_back, &
         !      porous_evals%water_content_back -  porous_evals%water_content_init, wC2, wC_lost
         ! write(*,'(a30, 20es18.10)') '%water_content  = ', porous_evals%water_content, &
         !      porous_evals%water_content - porous_evals%water_content_init, &
         !      porous_evals%water_content - porous_evals%water_content_old, WC1
         ! print*

         ! difference between water content and the flux
         !print*,'WWCDD:', porous_evals%water_content, porous_evals%water_content_init
         r1 = porous_evals%water_content - porous_evals%water_content_init
         r2 = porous_evals%flow_bound_Time
         r0 = porous_evals%flow_inlet_Time
         rel_diff = abs ( r1 - r2) / max(1E-15, max(abs(r1), abs(r2) ) )

         ! write(*,*)
         ! write(*,'(a6, 30es12.4)') 'PM_flo:', &
         !      porous_evals%water_content_init, &
         !      porous_evals%water_content_old, &
         !      porous_evals%water_content_back, &
         !      porous_evals%water_content, &
         !      porous_evals%flow_inlet, porous_evals%flow_bound
         ! write(*,*)


         !bilance among water contents, flow and violation of DBC
         bilance_loc = porous_evals%water_content -porous_evals%water_content_old &
              - porous_evals%flow_bound + porous_evals%DBC_inaccuracy
         bilance_glob = r1 - r2 + porous_evals%DBC_inaccuracy_Time

         write(*,'(a10, i5, i7, a4, es8.2, a2,2es10.2,a1,es10.2,a1,es9.2,a2,2es10.2)') &
              'PM chars:', gridA%nelem, state%nsize ,&
              ', t=', state%time%ttime, '|', r1, r2, '(',r0,')',rel_diff,'|', &
              bilance_loc, bilance_glob
              !r3,  porous_evals%water_content_losses/ r1

         ! print*,'__________________________________________________________________'
         ! print*,' porous_evals%water_content_ini =', porous_evals%water_content_init
         ! print*,' porous_evals%water_content_old =', porous_evals%water_content_old
         ! print*,' porous_evals%water_content_bak =', porous_evals%water_content_back, r3
         ! print*,' porous_evals%water_content_act =', porous_evals%water_content
         ! print*,' flow difference                =', r1
         ! print*
         ! print*,' porous_evals%flow_bound = ', porous_evals%flow_bound, r2
         ! print*,' BC _ violance           = ',  porous_evals%DBC_inaccuracy
         ! print*,' bilance_loc             = ', porous_evals%water_content -porous_evals%water_content_old&
         !      - porous_evals%flow_bound + porous_evals%DBC_inaccuracy
         ! print*,' bilance                 = ', r1 - r2 + porous_evals%DBC_inaccuracy_Time
         ! print*,'__________________________________________________________________'
         ! write(87,'(a12, 2es12.4)') 'bilances:', &
         !      porous_evals%water_content -porous_evals%water_content_old&
         !      - porous_evals%flow_bound + porous_evals%DBC_inaccuracy, &
         !      r1 - r2 + porous_evals%DBC_inaccuracy_Time

         !write(*,'(a30,20es12.4)') '  porous_evals%flows:  = ', &
         !     porous_evals%flow_inlet, porous_evals%flow_bound, &
         !     porous_evals%flow_inlet_Time, porous_evals%flow_bound_Time

         write(ifile, 998) state%time%iter, state%time%ttime, &                        ! 1..2
              porous_evals%water_content_init,  porous_evals%water_content, &          ! 3..4
              r1, r2, rel_diff, abs(r3), abs(porous_evals%water_content_losses), &     ! 5..9
              porous_evals%water_content_losses / max(1E-15, r1), &                    ! 10
              porous_evals%flow_inlet, porous_evals%flow_bound, &                      ! 11..12
              porous_evals%flow_inlet_Time, porous_evals%flow_bound_Time, &            ! 13..14
              abs(porous_evals%DBC_inaccuracy),  abs(porous_evals%DBC_inaccuracy_Time), & ! 15..16
              bilance_loc, bilance_glob,  gridA%nelem, state%nsize, &                  ! 17..20
              abs(porous_evals%algeb_resid),  abs(porous_evals%algeb_resid_tot), &        ! 21..22
              porous_evals%water_content - porous_evals%water_content_old, &   ! 23
              abs(  porous_evals%water_content - porous_evals%water_content_old - porous_evals%flow_bound)/ &
              max(1E-25, abs(porous_evals%flow_bound) )                        ! 24
         close(ifile)
         !stop "83y38ih3"

         ! storing the value for the next time step
         porous_evals%water_content_old  = porous_evals%water_content

      class default
         stop 'Eval_Characteristic_paramets, For TimeTDG_t only!!!'
      end select
    end associate ! time
998 format( i5, 17es12.4, i6, i8, 4es12.4 )

  end subroutine Eval_Forchheimer_paramets


  !> compute the water content on one element
  subroutine Elem_water_content_Forch(elem, val)
    class(element), intent(inout) :: elem
    real, intent(inout) :: val
    real, dimension(:,:), allocatable :: wi
    real, dimension(:), allocatable :: theta
    real :: compress 
    integer :: Qdof, l

    ! compressibility HAS to be in agreement with model.f90, compress = 5E-10
    compress = 5E-10
    
    Qdof = elem%Qdof
    allocate(wi(1:Qdof, 1:ndim), theta(1:Qdof) )

    ! w in integ nodes
    call Eval_w_Elem(elem, wi(1:Qdof, 1:ndim) )

    ! evaluation of the density content in integ nodes
    do l=1,Qdof
       theta(l) = exp(compress * wi(l,  1)  )

    enddo

    ! integration of the water content over elem
    call IntegrateFunction(elem, theta(1:Qdof),  val)
    !if(elem%i == 1 ) write(*,'(a8, i5, 300es12.4)') 'integ;',elem%i, val

    deallocate(wi, theta)
  end subroutine Elem_water_content_Forch

  
end module pm_fluxes
