read_my_data.f90

Go to the documentation of this file.

MODULE my_data_module
  USE def_type_mesh
  USE solve_petsc
  TYPE my_data
     LOGICAL                                 :: iformatted
     CHARACTER(LEN=200)                      :: directory, file_name
     LOGICAL                                 :: if_read_partition
     CHARACTER(LEN=3)                        :: type_pb
     INTEGER, DIMENSION(2)                   :: ndim
     LOGICAL                                 :: select_mode
     INTEGER                                 :: m_max
     INTEGER, DIMENSION(:), POINTER          :: list_mode_lect
     !===Data time iterations========================================================
     INTEGER                                 :: nb_iteration
     REAL(KIND=8)                            :: dt
     !===Data for LES================================================================
     LOGICAL                                 :: LES
     REAL(KIND=8)                            :: LES_coeff1, LES_coeff2, LES_coeff3, LES_coeff4
     LOGICAL                                 :: if_LES_in_momentum
     !===Data for precession=========================================================
     REAL(KIND=8)                            :: taux_precession, angle_s_pi
     LOGICAL                                 :: precession
     !===Data for NS penalization====================================================
     LOGICAL                                 :: if_ns_penalty
     LOGICAL                                 :: if_impose_vel_in_solids
     LOGICAL                                 :: if_compute_momentum_pseudo_force
     REAL(KIND=8)                            :: div_stab_in_ns
     !===Data for Navier-Stokes======================================================
     LOGICAL                                 :: if_navier_stokes_with_u
     LOGICAL                                 :: if_tensor_sym
     LOGICAL                                 :: if_moment_bdf2
     LOGICAL                                 :: irestart_u
     INTEGER                                 :: nb_dom_ns
     INTEGER, DIMENSION(:), POINTER          :: list_dom_ns
     REAL(KIND=8)                            :: Re
     TYPE(solver_param)                      :: my_par_vv, my_par_pp, my_par_mass
     INTEGER                                 :: pp_nb_dirichlet_sides
     INTEGER, DIMENSION(:), POINTER          :: pp_list_dirichlet_sides
     INTEGER, DIMENSION(3)                   :: vv_nb_dirichlet_sides
     TYPE(dyn_int_line), DIMENSION(3)        :: vv_list_dirichlet_sides
     INTEGER                                 :: vv_nb_dirichlet
     INTEGER                                 :: vv_nb_dirichlet_normal_velocity
     INTEGER, DIMENSION(:), POINTER          :: vv_list_dirichlet_normal_velocity_sides
     !===Data for Maxwell============================================================
     LOGICAL                                 :: if_maxwell_with_H
     LOGICAL                                 :: irestart_h
     INTEGER                                 :: nb_dom_H, nb_dom_phi, nb_inter, nb_inter_mu
     INTEGER                                 :: type_fe_H, type_fe_phi, nb_dirichlet_sides_H
     INTEGER, DIMENSION(:), POINTER          :: list_dom_H, list_dom_phi, list_inter_H_phi
     INTEGER, DIMENSION(:), POINTER          :: list_inter_mu, list_dirichlet_sides_H
     REAL(KIND=8)                            :: Rem, mu_phi
     LOGICAL                                 :: analytical_permeability
     LOGICAL                                 :: if_use_fem_integration_for_mu_bar
     LOGICAL                                 :: if_permeability_variable_in_theta
     REAL(KIND=8), DIMENSION(:), POINTER     :: mu_H, sigma
     REAL(KIND=8), DIMENSION(3)              :: stab
     TYPE(solver_param)                      :: my_par_H_p_phi
     INTEGER                                 :: phi_nb_dirichlet_sides
     INTEGER, DIMENSION(:), POINTER          :: phi_list_dirichlet_sides 
     LOGICAL                                 :: if_quasi_static_approx
     !===Data for temperature========================================================
     LOGICAL                                 :: if_temperature
     LOGICAL                                 :: irestart_T
     REAL(KIND=8)                            :: gravity_coefficient
     REAL(KIND=8), DIMENSION(:), POINTER     :: temperature_diffusivity
     REAL(KIND=8), DIMENSION(:), POINTER     :: vol_heat_capacity     
     TYPE(solver_param)                      :: my_par_temperature
     INTEGER                                 :: temperature_nb_dirichlet_sides
     INTEGER, DIMENSION(:), POINTER          :: temperature_list_dirichlet_sides
     INTEGER                                 :: nb_dom_temp 
     INTEGER, DIMENSION(:), POINTER          :: list_dom_temp
     INTEGER                                 :: nb_inter_v_T
     INTEGER, DIMENSION(:), POINTER          :: list_inter_v_T
     !===Data for level set==========================================================
     LOGICAL                                 :: if_level_set
     TYPE(solver_param)                      :: my_par_level_set
     INTEGER                                 :: level_set_nb_dirichlet_sides
     INTEGER, DIMENSION(:), POINTER          :: level_set_list_dirichlet_sides
     REAL(KIND=8)                            :: level_set_cmax
     REAL(KIND=8)                            :: level_set_comp_factor
     CHARACTER(LEN=3)                        :: level_set_reconstruction_type
     REAL(KIND=8)                            :: level_set_tanh_coeff_reconstruction
     INTEGER                                 :: nb_fluid
     REAL(KIND=8), DIMENSION(:), POINTER     :: density_fluid
     REAL(KIND=8), DIMENSION(:), POINTER     :: dyna_visc_fluid
     LOGICAL                                 :: variation_sigma_fluid
     REAL(KIND=8), DIMENSION(:), POINTER     :: sigma_fluid
     LOGICAL                                 :: if_surface_tension
     REAL(KIND=8), DIMENSION(:), POINTER     :: coeff_surface
     LOGICAL                                 :: if_mass_correction
     LOGICAL                                 :: if_kill_overshoot
     REAL(KIND=8)                            :: multiplier_for_h_min_distance
     !===Computed data for level set
     REAL(KIND=8)                            :: h_min_distance
     LOGICAL                                 :: if_level_set_P2
     LOGICAL                                 :: if_compression_mthd_JLG
     !===Data for periodicity========================================================
     TYPE(periodic_data)                     :: my_periodic
     !===Data for convergence tests==================================================
     LOGICAL                                 :: test_de_convergence
     CHARACTER(len=200)                      :: data_directory_debug
     INTEGER                                 :: numero_du_test_debug
     REAL(KIND=8), DIMENSION(4)              :: norm_ref
     !===Data for Arpack=============================================================
     LOGICAL                                 :: if_arpack
     CHARACTER(len=2)                        :: arpack_which
     INTEGER                                 :: nb_eigenvalues, arpack_iter_max
     REAL(KIND=8)                            :: arpack_tolerance
     LOGICAL                                 :: if_arpack_vtu_2d
     !===Data for stress bc==========================================================
     REAL(KIND=8)                            :: stab_bdy_ns
     !===Data for postprocessing=====================================================
     INTEGER                                 :: number_of_planes_in_real_space
     LOGICAL                                 :: check_numerical_stability
     LOGICAL                                 :: is_mesh_symmetric
     LOGICAL                                 :: if_zero_out_modes
     INTEGER                                 :: nb_select_mode_ns, nb_select_mode_mxw
     INTEGER, DIMENSION(:), POINTER          :: list_select_mode_ns, list_select_mode_mxw
     INTEGER                                 :: freq_restart, freq_en, freq_plot
     LOGICAL                                 :: verbose_timing, verbose_divergence, verbose_CFL
     LOGICAL                                 :: if_just_processing
     CONTAINS
     procedure, PUBLIC                       :: init
  END TYPE my_data
