!>  basis geometrical entities which cannot be stored in parallel mode
!> (should shared by each node)
module define_state


  use adaptation_mod
  use boundary_mod
  use cpu_time_mod
  use integration
  use linear_mod
  use model_mod
  use modelNS_mod
  use porous_mod
  use nonlinear_mod
  use loc_rav_tho_ned
  use rav_tho_nedelec
  use scalar_mod
  use space_mod
  use time_mod

  implicit none
  public

  !> shape of matrices generated by DGFEM
  type :: MatrixShape
     integer :: nsize                         ! matrix dimension
     integer :: nonzero                    !!!! number of all nonzero elements
     integer, dimension(:), pointer :: irwst  ! begins of rows
     integer, dimension(:), pointer :: idx    ! column of elements in rows
     integer, dimension(:,:,:), pointer ::nij   ! links to matrix element from triangle

     !integer, dimension(:,:), pointer ::ivc   ! begins of a vector sequence
     !integer, dimension(:,:), pointer ::nivc  ! see InitMatrixShape
  end type MatrixShape

  !> type of the scalar problem
  type :: scalar_type
     integer :: icase                        ! index of the scalar problem
     integer :: idiff                        ! index of diffusion
     integer :: iconv                        ! index of convection
     integer :: ireac                        ! index of reaction
     integer :: iexact                       ! index of the exact solution + RHS
     real :: conv_rat                        ! parameter for the extrapolation of Jacobi
     real :: param1                          ! first parameter
  end type scalar_type

  !> main data structure
  type, public :: solver_state
     type(rtn_fe), dimension(:), pointer :: RTN => NULL()            ! RTN elements
     type(basis_rtn_fe), dimension(:), pointer :: loc_RTN => NULL()  ! local RTN elements
     real, dimension(:), pointer :: A => NULL()                      ! matrix for conforming FEM
     ! logical :: implicitly                ! .true. = implicit performace ! OOP in nlSolver !NOT YET
     logical :: init_only_matrix_blocks   ! .true. = initilize only matrix blocks
     logical :: only_diag                 ! .true. = only diagonal blocks (for local problems)
     logical :: local_problem             ! .true. = solution of local problems
     logical :: RHS_presented             ! .true. = RHS is presented
     logical :: homogenDirichlet          ! .true. = a priori homogenous Dirichlet BC
     logical :: num_flux                  ! .true. = numerical flux, .false. = physical
     logical :: time_dependent            ! .true. = really time dependent problem , .false. otherwise

     !OOP should be in linSolver
     logical :: MGsolver
     logical :: MGLinSolver
     logical :: MGNlnSolver
     integer :: MaxMGlevel                ! maximal level for multigrid method
     integer :: CurMGlevel                ! current level for multigrid method
     integer, dimension(9) :: MGlvlSize
     integer, dimension(9) :: MGlvlDof
     integer, dimension(9) :: MGlvlDeg
     integer :: MGnsize                   ! vector size for MG

     integer :: type_IC                   ! type of IC
     integer :: itime                     ! = 0, no CPU time outputs

     integer :: isol, isol_save, imat         ! actual index of sol* files

     integer::  lin_solver_not_conv       ! linear solver did not converges

     real :: conv_rez                     ! stopping rezid for steady
     real :: EcD_tol                      ! stoppiing criterion for cD (drag)
     real :: EcL_tol                      ! stoppiing criterion for cL (lift)
     real :: EcM_tol                      ! stoppiing criterion for cM (moment)

     real :: ttime_save, timeprn_save     ! total time save for possible adaptation

     real, dimension(:,:), allocatable :: estim   ! total error estimator at current time level
     real, dimension(:,:), allocatable :: loc_estim   ! total error estimator
     real, dimension(:), allocatable :: T_estim   ! estim of the error for several time levels
     real, dimension(:), allocatable :: L_estim   ! estim of the error for 1 time level, since last successful control of the mesh

     real :: time_AD_start, time_AD_end   ! time levels for adaptations

     real :: estimTot_SC                  ! guaranteed upper bound on the energy error for u_h^i

     real, dimension(:), allocatable :: err   ! error
     real, dimension(:), allocatable :: errSTnorm   ! error in space-time norms
     real, dimension(:), allocatable :: errSTloc    ! error in space-time norms
     real, dimension(:,:), allocatable :: errSnorm  ! L2/H1,... norms in different time nodes 1-endpoint,2-a quadr. node, 3-1/2

     integer :: num_limits                ! number of limiting
     integer :: num_call_F                ! number of calling vector F
     integer :: num_call_C                ! number of calling matrix C

     real :: max_eigenvals                ! maximal eigenvals|Gi|/|Ki|

     integer ::  nsize                        ! # of unknowns = Mshape%nsize
     integer ::  nsizeP                       ! # of unknowns for pressure = Mshape%nsize
     integer ::  nonzero                      ! # of nonzero  elements
     integer ::  nonzeroVP                    ! # of nonzero  elements
     integer ::  nonzeroPV                    ! # of nonzero  elements

     integer :: numBC                         ! number of prescribed BC
     type(typeBC), dimension(:), allocatable:: BC ! prescribed boundary conditions !OOP used

     !TODO: FR what are the following used to? where tu put them?
     real :: rho_infty                        ! far field density
     real :: v_infty                          ! far field velocity
     real :: p_infty                          ! far field pressure
     real :: alpha_infty                      ! far field angle of attack
     real :: theta_infty                      ! far field temperature

     real, dimension(:, :), allocatable :: cDLM   !cDLM(k,1:3) = (cD, cL, cM, cDp, cLp) at t_k
     real, dimension(:), allocatable :: EcDLM            ! = cDLM_max - cDLM_min

     real :: ST_Vp, ST_Vc, ST_Ep, ST_Ec

     logical:: print

     real :: start_time                       ! initial CPU time
     real :: timeprn                          ! time for printing output
     real :: CPU_prepare, CPU_solve , CPU_estim  ! CPU time necessary for particular opeartions
     real :: CPU_estim2, CPU_adapt                      ! CPU time necessary for particular opeartions
     real :: CPU_constaint

     integer :: no_refused_steps              ! number of refused time steps

     !OOP doubled in mesh, where to put this? -> space%adapt
     !     integer, dimension(:,:), allocatable :: AMAhistory   !history of AMA computation
     logical :: tri_solA

     real, dimension(:,:), allocatable :: Schur_A

     real, dimension(:,:), allocatable :: cons  ! auxiliary array for verifying the code
     real :: LevelSet_Area


     integer, private :: p_mod ! actual polynomial degree with respect to space
     integer, private :: q_mod ! actual polynomial degree with respect to time

     integer :: p_mod_max ! max polynomial degree with respect to space
     integer :: q_mod_max ! max polynomial degree with respect to time

     integer, allocatable, dimension(:,:) :: bigNSize ! size of the big system (0:p_mod_max, 0:q_mod_max)
     integer, allocatable, dimension(:,:) :: bigNonzero ! size of the big system (0:p_mod_max, 0:q_mod_max)

     !!!!!ADGo DATA !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     character(len=20) :: modelName
     !  character(len=20) :: non_solver, lin_solver
     class(Model_t ) , allocatable :: model
     class(NonlinearSol_t ), allocatable :: nlSolver
     class(LinearSol_t), allocatable :: linSolver

     class(Space_t), allocatable :: space
     class(Time_t), allocatable :: time
     class(Cpu_t), allocatable :: cpuTime

     logical :: dwr_p ! DWR_P is used - to allocation of wSTplus, zSTplus and computation of bigger problems

     logical :: wActual != .true. !.true. ! TESTING for testing the new array wActual which should be used insted of the elem%w(0,


   contains

     procedure :: init_Problem ! new version of initProblem
     procedure :: initRTN
     procedure :: allocTimeDiscMethod
     procedure :: initErrorEstim
     procedure :: printInitialConditions
     procedure :: printSpaceErrorEstims
     procedure :: readModelData
     procedure :: readNonlinearSolData
     procedure :: readLinearSolData
     procedure :: setStateParameters
     procedure :: getBigNSize
     procedure :: getP_mod
     procedure :: getQ_mod
     procedure :: setP_mod
     procedure :: setQ_mod



  end type solver_state


contains


  subroutine readModelData ( this, Re, isca, t2)
    class( solver_state ), intent( inout ) :: this
    real, intent(in) :: Re
    integer, intent(in) :: isca
    real, intent( in ) :: t2

    associate ( problem => this%model )

    if ( this%modelName == 'sca' ) this%modelName = 'scalar'

    select case(this%modelName)
    case ( 'scalar' )

       ndim = 1
       allocate ( Scalar_t::this%model )

       select type ( problem )
       type is (Scalar_t)
          !FERROR what is isca used for?
          !   problem%isca = 5
          call problem%init( Re, isca,t2 )
       end select

    case ( 'NSe' )
       ndim = 4
       allocate( NavierStokes_t :: this%model)

       select type ( problem )
       type is ( NavierStokes_t )
          call problem%init( Re )
       end select

    case ( 'Eul' )
       ndim = 4
       allocate( Eulerian_t :: this%model)

       select type ( problem )
       type is ( Eulerian_t )
          call problem%init( Re, isca,t2)
       end select

    case ( 'RANS_2e' )
       ndim = 6
       allocate( RANS_2e_t :: this%model)

       select type ( problem )
       type is ( RANS_2e_t )
          call problem%init( Re )
       end select

    case ( 'pedes' )
       ndim = 3
       allocate( Pedestrian_t :: this%model)

       select type ( problem )
       type is ( Pedestrian_t )
          call problem%init( Re, isca, t2 )
       end select

    case ( 'swe' )
       ndim = 3
       allocate( ShallowWater_t :: this%model)

       select type ( problem )
       type is ( ShallowWater_t )
          call problem%init( Re, isca, t2 )
       end select

    case ( '2eqs')
       ndim = 2
       allocate ( TwoEqs_t::this%model )
       !print*, ' Model ',this%modelName,' is not implemented YET'
       select type ( problem )
       type is (TwoEqs_t)
          !FERROR what is isca used for?
          !   problem%isca = 5
          call problem%init( Re, isca,t2 )
       end select

    case ( 'porous')
       ndim = 1
       allocate ( Porous_t::this%model )
       !print*, ' Model ',this%modelName,' is not implemented YET'
       select type ( problem )
       type is (Porous_t)
          !FERROR what is isca used for?
          !   problem%isca = 5
          call problem%init( Re, isca,t2 )
       end select

    case ( 'wet_steam' )
       print*, ' Model ',this%modelName,' is not implemented YET'
       stop
    case ( 'incNS' )
       print*, ' Model ',this%modelName,' is not implemented YET'
       stop
    case ( 'visNS' )
       print*, ' Model ',this%modelName,' is not implemented YET'
       stop
    case default
       print*, ' Model ',this%modelName,' is not implemented'
       stop
    end select
  end associate

  this%imat = 0
  
end subroutine readModelData


subroutine printInitialConditions ( this , rsolfile )
  class( solver_state ), intent( in ) :: this
  character(len=*), intent( in) :: rsolfile

  if( this%type_IC == 1) then
     if( this%modelName == 'scalar' .or. this%modelName == '2eqs') then
        write(*,*) ' # IC taken from the exact solution'
     else
        write(*,*) ' # IC taken from BC'
     endif

  elseif( this%type_IC == 0) then
     write(*,*) ' # IC taken from file',rsolfile(1:len_trim(rsolfile))

  elseif( this%type_IC == -1) then
     write(*,*) ' # IC taken from file "resultsx" (piecewise constant)'

  elseif( this%type_IC == -2 .or. this%type_IC == -3) then
     write(*,*) ' # IC taken from file "dgm.sol" (DG for visualization)'

  elseif( this%type_IC > 0) then
     write(*,*) ' # IC given explicitly in subroutine SetOneElementIC(0) '

  else
     print*,'UNKNOWN type of IC'
     stop
  endif

end subroutine printInitialConditions


! subroutine initSpaceDiscMethod( this )
!   class( solver_state ), intent( inout ) :: this
!!   character(len=*), intent( in ) :: disc_space
!!   real, intent( in ) :: sigma
!!   integer, intent( in ) :: deg
!
!   allocate( Space_t :: this%space )
!
!   if(this%space%disc_space == 'NIPG') then
!      this%m_IPG = 1
!   elseif(this%space%disc_space == 'IIPG') then
!     this%m_IPG = 0
!   elseif(this%space%disc_space == 'SIPG') then
!     this%m_IPG = -1
!   elseif(this%space%disc_space == 'FEM') then
!     this%space%m_IPG = 5
!     print*,'FEM NOT TESTED !!!'
!   else
!     print*,'Bad space discretization, only SIPG, NIPG, IIPG are implemented'
!     stop
!   endif
!
!   if(this%deg < 0) then
!     print*,'Default approximation degree has to be nonnegative ',this%deg
!     stop
!   elseif( this%deg > MaxDegreeImplemented ) then
!       write(*,*) 'The default this%deg=',this%deg,&
!            ' is higher than MaxDegreeImplemented =', MaxDegreeImplemented
!       stop
!
!   else
!     write(*,'(a4,a24, a4,a7,es12.4, a4,i1)') &
!          '  # ','Space discretization by ', &
!          this%disc_space, ', c_W =',this%sigma, ", P_",this%deg
!   endif
!
! end subroutine initSpaceDiscMethod

subroutine allocTimeDiscMethod( this, disc_time, tdeg , FinTime)
  class( solver_state ), intent( inout ) :: this
  character(len=20), intent(in) :: disc_time
  integer, intent(in) :: tdeg
  real, intent(in) :: FinTime

  !print*,'subroutine allocTimeDiscMethod( '
  
  ! this%cn = .false.  ! OLD techniques
  ! this%tdg = .false. ! OLD techniques
  ! this%stdgm  = .false.

  select case(disc_time)

  case('RK')
     !     if(this%Tdeg /= 1 ) then
     !        print*,' Runge-Kutta method only with Tdeg  = 1 '
     !        stop
     !     endif
     !     this%time_method = 'E'    ! explicit
     !     if(this%Tdeg <= 0) then
     !        print*,'Degree of time degree of RK has to be >0 ',this%Tdeg
     !        stop
     !     endif
     !     write(*,'(a45,i1)') '# Time discretization by RK method order = ',this%Tdeg

     stop 'RK Not implemented in ADGo'

  case('SEMI')

     allocate(TimeBDF_t :: this%time)
     this%time%time_method = 'S'    ! semi-implicit
     !this%time%stdgm  = .false.
     this%time%disc_time = disc_time
     this%time%deg = tdeg
     this%time%FinTime = FinTime
     write(*,'(a53,i1)') &
          '# Semi-implicit time discretization by BDF, order = ', this%time%deg

  case('BDF')

     allocate(TimeBDF_t :: this%time)
     this%time%time_method = 'I'    ! implicit
     !this%time%stdgm  = .false.
     this%time%disc_time = disc_time
     this%time%deg = tdeg
     this%time%FinTime = FinTime

     if(this%time%deg <= 0) then
        print*,'Degree of time degree of BDF has to be >0 ',this%time%deg
        stop
     endif
     write(*,'(a55,i1)') &
          ' # Fully implicit time discretization by BDF, order = ',this%time%deg

  case('STDG')
     allocate( TimeTDG_t :: this%time )
     this%time%time_method = 'I'
     !this%time%stdgm  = .true.
     this%time%disc_time = disc_time
     this%time%deg = tdeg
     this%time%FinTime = FinTime

     if(this%time%deg < 0) then
        print*,'Degree of time degree of ST DG  has to be >=0 ',this%time%deg
        stop
     endif

     write(*,'(a57,i1)') &
          ' # Fully implicit time discretization by ST DG, order = ',this%time%deg

  case('CN') ! Crank - Nicolson
     !     this%cn = .true.
     !     print*,'NOT TESTED'
     stop 'CN Not implemented in ADGo'

  case('TDG') ! extrapolated STDGM
     !     this%tdg = .true.
     !     print*,'NOT TESTED'
     stop 'TDG Not implemented in ADGo'

  case default
     print*,'Unknown time integration method: ',this%time%disc_time,' ??'
     stop
  end select

end subroutine allocTimeDiscMethod

subroutine printSpaceErrorEstims( this )
  class( solver_state ), intent( in ) :: this

  if(this%space%estim_space == 'inter') then
     write(*,'(a30,a20, a14, 2es9.2, a4, f4.1)') '  # Type of error estimation: ',  &
          ' interpolation eror', &
       ', tolerances:', this%space%adapt%tol_max, this%space%adapt%tol_min,', L^q=', this%space%adapt%Lq
  else
     write(*,'(a30,a6, a14, 3es9.2)') '  # Type of error estimation: ', this%space%estim_space, &
          ', tolerances:', this%space%adapt%tol_max, this%space%adapt%tol_min! , this%space%adapt%Lq
  endif


end subroutine printSpaceErrorEstims

!> initialization of the error estims should be called in init_Problem
subroutine initErrorEstim( this )
  class( solver_state ), intent( inout ) :: this

  write(debug,*)  'initErrorEstims + initSpaceEstims are now done -> should be joined'

  ! Space errors
  allocate ( this%err(1:max_errS) )
  this%err(:) = 0.
  this%err(interLq) = 1E+30
  this%err(interL8) = 1E+30

  ! Space-time errors
  allocate(this%errSTnorm(1:max_errSTnorm), this%errSTloc(1:max_errSTnorm))

  this%errSTnorm(:) = 0.
  this%errSTloc(:)  = 0.

  if (this%time%disc_time == 'STDG') then
     allocate( this%errSnorm( 1:2,1:3) )     !L2/H1 norms in 3 time nodes
     this%errSnorm(:,:) = 0.
  endif



end subroutine initErrorEstim


subroutine initRTN( this )
  class( solver_state ), intent( inout ) :: this
  integer :: i

  write(debug,*) 'initRTN called '


  ! for local RTN functions
  allocate( this%loc_RTN(0: MaxDegreeImplemented ) )
  this%loc_RTN(:)%defined = .false.  ! will be defined during the computation


  ! standard RTN functions using reference element
  allocate( this%RTN(0: MaxRTNImplemented) )
  do i = 0, MaxRTNImplemented
     this%RTN(i)%deg = i
     this%RTN(i)%Qnum = this%space%Qdeg(i, 1) !!! OLD max(1, 2*i)
     this%RTN(i)%Gnum = i+1
     call Init_RTN_FE(this%RTN(i), &
          this%space%V_rule(this%RTN(i)%Qnum), this%space%G_rule(this%RTN(i)%Gnum) )

     allocate(this%RTN(i)%Vnode(1:maxVrule) )
     this%RTN(i)%Vnode(1:maxVrule)%def = .false.
     this%RTN(i)%Vnode(1:maxVrule)%defder = .false.
  enddo

  !print*, 'initRTN done!'


end subroutine initRTN


subroutine readNonlinearSolData ( this, non_solver, non_alg_stop, tol, max_iter, min_update )
  class( solver_state ), intent( inout ) :: this
  character(len = 20), intent( inout ) :: non_solver
  character(len = 20), intent( in ) :: non_alg_stop
  real, intent( in ) :: tol
  integer, intent( in ) :: max_iter
  integer, intent( in ) :: min_update

  if ( non_solver == 'newton') non_solver = 'Newton'

  select case( non_solver )
  case ( 'Newton' )
     allocate ( Newton_t::this%nlSolver )
     call this%nlSolver%init( non_solver, non_alg_stop, tol, max_iter, min_update )

  ! Newton with the Anderson acceleration,
  case ( 'NewtonA' )
     allocate ( Newton_t::this%nlSolver )
     call this%nlSolver%init( non_solver, non_alg_stop, tol, max_iter, min_update )

  ! Newton with the pseudo-time stepping
  case ( 'NewtonP' )
     allocate ( Newton_t::this%nlSolver )
     call this%nlSolver%init( non_solver, non_alg_stop, tol, max_iter, min_update )

  case default
     stop ' Unknown nonlinear solver! Only "Newton" or "NewtonA" nolinear solvers are implemented.'
  end select

end subroutine readNonlinearSolData

subroutine readLinearSolData( this, name, tol, mg_type )
  class( solver_state ), intent ( inout ) :: this
  character(len=20), intent( in ) :: name
  character(len=10), intent( in ) :: mg_type
  real, intent( in ) :: tol

  if ( mg_type == 'none' ) then
     ! temporarily MGsolver
     this%MGsolver = .false.
     allocate ( Gmres_t :: this%linSolver )

  else
     this%MGsolver = .true.
     !this%lin_solver = mg_type -> in linSolver%name
     allocate( MultiGrid_t :: this%linSolver )

     stop 'readLinearSolData - MG not implemented yet.'

  endif

  call this%linSolver%init( name, tol , mg_type)
  ! multigrid methods set %lin_solver from lin



end subroutine readLinearSolData

subroutine setStateParameters( this )
  class( solver_state ), intent( inout ) :: this
  integer :: i

  this%init_only_matrix_blocks = .false.
  this%local_problem = .false.
  this%RHS_presented = .true.
  this%homogenDirichlet = .false.
  this%num_flux = .true.
  this%print = .false.

  ! initialization of the time step when Stopping Criteria based on AEE are satisfied
  this%time%iter_SC = -1

  !print*, 'adaptType === ',this%space%adapt%adapt_type

  !if (this%space%estim_space == "DWR" .and. &
  !    (this%space%adapt%adapt_type == "Ahp" .or. &
  !     this%space%adapt%adapt_type == "HGhp") ) then
  !     ! FR_dofP - there is nothing like HGhp - what to do?

  this%p_mod_max = 1
  !!!this%p_mod_max = 2

  !else
  !    this%p_mod_max = 1
  !end if


  if ( this%time_dependent ) then
      this%q_mod_max = 1
      !!!this%q_mod_max = 2
  else
      this%q_mod_max = 0
  end if

  call this%setP_mod( 0 )
  call this%setQ_mod( 0 )

  write(*,*) '#  State p_mod_max (based on hp/h): ', this%p_mod_max, ' q_mod_max: ', this%q_mod_max

end subroutine setStateParameters


!> initialization of the problem - OOP version,
!> should be used one time at the beginning of computation
!>
!> creation of volume and edge quadrature rules
!> evaluation of reference test functions at quadrature rules
subroutine init_Problem( this )
  class( solver_state ) :: this

  this%wActual = .true.
  !this%wActual = .false.
  print*, "State wActual (w / wActual is used in ComputeTerms) =", this%wActual

  call this%setStateParameters()

  !init space part of the structure - quadratures etc.
  call this%space%init()

  call this%time%init()

  call this%initErrorEstim()

  call this%initRTN()

  this%print = .false.

  allocate( this%nlSolver%b(1:1), this%nlSolver%b1(1:1), &
       this%nlSolver%x(1:1), this%nlSolver%rr(1:1) )

  this%time%iter_SC = -1

  write(debug, *) 'should be AMA history in init_Problem ? '
  !AMA
  !history of AMA adaptations
  !    if(state%space%adapt%adapt_type == 'ANI' .or. state%space%adapt%adapt_type == 'Ahp'.or. state%space%adapt%adapt_type == 'Ihp') &
  !         allocate ( state%space%adapt%AMAhistory(0:max(1, state%space%adapt%max_adapt_level), 1:5) )

  this%num_call_F = 0
  this%num_call_C = 0

  ! auxiliary array for verifying the code
  allocate( this%cons(0:20, 1:2) )

  ! file with the convergence of the Newton-like method
  open( 48, file = 'criter', action="write", status="replace" )
  close( 48 )

  ! file with the convergence of the characteristic parameters of the porous media flow
  if(this%modelName == 'porous' ) then
     open( 48, file = 'PM_flow', action="write", status="replace" )
     close( 48 )
     open( 48, file = 'pm_data_IN', action="write", status="replace" )
     close( 48 )
  endif
  
  !print*,'# Problem initialized'

end subroutine init_Problem

   function getBigNSize( this ) result (nsize)
    class( solver_state ) :: this
    integer :: nsize, i, j

    if (.not. allocated( this%bigNSize ) ) &
      stop "Problem in getBigNsize - not allocated!"

    nsize = this%bigNSize( this%getP_mod(), this%getQ_mod() )

   end function getBigNSize


   subroutine setP_mod( this, i )
      class( solver_state ) :: this
      integer, intent(in) :: i

      if (i < 0 .or. i > this%p_mod_max) then
         stop "Value for p_mod out of limits!"
      else if (i == this%p_mod) then
         ! do nothing
      else
         this%p_mod = i

         !! set all flags which indicate calculations which have to be done, when p_mod is changed
         ! compute precond is needed again
         this%linSolver%precond_update = .true.
      end if

   end subroutine setP_mod

   subroutine setQ_mod( this, i )
      class( solver_state ) :: this
      integer, intent(in) :: i

      if (i < 0 .or. i > this%q_mod_max) then
         stop "Value for q_mod out of limits!"
      else if (i == this%q_mod) then
         ! do nothing
      else
         this%q_mod = i

         !! set all flags which indicate calculations which have to be done, when q_mod is changed
         ! compute precond is needed again
         this%linSolver%precond_update = .true.
      end if

   end subroutine setQ_mod

   function getP_mod( this ) result (p_mod)
      class( solver_state ) :: this
      integer :: p_mod

      p_mod = this%p_mod
   end function getP_mod

   function getQ_mod( this ) result (q_mod)
      class( solver_state ) :: this
      integer :: q_mod

      q_mod = this%q_mod
   end function getQ_mod



end module define_state