CONTAINS
  SUBROUTINE init(a)
    class(my_data), INTENT(INOUT) :: a
    !===Logicals
    a%iformatted=.false.
    a%if_read_partition=.false.
    a%select_mode=.false.
    a%LES=.false.
    a%if_LES_in_momentum=.false.
    a%precession=.false.
    a%if_ns_penalty=.false.
    a%if_impose_vel_in_solids=.false.
    a%if_compute_momentum_pseudo_force=.false.
    a%if_navier_stokes_with_u=.false.
    a%if_tensor_sym=.false.
    a%if_moment_bdf2=.false.
    a%irestart_u=.false.
    a%if_maxwell_with_H=.false.
    a%irestart_h=.false.
    a%irestart_T=.false.
    a%analytical_permeability=.false.
    a%if_use_fem_integration_for_mu_bar=.false.
    a%if_permeability_variable_in_theta=.false.
    a%if_quasi_static_approx=.false.
    a%if_temperature=.false.
    a%if_level_set=.false.
    a%variation_sigma_fluid=.false.
    a%if_surface_tension=.false.
    a%if_mass_correction=.false.
    a%if_kill_overshoot=.false.
    a%if_level_set_P2=.false.
    a%if_compression_mthd_JLG=.false.
    !===Done in sfemansinitialize. Do not touch. a%test_de_convergence
    a%if_arpack=.false.
    a%if_arpack_vtu_2d=.false.
    a%check_numerical_stability=.false.
    a%is_mesh_symmetric=.false.
    a%if_zero_out_modes=.false.
    a%verbose_timing =.false.
    a%verbose_divergence=.false. 
    a%verbose_CFL=.false.
    a%if_just_processing=.false.
    a%my_par_vv%verbose=.false.
    a%my_par_pp%verbose=.false.
    a%my_par_mass%verbose=.false.
    a%my_par_H_p_phi%verbose=.false.
    a%my_par_temperature%verbose=.false.
    a%my_par_level_set%verbose=.false.
    !===Reals
    a%LES_coeff1=0.d0
    a%LES_coeff2=0.d0 
    a%LES_coeff3=0.d0 
    a%LES_coeff4=0.d0    
    a%taux_precession=0.d0
    a%angle_s_pi=0.d0
    a%div_stab_in_ns=0.d0
    a%stab=0.d0
    a%gravity_coefficient=0.d0 
    a%level_set_cmax=0.d0
    a%level_set_comp_factor=0.d0
    a%level_set_tanh_coeff_reconstruction=0.d0
    a%norm_ref=0.d0
    a%arpack_tolerance=0.d0
    a%stab_bdy_ns=0.d0
    a%h_min_distance=0.d0
    a%multiplier_for_h_min_distance=0.d0
    !===Integers
    a%vv_nb_dirichlet=0
    a%vv_nb_dirichlet_normal_velocity=0
    a%freq_restart = 10000000
    a%freq_en =  10000000
    a%freq_plot = 10000000
  END SUBROUTINE init
END MODULE my_data_module

MODULE input_data
  USE my_data_module
  IMPLICIT NONE
  PUBLIC :: read_my_data
  TYPE(my_data), PUBLIC  :: inputs
  PRIVATE 
  
CONTAINS

  SUBROUTINE read_my_data(data_fichier)
    USE chaine_caractere
    USE my_util
    IMPLICIT NONE
    LOGICAL                        :: test
    CHARACTER(len=200), INTENT(IN) :: data_fichier
    INTEGER                        :: k
    CHARACTER(len=8), DIMENSION(3) :: vel_component=(/'uradial','utheta ','uzaxis '/)

    !===Initialize data to zero and false by default===============================
    CALL inputs%init

    !===Open data file=============================================================
    OPEN(unit = 21, file = data_fichier, form = 'formatted', status = 'unknown')

    !===Location of mesh============================================================
    CALL read_until(21, '===Is mesh file formatted (true/false)?')
    READ(21,*) inputs%iformatted
    CALL read_until(21, '===Directory and name of mesh file') 
    READ(21,*) inputs%directory, inputs%file_name

    !===Processor distribution======================================================
    CALL read_until(21, '===Number of processors in meridian section')
    READ (21,*) inputs%ndim(1)
    CALL read_until(21, '===Number of processors in Fourier space')
    READ (21,*) inputs%ndim(2)

    !===Fourier modes===============================================================
    CALL read_until(21, '===Number of Fourier modes')
    READ(21,*) inputs%m_max
    IF (mod(inputs%m_max,inputs%ndim(2))/= 0) THEN
       CALL error_petsc('BUG in read_my_data, MOD(nb_select_mode,nb_procs_F)/= 0')
    END IF
    CALL find_string(21,'===Select Fourier modes? (true/false)', test)
    IF (test) THEN
       READ (21,*) inputs%select_mode
    ELSE
       inputs%select_mode = .false.
    END IF
    IF (inputs%select_mode) THEN
       ALLOCATE(inputs%list_mode_lect(inputs%m_max))
       CALL read_until(21, '===List of Fourier modes (if select_mode=.TRUE.)')
       READ(21,*) inputs%list_mode_lect
    END IF

    !===Type of problem to be solved================================================
    CALL read_until(21, '===Problem type: (nst, mxw, mhd, fhd)')
    READ(21,*) inputs%type_pb
    IF (inputs%type_pb/='nst' .AND. inputs%type_pb/='mhd' .AND. inputs%type_pb/='mxw' .AND. inputs%type_pb/='fhd') THEN
       CALL error_petsc('BUG in read_my_data, type_pb of probleme not yet defined')
    END IF

    !===Restarts====================================================================
    CALL read_until(21, '===Restart on velocity (true/false)')
    READ (21 ,*) inputs%irestart_u
    IF (inputs%type_pb=='mhd' .OR. inputs%type_pb=='mxw' .OR. inputs%type_pb=='fhd') THEN 
       CALL read_until(21, '===Restart on magnetic field (true/false)')
       READ (21 ,*) inputs%irestart_h
    ELSE
       inputs%irestart_h = .false.
    END IF
    CALL find_string(21, '===Restart on temperature (true/false)', test)
    IF (test) THEN
       READ (21 ,*) inputs%irestart_T
    ELSE
       inputs%irestart_T = .false.
    END IF

    !===Mesh partitioning===========================================================
    CALL find_string(21, '===Do we read metis partition? (true/false)', test)
    IF (test) THEN
       READ (21,*) inputs%if_read_partition
    ELSE
       inputs%if_read_partition = .false.
    END IF
    IF (.NOT.inputs%if_read_partition .AND. (inputs%irestart_u .OR. inputs%irestart_h .OR. inputs%irestart_T)) THEN
       call error_petsc('Possibility of bug: set "===Do we read metis partition? (true/false) "' // &
            'parameter to .true. when restarting a computation since number of procs, ' // &
            'or machine, or type_pb may have changed. Make sure your mesh_part file is the correct one.')
    END IF

    !===Time integration============================================================
    CALL read_until(21, '===Time step and number of time iterations')
    READ(21,*) inputs%dt, inputs%nb_iteration

    !===Check numerical stability===================================================
    CALL find_string(21, '===Check numerical stability (true/false)', test)
    IF (test) THEN
       READ (21,*) inputs%check_numerical_stability
    ELSE
       inputs%check_numerical_stability = .false.
    END IF

    !===Periodicity=================================================================
    CALL find_string(21, '===How many pieces of periodic boundary?', test)
    IF (test) THEN
       READ (21,*) inputs%my_periodic%nb_periodic_pairs
    ELSE
       inputs%my_periodic%nb_periodic_pairs = 0
    END IF
    IF (inputs%my_periodic%nb_periodic_pairs.GE.1) THEN
       ALLOCATE(inputs%my_periodic%list_periodic(2,inputs%my_periodic%nb_periodic_pairs))
       ALLOCATE(inputs%my_periodic%vect_e(2,inputs%my_periodic%nb_periodic_pairs))
       CALL read_until(21, '===Indices of periodic boundaries and corresponding vectors')
       DO k = 1, inputs%my_periodic%nb_periodic_pairs
          READ(21,*) inputs%my_periodic%list_periodic(:,k), inputs%my_periodic%vect_e(:,k)
       END DO
    END IF

    !===Mesh symmetry===============================================================
    CALL find_string(21, '===Is the mesh symmetric (true/false)?', test)
    IF (test) THEN
       READ (21,*) inputs%is_mesh_symmetric
    ELSE
       inputs%is_mesh_symmetric = .false.
    END IF

    !===Navier Stokes data==========================================================
    IF (inputs%type_pb=='nst' .OR. inputs%type_pb=='mhd' .OR. inputs%type_pb=='fhd' .OR. inputs%irestart_u) THEN
       !==========Navier Stokes with u or m===============!
       CALL find_string(21, '===Solve Navier-Stokes with u (true) or m (false)?', test)
       IF (test) THEN
          READ(21,*) inputs%if_navier_stokes_with_u
       ELSE
          inputs%if_navier_stokes_with_u = .true.
       END IF

!!$       CALL find_string(21, '===Is tensor symmetric (true/false)?', test)
!!$       IF (test) THEN
!!$          READ(21,*) inputs%if_tensor_sym
!!$       ELSE
!!$          IF(inputs%if_navier_stokes_with_u) THEN
!!$             inputs%if_tensor_sym = .FALSE.
!!$          ELSE
!!$             inputs%if_tensor_sym = .TRUE.
!!$          END IF
!!$       END IF
       !===Tensor symmetric always used.
       inputs%if_tensor_sym = .true.

       IF (.NOT.inputs%if_navier_stokes_with_u) THEN
          CALL find_string(21, '===Do we solve momentum with bdf2 (true/false)?', test)
          IF (test) THEN
             READ(21,*) inputs%if_moment_bdf2
          ELSE
             inputs%if_moment_bdf2=.false.
          END IF
       END IF

       CALL read_until(21, '===Number of subdomains in Navier-Stokes mesh')
       READ(21,*) inputs%nb_dom_ns 
       ALLOCATE(inputs%list_dom_ns(inputs%nb_dom_ns))
       CALL read_until(21, '===List of subdomains for Navier-Stokes mesh')
       READ(21,*) inputs%list_dom_ns

       !==========Dirichlet BCs for velocity==============!
       DO k = 1, 3 
          CALL find_string(21, '===How many boundary pieces for Dirichlet BCs on '//trim(vel_component(k))//'?', test)
          IF (test) THEN
             CALL error_petsc( '===How many boundary pieces for Dirichlet BCs on '//trim(vel_component(k))//'? is disabled')
             !READ(21,*) inputs%vv_nb_dirichlet_sides(k)
          ELSE
             inputs%vv_nb_dirichlet_sides(k) = 0
          END IF
          IF (inputs%vv_nb_dirichlet_sides(k)>0) THEN
             ALLOCATE(inputs%vv_list_dirichlet_sides(k)%DIL(inputs%vv_nb_dirichlet_sides(k)))
             CALL read_until(21, '===List of boundary pieces for Dirichlet BCs on '//trim(vel_component(k)))
             READ(21,*) inputs%vv_list_dirichlet_sides(k)%DIL
          ELSE
             ALLOCATE(inputs%vv_list_dirichlet_sides(k)%DIL(0))
          END IF
       END DO
       CALL find_string(21, '===How many boundary pieces for full Dirichlet BCs on velocity?', test)
       IF (test) THEN
          READ(21,*) inputs%vv_nb_dirichlet
       ELSE
          inputs%vv_nb_dirichlet=0
       END IF
       IF (inputs%vv_nb_dirichlet>0) THEN
          DO k = 1, 3
             ALLOCATE(inputs%vv_list_dirichlet_sides(k)%DIL(inputs%vv_nb_dirichlet))
             CALL read_until(21, '===List of boundary pieces for full Dirichlet BCs on velocity')
             READ(21,*) inputs%vv_list_dirichlet_sides(k)%DIL
          END DO
       ELSE
          DO k = 1, 3
             ALLOCATE(inputs%vv_list_dirichlet_sides(k)%DIL(0))
          END DO
       END IF

       !===Homogeneous normal velocity====================!
       CALL find_string(21, '===How many boundary pieces for homogeneous normal velocity?', test)
       IF (test) THEN
          READ(21,*) inputs%vv_nb_dirichlet_normal_velocity
       ELSE
          inputs%vv_nb_dirichlet_normal_velocity=0
       END IF
       IF (inputs%vv_nb_dirichlet_normal_velocity>0) THEN
          ALLOCATE(inputs%vv_list_dirichlet_normal_velocity_sides(inputs%vv_nb_dirichlet_normal_velocity))
          CALL read_until(21, '===List of boundary pieces for homogeneous normal velocity')
          READ(21,*) inputs%vv_list_dirichlet_normal_velocity_sides
          CALL find_string(21, '===stab_bdy_ns', test)
          IF (test) THEN
             READ (21,*) inputs%stab_bdy_ns
          ELSE
             inputs%stab_bdy_ns=1.d0
          END IF
          DO k = 1, inputs%vv_nb_dirichlet_normal_velocity
             IF (minval(abs(inputs%vv_list_dirichlet_normal_velocity_sides(k) &
                  - inputs%vv_list_dirichlet_sides(1)%DIL))==0) THEN
                CALL error_petsc('Boundary conditions are mixed up')
             END IF
          END DO
       END IF

       !===Stress boundary conditions=====================! !Disabled
       CALL find_string(21, '===Stress boundary conditions? (true/false)', test) 
       IF (test) THEN
          CALL error_petsc('===Stress boundary conditions Does not work any more')
       END IF

       !==========Dirichlet BCs for pressure==============! !Disabled
       !CALL find_string(21, '===How many boundary pieces for Dirichlet BCs on pressure?', test)
       test=.false.
       IF (test) THEN
          READ(21,*) inputs%pp_nb_dirichlet_sides
       ELSE
          inputs%pp_nb_dirichlet_sides = 0
       END IF
       IF (inputs%pp_nb_dirichlet_sides>0) THEN
          ALLOCATE(inputs%pp_list_dirichlet_sides(inputs%pp_nb_dirichlet_sides))
          CALL read_until(21, '===List of boundary pieces for Dirichlet BCs on pressure')
          READ(21,*) inputs%pp_list_dirichlet_sides
       ELSE
          ALLOCATE(inputs%pp_list_dirichlet_sides(0))
       END IF

       !==========Reynolds number=========================!
       CALL read_until(21, '===Reynolds number')
       READ(21,*) inputs%Re

       !==========DIV penalty=============================!
       IF (inputs%if_navier_stokes_with_u) THEN
          CALL find_string(21, '===Coefficient for penalty of divergence in NS?', test)
          IF (test) THEN
             READ(21,*) inputs%div_stab_in_ns
          ELSE
             inputs%div_stab_in_ns = 0.d0
          END IF
       END IF

       !==========Precession==============================!
       CALL find_string(21, '===Is there a precession term (true/false)?', test)
       IF (test) THEN
          READ(21,*) inputs%precession
       ELSE
          inputs%precession = .false.
       END IF
       IF (inputs%precession) THEN
          CALL read_until(21, '===Precession rate')
          READ(21,*) inputs%taux_precession
          CALL read_until(21, '===Precession angle over pi')
          READ(21,*) inputs%angle_s_pi
       ELSE
          inputs%taux_precession = 0.d0
          inputs%angle_s_pi = 0.d0
       END IF
   
       !==========NS penalty==============================!
       CALL find_string(21, '===Use penalty in NS domain (true/false)?', test)
       IF (test) THEN
          READ(21,*) inputs%if_ns_penalty
       ELSE
          inputs%if_ns_penalty = .false.
       END IF
       IF (inputs%if_ns_penalty) THEN
          CALL find_string(21, '===Use nonzero velocity in solids (true/false)?', test)
          IF (test) THEN
             READ(21,*) inputs%if_impose_vel_in_solids
          ELSE
             inputs%if_impose_vel_in_solids = .false.
          END IF
          CALL find_string(21, '===Compute z momentum (true/false)?', test)
          IF (test) THEN
             READ(21,*) inputs%if_compute_momentum_pseudo_force
          ELSE
             inputs%if_compute_momentum_pseudo_force = .false.
          END IF
       ELSE
          inputs%if_impose_vel_in_solids = .false.
          inputs%if_compute_momentum_pseudo_force = .false.
       END IF

       !==========Solver parameters for velocity==========!
       CALL read_until(21, '===Maximum number of iterations for velocity solver')
       READ(21,*) inputs%my_par_vv%it_max
       CALL read_until(21, '===Relative tolerance for velocity solver')
       READ(21,*) inputs%my_par_vv%rel_tol
       CALL read_until(21, '===Absolute tolerance for velocity solver')
       READ(21,*) inputs%my_par_vv%abs_tol
       CALL find_string(21, '===Velocity solver verbose? (true/false)', test)
       IF (test) THEN
          READ(21,*) inputs%my_par_vv%verbose
       END IF
       CALL read_until(21, '===Solver type for velocity (FGMRES, CG, ...)')
       READ(21,*) inputs%my_par_vv%solver
       CALL read_until(21, '===Preconditionner type for velocity solver (HYPRE, JACOBI, MUMPS...)')
       READ(21,*) inputs%my_par_vv%precond

       !==========Solver parameters for pressure==========!
       CALL read_until(21, '===Maximum number of iterations for pressure solver')
       READ(21,*) inputs%my_par_pp%it_max
       CALL read_until(21, '===Relative tolerance for pressure solver')
       READ(21,*) inputs%my_par_pp%rel_tol
       CALL read_until(21, '===Absolute tolerance for pressure solver')
       READ(21,*) inputs%my_par_pp%abs_tol
       CALL find_string(21, '===Pressure solver verbose? (true/false)',test)
       IF (test) THEN
          READ(21,*) inputs%my_par_pp%verbose
       END IF
       CALL read_until(21, '===Solver type for pressure (FGMRES, CG, ...)')
       READ(21,*) inputs%my_par_pp%solver
       CALL read_until(21, '===Preconditionner type for pressure solver (HYPRE, JACOBI, MUMPS...)')
       READ(21,*) inputs%my_par_pp%precond

       !==========Solver parameters for mass matrix=======!
       CALL read_until(21, '===Maximum number of iterations for mass matrix solver')
       READ(21,*) inputs%my_par_mass%it_max
       CALL read_until(21, '===Relative tolerance for mass matrix solver')
       READ(21,*) inputs%my_par_mass%rel_tol
       CALL read_until(21, '===Absolute tolerance for mass matrix solver')
       READ(21,*) inputs%my_par_mass%abs_tol
       CALL find_string(21, '===Mass matrix solver verbose? (true/false)',test)
       IF (test) THEN
          READ(21,*) inputs%my_par_mass%verbose
       END IF
       CALL read_until(21, '===Solver type for mass matrix (FGMRES, CG, ...)')
       READ(21,*) inputs%my_par_mass%solver
       CALL read_until(21, '===Preconditionner type for mass matrix solver (HYPRE, JACOBI, MUMPS...)')
       READ(21,*) inputs%my_par_mass%precond

       !==========LES coefficients========================!
       CALL find_string(21, '===Use LES? (true/false)', test)
       IF (test) THEN
          READ(21,*) inputs%LES
       ELSE
          inputs%LES = .false.
       END IF

       IF (inputs%LES) THEN
          CALL find_string(21, '===Coefficients for LES', test)
          IF (test) THEN
             CALL error_petsc('===Coefficients for LES is disabled')
          END IF
          !===LES_coeff4 is fixed and coeff3 is disabled
          inputs%LES_coeff4 = 0.125d0  
          inputs%LES_coeff3 = 0.d0
          !===LES_coeff2 in [0.15,1.0] for viscosity entropy
          !===LES_coeff2 is equal to 1.d10 for test with first order viscosity
          CALL read_until(21, '===Coefficient multiplying residual')
          READ(21,*) inputs%LES_coeff2
          !===LES_coeff1 is equal to LES_coeff4*MAX(velocity)=0.125d0*MAX(velocity)
          CALL find_string(21, '===Coefficient for explicit LES', test)
          IF (test) THEN
             READ(21,*) inputs%LES_coeff1
          ELSE
             inputs%LES_coeff1 = 0.d0
          END IF
        ELSE
          inputs%LES_coeff1 = 0.d0
          inputs%LES_coeff2 = 0.d0
          inputs%LES_coeff3 = 0.d0
          inputs%LES_coeff4 = 0.d0
       END IF
    END IF
    IF (.NOT.inputs%if_navier_stokes_with_u) THEN
       CALL find_string(21, '===Use LES in momentum? (true/false)', test)
       IF (test) THEN
          READ(21,*) inputs%if_LES_in_momentum
       ELSE
          inputs%if_LES_in_momentum=.true.
       END IF
    END IF

    !===Maxwell data================================================================
    IF (inputs%type_pb=='mxw' .OR. inputs%type_pb=='mhd' .OR. inputs%type_pb=='fhd') THEN
       !==========Maxwell with H or B=====================!
       CALL find_string(21, '===Solve Maxwell with H (true) or B (false)?', test)
       IF (test) THEN
          READ(21,*) inputs%if_maxwell_with_H
       ELSE
          inputs%if_maxwell_with_H = .false.
       END IF

       !==========Quasi-static approximation==============!
       CALL find_string(21, '===Quasi-static approximation (true) or (false)?', test)
       IF (test) THEN
          READ(21,*) inputs%if_quasi_static_approx
       ELSE
          inputs%if_quasi_static_approx = .false.
       END IF

       !==========Subdomains for H========================!
       CALL read_until(21, '===Number of subdomains in magnetic field (H) mesh')
       READ(21,*) inputs%nb_dom_H  ! number of sub_domains for H
       ALLOCATE(inputs%list_dom_H(inputs%nb_dom_H))
       CALL read_until(21, '===List of subdomains for magnetic field (H) mesh')
       READ(21,*) inputs%list_dom_H

       !==========Interfaces H/H==========================!
       CALL read_until(21, '===Number of interfaces in H mesh')
       READ(21,*) inputs%nb_inter_mu
       ALLOCATE(inputs%list_inter_mu(inputs%nb_inter_mu))
       IF (inputs%nb_inter_mu>0) THEN
          CALL read_until(21, '===List of interfaces in H mesh')
          READ(21,*) inputs%list_inter_mu
       END IF

       !==========Dirichlet BCs for H=====================!
       CALL read_until(21, '===Number of Dirichlet sides for Hxn')
       READ(21,*) inputs%nb_dirichlet_sides_H
       ALLOCATE(inputs%list_dirichlet_sides_H(inputs%nb_dirichlet_sides_H))
       IF (inputs%nb_dirichlet_sides_H>0) THEN
          CALL read_until(21, '===List of Dirichlet sides for Hxn')
          READ(21,*) inputs%list_dirichlet_sides_H
       END IF

       !==========Permeability for H======================!
       CALL find_string(21, '===Is permeability defined analytically (true/false)?', test)
       IF (test) THEN
          READ(21,*) inputs%analytical_permeability
       ELSE
          inputs%analytical_permeability = .false.
       END IF
       IF (.NOT.inputs%analytical_permeability) THEN
          CALL read_until(21, '===Permeability in the conductive part (1:nb_dom_H)')
          ALLOCATE(inputs%mu_H(inputs%nb_dom_H))
          READ(21,*) inputs%mu_H
       END IF
       inputs%if_permeability_variable_in_theta = .false.
       IF (inputs%analytical_permeability) THEN
          CALL find_string(21, '===Is permeability variable in theta (true/false)?', test)
          IF (test) THEN
             READ(21,*) inputs%if_permeability_variable_in_theta
          END IF
       END IF

       !==========FEM or Gaussian integration for mu_bar==!
       inputs%if_use_fem_integration_for_mu_bar = .true.
       IF (inputs%analytical_permeability) THEN
          CALL find_string(21, '===Use FEM Interpolation for magnetic permeability (true/false)?', test)
          IF (test) THEN
             READ(21,*) inputs%if_use_fem_integration_for_mu_bar
          END IF
       END IF

       !==========Conductivity for H======================!
       CALL read_until(21, '===Conductivity in the conductive part (1:nb_dom_H)')
       ALLOCATE(inputs%sigma(inputs%nb_dom_H))
       READ(21,*) inputs%sigma

       !==========Finite element type=====================!
       CALL read_until(21, '===Type of finite element for magnetic field')
       READ(21,*) inputs%type_fe_H

       !==========Magnetic Reynolds number================!
       CALL read_until(21, '===Magnetic Reynolds number')
       READ(21,*) inputs%Rem

       !==========Stabilization parameters================!
       CALL read_until(21, '===Stabilization coefficient (divergence)')
       READ(21,*) inputs%stab(1)
       IF (inputs%nb_inter_mu>0 .OR. inputs%nb_dirichlet_sides_H>0) THEN
          CALL read_until(21, '===Stabilization coefficient for Dirichlet H and/or interface H/H')
          READ(21,*) inputs%stab(3)
       ELSE
          inputs%stab(3) = 0.d0
       END IF

       !==========Subdomains for phi======================!
       CALL find_string(21, '===Number of subdomains in magnetic potential (phi) mesh', test)
       IF (test) THEN
          READ (21,*) inputs%nb_dom_phi
       ELSE
          inputs%nb_dom_phi = 0
       END IF
       ALLOCATE(inputs%list_dom_phi(inputs%nb_dom_phi))
       IF (inputs%nb_dom_phi>0) THEN
          !==========List of subdomains for phi======================!
          CALL read_until(21, '===List of subdomains for magnetic potential (phi) mesh')
          READ(21,*) inputs%list_dom_phi

          !==========Dirichlet BCs on phi====================!
          CALL read_until(21, '===How many boundary pieces for Dirichlet BCs on phi?')
          READ(21,*) inputs%phi_nb_dirichlet_sides
          IF (inputs%phi_nb_dirichlet_sides>0) THEN
             ALLOCATE(inputs%phi_list_dirichlet_sides(inputs%phi_nb_dirichlet_sides))
             CALL read_until(21, '===List of boundary pieces for Dirichlet BCs on phi')
             READ(21,*) inputs%phi_list_dirichlet_sides
          ELSE
             ALLOCATE(inputs%phi_list_dirichlet_sides(0))
          END IF

          !==========H/phi interface=========================!
          CALL read_until(21, '===Number of interfaces between H and phi')
          READ(21,*) inputs%nb_inter  ! number of interfaces between H and phi
          ALLOCATE(inputs%list_inter_H_phi(inputs%nb_inter))
          CALL read_until(21, '===List of interfaces between H and phi')
          IF (inputs%nb_inter>0) THEN
             READ(21,*) inputs%list_inter_H_phi
          END IF

          !==========Permeability in vacuum==================!
          CALL read_until(21, '===Permeability in vacuum')
          READ(21,*) inputs%mu_phi

          !==========Finite element type=====================!
          CALL read_until(21, '===Type of finite element for scalar potential')
          READ(21,*) inputs%type_fe_phi

          !==========Stabilization parameters================!
          CALL read_until(21, '===Stabilization coefficient (interface H/phi)')
          READ(21,*) inputs%stab(2)
       END IF

       !==========Solver parameters=======================!
       CALL read_until(21, '===Maximum number of iterations for Maxwell solver')
       READ(21,*) inputs%my_par_H_p_phi%it_max
       CALL read_until(21, '===Relative tolerance for Maxwell solver')
       READ(21,*) inputs%my_par_H_p_phi%rel_tol
       CALL read_until(21, '===Absolute tolerance for Maxwell solver')
       READ(21,*) inputs%my_par_H_p_phi%abs_tol
       CALL find_string(21, '===Maxwell solver verbose? (true/false)',test)
       IF (test) THEN
          READ(21,*) inputs%my_par_H_p_phi%verbose
       END IF
       CALL read_until(21, '===Solver type for Maxwell (FGMRES, CG, ...)')
       READ(21,*) inputs%my_par_H_p_phi%solver
       CALL read_until(21, '===Preconditionner type for Maxwell solver (HYPRE, JACOBI, MUMPS...)')
       READ(21,*) inputs%my_par_H_p_phi%precond
    END IF

    !===Data temperature==============================================================
    CALL find_string(21, '===Is there a temperature field?', test)
    IF (test) THEN
       READ (21,*) inputs%if_temperature
    ELSE
       inputs%if_temperature = .false.
    END IF
    IF (inputs%if_temperature) THEN
       !==========Gravity coefficient for temperature=====!
       CALL find_string(21, '===Non-dimensional gravity coefficient', test)
       IF (test) THEN
          READ (21,*) inputs%gravity_coefficient
       ELSE
          inputs%gravity_coefficient = 0.d0
       END IF
       !==========Subdomains for temp=====================!
       CALL read_until(21, '===Number of subdomains in temperature mesh')
       READ(21,*) inputs%nb_dom_temp  ! number of sub_domains for temp
       ALLOCATE(inputs%list_dom_temp(inputs%nb_dom_temp))
       CALL read_until(21, '===List of subdomains for temperature mesh')
       READ(21,*) inputs%list_dom_temp
       ALLOCATE(inputs%vol_heat_capacity(inputs%nb_dom_temp)) ! Two choices for the user, indicate 1) the heat capacity and the conductivity (necessary if non uniform heat capacity) or 2) the diffusivity only
       ALLOCATE(inputs%temperature_diffusivity(inputs%nb_dom_temp)) ! In both cases, temperature_diffusivity is used, it contains the conductivities in 1) and the diffusivities in 2)
       CALL find_string(21, '===Volumetric heat capacity (1:nb_dom_temp)', test)
       IF (test) THEN ! Case 1)
          READ(21,*) inputs%vol_heat_capacity
          CALL read_until(21, '===Thermal conductivity (1:nb_dom_temp)')
          READ(21,*) inputs%temperature_diffusivity
       ELSE ! Case 2)        
          inputs%vol_heat_capacity = 1.d0 ! Heat capacity is equalized to one so that it does not impact the calculus
          CALL read_until(21, '===Diffusivity coefficient for temperature (1:nb_dom_temp)')
          READ(21,*) inputs%temperature_diffusivity
       END IF
       !==========Dirichlet BCs on temperature============!
       CALL read_until(21, '===How many boundary pieces for Dirichlet BCs on temperature?')
       READ(21,*) inputs%temperature_nb_dirichlet_sides
       IF (inputs%temperature_nb_dirichlet_sides>0) THEN
          ALLOCATE(inputs%temperature_list_dirichlet_sides(inputs%temperature_nb_dirichlet_sides))
          CALL read_until(21, '===List of boundary pieces for Dirichlet BCs on temperature')
          READ(21,*) inputs%temperature_list_dirichlet_sides
       ELSE
          ALLOCATE(inputs%temperature_list_dirichlet_sides(0))
       END IF
       !==========Interfaces between vel and temp=========!
       CALL find_string(21, '===Number of interfaces between velocity and temperature only domains (for nst applications)', test)
       IF (test) THEN
          READ(21,*) inputs%nb_inter_v_T
          ALLOCATE(inputs%list_inter_v_T(inputs%nb_inter_v_T))
          CALL read_until(21, '===List of interfaces between velocity and temperature only domains (for nst applications)')
          READ(21,*) inputs%list_inter_v_T
       ELSE
          ALLOCATE(inputs%list_inter_v_T(0))
       END IF
       !==========Solver parameters=======================!
       CALL read_until(21, '===Maximum number of iterations for temperature solver')
       READ(21,*) inputs%my_par_temperature%it_max
       CALL read_until(21, '===Relative tolerance for temperature solver')
       READ(21,*) inputs%my_par_temperature%rel_tol
       CALL read_until(21, '===Absolute tolerance for temperature solver')
       READ(21,*) inputs%my_par_temperature%abs_tol
       CALL find_string(21, '===Temperature solver verbose? (true/false)',test)
       IF (test) THEN
          READ(21,*) inputs%my_par_temperature%verbose
       END IF
       CALL read_until(21, '===Solver type for temperature (FGMRES, CG, ...)')
       READ(21,*) inputs%my_par_temperature%solver
       CALL read_until(21, '===Preconditionner type for temperature solver (HYPRE, JACOBI, MUMPS...)')
       READ(21,*) inputs%my_par_temperature%precond
    END IF

    !===Data Level set================================================================
    CALL find_string(21, '===Is there a level set?', test)
    IF (test) THEN
       READ (21,*) inputs%if_level_set
    ELSE
       inputs%if_level_set = .false.
    END IF
    IF (inputs%if_level_set) THEN
       !==========Number of fluids========================!
       CALL read_until(21, '===How many fluids?')
       READ(21,*) inputs%nb_fluid
       !==========Level_set multiplier for h_min==========!
       CALL read_until(21, '===multiplier for h_min for level set')
       READ(21,*) inputs%multiplier_for_h_min_distance
       !==========Level_set c_max=========================!
       !===Disabled functionality
       inputs%level_set_cmax=0.d0
       !==========Level_set compression factor============!
       CALL read_until(21, '===Compression factor for level set')
       READ(21,*) inputs%level_set_comp_factor
       !==========Densities of fluids=====================!
       ALLOCATE(inputs%density_fluid( inputs%nb_fluid))
       CALL read_until(21, '===Density of fluid 0, fluid 1, ...')
       READ(21,*) inputs%density_fluid
       !==========Dynamic viscosities of fluids===========!
       ALLOCATE(inputs%dyna_visc_fluid( inputs%nb_fluid))
       CALL read_until(21, '===Dynamic viscosity of fluid 0, fluid 1, ...')
       READ(21,*) inputs%dyna_visc_fluid
       !==========Conductivities of fluids================!
       IF (inputs%type_pb=='mhd' .OR. inputs%type_pb=='fhd') THEN
          CALL find_string(21, '===Conductivity of fluid 0, fluid 1, ...', test)
          IF (test) THEN
             inputs%variation_sigma_fluid =.true.
             ALLOCATE(inputs%sigma_fluid(inputs%nb_fluid))
             READ(21,*) inputs%sigma_fluid
          ELSE
             inputs%variation_sigma_fluid =.false.
          END IF
       END IF
       !==========Surface tension=========================!
       CALL find_string(21, '===Is there a surface tension?', test)
       IF (test) THEN
          READ (21,*) inputs%if_surface_tension
       ELSE
          inputs%if_surface_tension = .false.
       END IF
       IF (inputs%if_surface_tension) THEN
          !==========Coefficient for surface tension======!
          CALL read_until(21, '===Coefficients of surface tension for level set 0, level set 1, ...')
          ALLOCATE(inputs%coeff_surface(inputs%nb_fluid-1))
          READ(21,*) inputs%coeff_surface
       END IF
       !==========Mass correction=========================!
       CALL find_string(21, '===Do we apply mass correction? (true/false)', test)
       IF (test) THEN
          READ (21,*) inputs%if_mass_correction
       ELSE
          inputs%if_mass_correction=.true.
       END IF
       !==========Dirichlet BCs on level_set==============!
       CALL read_until(21, '===How many boundary pieces for Dirichlet BCs on level set?')
       READ(21,*) inputs%level_set_nb_dirichlet_sides
       IF (inputs%level_set_nb_dirichlet_sides>0) THEN
          ALLOCATE(inputs%level_set_list_dirichlet_sides(inputs%level_set_nb_dirichlet_sides))
          CALL read_until(21, '===List of boundary pieces for Dirichlet BCs on level set')
          READ(21,*) inputs%level_set_list_dirichlet_sides
       ELSE
          ALLOCATE(inputs%level_set_list_dirichlet_sides(0))
       END IF
       !==========Solver parameters=======================!
       CALL read_until(21, '===Maximum number of iterations for level set solver')
       READ(21,*) inputs%my_par_level_set%it_max
       CALL read_until(21, '===Relative tolerance for level set solver')
       READ(21,*) inputs%my_par_level_set%rel_tol
       CALL read_until(21, '===Absolute tolerance for level set solver')
       READ(21,*) inputs%my_par_level_set%abs_tol
       CALL find_string(21, '===Level set solver verbose? (true/false)',test)
       IF (test) THEN
          READ(21,*) inputs%my_par_level_set%verbose
       END IF
       CALL read_until(21, '===Solver type for level set (FGMRES, CG, ...)')
       READ(21,*) inputs%my_par_level_set%solver
       CALL read_until(21, '===Preconditionner type for level set solver (HYPRE, JACOBI, MUMPS...)')
       READ(21,*) inputs%my_par_level_set%precond

       !==========Reconstruction parameters===============!
       CALL find_string(21, '===How are the variables reconstructed from the level set function? (lin, reg)', test)
       IF (test) THEN
          READ (21,*) inputs%level_set_reconstruction_type
       ELSE
          inputs%level_set_reconstruction_type = 'reg'
       END IF
       IF (inputs%level_set_reconstruction_type=='reg') THEN
          CALL find_string(21, '===Value of the regularization coefficient in (0,0.5]', test)
          IF (test) THEN
             READ (21,*) inputs%level_set_tanh_coeff_reconstruction
             IF (inputs%level_set_tanh_coeff_reconstruction.LE.1d-2 .OR. &
                  .5d0 .LE. inputs%level_set_tanh_coeff_reconstruction) THEN
                inputs%level_set_tanh_coeff_reconstruction = 0.45d0
             END IF
          ELSE
             inputs%level_set_tanh_coeff_reconstruction = 0.45d0
          END IF
       ELSE IF (inputs%level_set_reconstruction_type/='lin') THEN
          CALL error_petsc('BUG in read_my_data: variable reconstruction type not correct')
       END IF
       !==========Level set overshoot=====================!
       IF(inputs%level_set_reconstruction_type/='reg') THEN
          CALL find_string(21, '===Do we kill level set overshoot? (true/false)', test)
          IF (test) THEN
             READ (21,*) inputs%if_kill_overshoot
          ELSE
             inputs%if_kill_overshoot=.false.
          END IF
       ELSE
          inputs%if_kill_overshoot=.false.
       END IF

       !==========Level set type finite element===========!
       CALL find_string(21, '===Do we use P2 finite element for level_set? (true/false)', test)
       IF (test) THEN
          READ(21,*) inputs%if_level_set_P2
       ELSE
          inputs%if_level_set_P2=.false.
       END IF

       !==========Level set compression method============!
       CALL find_string(21, '===Do we use JLG compression method for level_set? (true/false)', test)
       IF (test) THEN
          READ(21,*) inputs%if_compression_mthd_JLG
       ELSE
          inputs%if_compression_mthd_JLG=.true.
       END IF

       !==========Check coherence=========================!
       IF (inputs%type_pb=='mxw') THEN
          CALL error_petsc('BUG in read_my_data: Level_set without Navier-Stokes')
       END IF
    END IF

    !===Data for arpack (eigenvalue problems)=======================================
    !==========Frequency parameters====================!
    IF (inputs%type_pb=='mxw') THEN
       CALL find_string(21, '===Do we use Arpack?', test)
       IF (test) THEN
          READ (21,*) inputs%if_arpack
       ELSE
          inputs%if_arpack = .false.
       END IF
       IF (inputs%if_arpack)THEN 
          CALL read_until(21, '===Number of eigenvalues to compute')
          READ (21,*) inputs%nb_eigenvalues
          CALL read_until(21, '===Maximum number of Arpack iterations')
          READ (21,*) inputs%arpack_iter_max
          CALL read_until(21, '===Tolerance for Arpack')
          READ (21,*) inputs%arpack_tolerance
          CALL read_until(21, &
               '===Which eigenvalues (''LM'', ''SM'', ''SR'', ''LR'' ''LI'', ''SI'')')
          READ (21,*) inputs%arpack_which
          CALL find_string(21, '===Create 2D vtu files for Arpack? (true/false)', test)
          IF (test) THEN
             READ (21,*) inputs%if_arpack_vtu_2d
          ELSE
             inputs%if_arpack_vtu_2d = .false.
          END IF
       END IF
    ELSE
       !==========Arpack currently works only for Maxwell=!
       inputs%if_arpack = .false.
    END IF

    !===Data for post processing====================================================
    CALL find_string(21, '===Number of planes in real space for Visualization', test)
    IF (test) THEN
       READ (21,*) inputs%number_of_planes_in_real_space
    ELSE
       inputs%number_of_planes_in_real_space = 10
    END IF

    !==========Frequency parameters====================!
    CALL find_string(21, '===Frequency to write restart file', test)
    IF (test) THEN
       READ (21,*) inputs%freq_restart
    ELSE
       inputs%freq_restart = 100000000
    END IF
    CALL find_string(21, '===Frequency to write energies', test)
    IF (test) THEN
       READ (21,*) inputs%freq_en
    ELSE
       inputs%freq_en = 100000000
    END IF
    CALL find_string(21, '===Frequency to create plots', test)
    IF (test) THEN
       READ (21,*) inputs%freq_plot
    ELSE
       inputs%freq_plot = 100000000
    END IF
    CALL find_string(21, '===Just postprocessing without computing? (true/false)', test)
    IF (test) THEN
       READ (21,*) inputs%if_just_processing
    ELSE
       inputs%if_just_processing = .false.
    END IF


    !==========Modes to be zeroed out==================!
    IF (inputs%type_pb=='mhd') THEN 
       CALL find_string(21, '===Should some modes be zeroed out?', test)
       IF (test) THEN
          READ (21,*) inputs%if_zero_out_modes
       ELSE
          inputs%if_zero_out_modes = .false.
       END IF
       IF (inputs%if_zero_out_modes) THEN
          CALL read_until(21, '===How many Navier-Stokes modes to zero out?')
          READ(21,*) inputs%nb_select_mode_ns
          ALLOCATE(inputs%list_select_mode_ns(inputs%nb_select_mode_ns))
          CALL read_until(21, '===List of Navier-Stokes modes to zero out?')
          READ(21,*) inputs%list_select_mode_ns
          CALL read_until(21, '===How Maxwell modes to zero out?')
          READ(21,*) inputs%nb_select_mode_mxw
          ALLOCATE(inputs%list_select_mode_mxw(inputs%nb_select_mode_mxw))
          CALL read_until(21, '===List of Maxwell modes to zero out?')
          READ(21,*) inputs%list_select_mode_mxw
       END IF
    END IF

    !==========Verbose=================================!
    CALL find_string(21, '===Verbose timing? (true/false)', test)
    IF (test) THEN
       READ (21,*) inputs%verbose_timing
    ELSE
       inputs%verbose_timing = .false.
    END IF
    CALL find_string(21, '===Verbose divergence? (true/false)', test)
    IF (test) THEN
       READ (21,*) inputs%verbose_divergence
    ELSE
       inputs%verbose_divergence = .false.
    END IF
    CALL find_string(21, '===Verbose CFL? (true/false)', test)
    IF (test) THEN
       READ (21,*) inputs%verbose_CFL
    ELSE
       inputs%verbose_CFL = .false.
    END IF

    !===Norms for reference tests===================================================
    IF (inputs%test_de_convergence) THEN
       CALL read_until(21, '===Reference results')
       DO k = 1, 4
          READ(21,*) inputs%norm_ref(k)
       END DO
    ELSE
       inputs%norm_ref = 1.d0
    END IF

    !===mxw becomes mxx if restart_u and mxw========================================
    IF (inputs%type_pb == 'mxw' .AND. inputs%irestart_u) THEN
       inputs%type_pb = 'mxx'
    END IF

    CLOSE(21)

    !===Check coherence of data=====================================================
    CALL check_coherence_of_data
    RETURN

  END SUBROUTINE read_my_data

  SUBROUTINE check_coherence_of_data
    USE my_util
    IMPLICIT NONE
    INTEGER :: k, mode_min, mode_max, delta_mode
    LOGICAL :: test

    !===Dirichlet BCs for Navier-Stokes=============================================
    IF (inputs%my_periodic%nb_periodic_pairs>0) THEN
       IF (inputs%type_pb=='nst' .OR. inputs%type_pb=='mhd' .OR. inputs%irestart_u) THEN
          !==========Velocity================================!
          DO k = 1, 3
             IF (inputs%vv_nb_dirichlet_sides(k)<1) cycle
             test = check_coherence_with_periodic_bcs(inputs%vv_list_dirichlet_sides(k)%DIL)
             IF (test) THEN
                CALL error_petsc(' BUG in read_my_data: Incompatible Dirichlet'// &
                     ' and periodic BCs on velocity')
             END IF
          END DO
          !==========Pressure================================!
          IF (inputs%pp_nb_dirichlet_sides>0) THEN
             test = check_coherence_with_periodic_bcs(inputs%pp_list_dirichlet_sides)
             IF (test) THEN
                CALL error_petsc(' BUG in read_my_data: Incompatible Dirichlet'// &
                     ' and periodic BCs on pressure')
             END IF
          END IF
       END IF
    END IF

    !===No penalty with momentum formulation========================================
    IF (.NOT.inputs%if_navier_stokes_with_u .AND. inputs%if_ns_penalty) THEN
        CALL error_petsc('BUG in read_my_data: No penalty with momentum formulation')
     END IF

    !===Dirichlet BCs for temperature for Navier-Stokes=============================
    IF (inputs%if_temperature.AND. inputs%my_periodic%nb_periodic_pairs>0) THEN
       IF (inputs%temperature_nb_dirichlet_sides>0) THEN
          test = check_coherence_with_periodic_bcs(inputs%temperature_list_dirichlet_sides)
          IF (test) THEN
             CALL error_petsc(' BUG in read_my_data: Incompatible Dirichlet'// &
                  ' and periodic BCs on temperature')
          END IF
       END IF
    END IF

    !===Dirichlet BCs for Level_set for Navier-Stokes=============================
    IF (inputs%if_level_set.AND. inputs%my_periodic%nb_periodic_pairs>0) THEN
       IF (inputs%level_set_nb_dirichlet_sides>0) THEN
          test = check_coherence_with_periodic_bcs(inputs%level_set_list_dirichlet_sides)
          IF (test) THEN
             CALL error_petsc(' BUG in read_my_data: Incompatible Dirichlet'// &
                  ' and periodic BCs on level_set')
          END IF
       END IF
    END IF

    !===Dirichlet BCs magnetic field for Maxwell====================================
    IF (inputs%my_periodic%nb_periodic_pairs>0) THEN
       IF (inputs%type_pb=='mxw' .OR. inputs%type_pb=='mhd' .OR. inputs%type_pb=='fhd') THEN
          !==========Magnetic field==========================!
          IF (inputs%nb_dirichlet_sides_H>0) THEN
             test = check_coherence_with_periodic_bcs(inputs%list_dirichlet_sides_H)
             IF (test) THEN
                CALL error_petsc(' BUG in read_my_data: Incompatible Dirichlet'// &
                     ' and periodic BCs on magnetic field')
             END IF
          END IF

          !==========Scalar potential========================!
          IF (inputs%nb_dom_phi>0 .AND. inputs%phi_nb_dirichlet_sides>0) THEN
             test = check_coherence_with_periodic_bcs(inputs%phi_list_dirichlet_sides)
             IF (test) THEN
                CALL error_petsc(' BUG in read_my_data: Incompatible Dirichlet'// &
                     ' and periodic BCs on scalar potential')
             END IF
          END IF
       END IF
    END IF

    !===Check temperature with Maxwell==============================================
    IF (inputs%my_periodic%nb_periodic_pairs>0) THEN
       IF (inputs%if_temperature .AND. inputs%type_pb=='mxw') THEN
          CALL error_petsc('Bug in read_my_data: incompatible temperature with maxwell')
       END IF
    END IF

    !===Check Arpack================================================================
    IF (inputs%if_arpack) THEN
       IF (inputs%ndim(2) /= inputs%m_max) THEN
          CALL error_petsc('Bug in read_my_data: #Fourier modes'// &
               ' not equal to #processors in Fourier direction')
       END IF
    END IF

    !===Check Fourier modes=========================================================
    IF (inputs%select_mode .AND. .NOT.inputs%if_arpack) THEN
       IF (SIZE(inputs%list_mode_lect)/=1) THEN
          mode_max = maxval(inputs%list_mode_lect)
          mode_min = minval(inputs%list_mode_lect)
          IF (mod(mode_max-mode_min,SIZE(inputs%list_mode_lect)-1)/=0) THEN
             CALL error_petsc('Bug in read_my_data: Fourier modes not equally spaced ')
          END IF
          delta_mode = (mode_max-mode_min)/(SIZE(inputs%list_mode_lect)-1)
          DO k = 0, SIZE(inputs%list_mode_lect)-1
             IF (minval(abs(inputs%list_mode_lect-(delta_mode*k+mode_min)))/=0) THEN
                CALL error_petsc('Bug in read_my_data: Fourier modes not equally spaced ')
             END IF
          END DO
          DO k = 1, SIZE(inputs%list_mode_lect)-1
             IF (inputs%list_mode_lect(k+1).LE.inputs%list_mode_lect(k)) THEN
                CALL error_petsc('Bug in read_my_data: Fourier modes not in increasing order ')
             END IF
          END DO
       END IF
    END IF

    !===Irestart for postprocessing=================================================
    IF (inputs%if_just_processing) THEN   
       inputs%irestart_u = .false.
       inputs%irestart_h = .false.
       IF (inputs%type_pb/='mxw') inputs%irestart_u = .true.
       IF (inputs%type_pb/='nst') inputs%irestart_h = .true.
    END IF

    !===Allocate dummy list_dom_* for new partitioning==============================
    IF (.NOT. ASSOCIATED(inputs%list_dom_ns))  ALLOCATE(inputs%list_dom_ns(0))
    IF (.NOT. ASSOCIATED(inputs%list_dom_H))   ALLOCATE(inputs%list_dom_H(0))
    IF (.NOT. ASSOCIATED(inputs%list_dom_phi)) ALLOCATE(inputs%list_dom_phi(0))
    IF (.NOT. ASSOCIATED(inputs%list_dom_temp))   ALLOCATE(inputs%list_dom_temp(0))

  END SUBROUTINE check_coherence_of_data

  FUNCTION check_coherence_with_periodic_bcs(list) RESULT(test)
    IMPLICIT NONE
    INTEGER, DIMENSION(:) :: list
    LOGICAL               :: test
    INTEGER               :: k, n

    test = .false.
    DO k = 1, inputs%my_periodic%nb_periodic_pairs
       DO n = 1, SIZE(list)
          IF (minval(abs(list(n)-inputs%my_periodic%list_periodic(:,k)))==0) THEN
             test = .true.
             RETURN
          END IF
       END DO
    END DO
  END FUNCTION check_coherence_with_periodic_bcs

END MODULE input_data