






















































































































 


































































































































































 

















































































































































































































































































































































































































































 






































































































 










 








































































































































































































































 










                                
                                

                                
                                

                                
                                

                                





















































































 

subroutine eigenproblem_marquina(handle,rhor,velxr,velyr,&
     velzr,epsr,rhol,velxl,velyl,velzl,epsl,gxx,gxy,gxz,&
     gyy,gyz,gzz,u,det,alp,beta,densl,sxl,syl,szl,taul,&
     densr,sxr,syr,szr,taur,flux1,flux2,flux3,flux4,flux5)

  USE GRHydro_Scalars

  implicit none

  REAL*8 GRHydro_Y_e_max
  integer, parameter :: cctki_use_GRHydro_Y_e_max = kind(GRHydro_Y_e_max)
  REAL*8 GRHydro_Y_e_min
  integer, parameter :: cctki_use_GRHydro_Y_e_min = kind(GRHydro_Y_e_min)
  REAL*8 GRHydro_atmo_tolerance
  integer, parameter :: cctki_use_GRHydro_atmo_tolerance = kind(GRHydro_atmo_tolerance)
  REAL*8 GRHydro_del_ptol
  integer, parameter :: cctki_use_GRHydro_del_ptol = kind(GRHydro_del_ptol)
  REAL*8 GRHydro_eos_rf_prec
  integer, parameter :: cctki_use_GRHydro_eos_rf_prec = kind(GRHydro_eos_rf_prec)
  REAL*8 GRHydro_eps_min
  integer, parameter :: cctki_use_GRHydro_eps_min = kind(GRHydro_eps_min)
  REAL*8 GRHydro_hot_atmo_Y_e
  integer, parameter :: cctki_use_GRHydro_hot_atmo_Y_e = kind(GRHydro_hot_atmo_Y_e)
  REAL*8 GRHydro_hot_atmo_temp
  integer, parameter :: cctki_use_GRHydro_hot_atmo_temp = kind(GRHydro_hot_atmo_temp)
  REAL*8 GRHydro_lorentz_overshoot_cutoff
  integer, parameter :: cctki_use_GRHydro_lorentz_overshoot_cutoff = kind(GRHydro_lorentz_overshoot_cutoff)
  REAL*8 GRHydro_max_temp
  integer, parameter :: cctki_use_GRHydro_max_temp = kind(GRHydro_max_temp)
  REAL*8 GRHydro_perc_ptol
  integer, parameter :: cctki_use_GRHydro_perc_ptol = kind(GRHydro_perc_ptol)
  REAL*8 GRHydro_rho_central
  integer, parameter :: cctki_use_GRHydro_rho_central = kind(GRHydro_rho_central)
  REAL*8 Tmunu_damping_radius_max
  integer, parameter :: cctki_use_Tmunu_damping_radius_max = kind(Tmunu_damping_radius_max)
  REAL*8 Tmunu_damping_radius_min
  integer, parameter :: cctki_use_Tmunu_damping_radius_min = kind(Tmunu_damping_radius_min)
  REAL*8 atmo_falloff_power
  integer, parameter :: cctki_use_atmo_falloff_power = kind(atmo_falloff_power)
  REAL*8 atmo_falloff_radius
  integer, parameter :: cctki_use_atmo_falloff_radius = kind(atmo_falloff_radius)
  REAL*8 atmo_tolerance_power
  integer, parameter :: cctki_use_atmo_tolerance_power = kind(atmo_tolerance_power)
  REAL*8 atmo_tolerance_radius
  integer, parameter :: cctki_use_atmo_tolerance_radius = kind(atmo_tolerance_radius)
  REAL*8 c2p_reset_pressure_to_value
  integer, parameter :: cctki_use_c2p_reset_pressure_to_value = kind(c2p_reset_pressure_to_value)
  REAL*8 comoving_factor
  integer, parameter :: cctki_use_comoving_factor = kind(comoving_factor)
  REAL*8 comoving_tanh_factor
  integer, parameter :: cctki_use_comoving_tanh_factor = kind(comoving_tanh_factor)
  REAL*8 comoving_tanh_offset
  integer, parameter :: cctki_use_comoving_tanh_offset = kind(comoving_tanh_offset)
  REAL*8 enhanced_ppm_C2
  integer, parameter :: cctki_use_enhanced_ppm_C2 = kind(enhanced_ppm_C2)
  REAL*8 enhanced_ppm_C3
  integer, parameter :: cctki_use_enhanced_ppm_C3 = kind(enhanced_ppm_C3)
  REAL*8 initial_Gamma
  integer, parameter :: cctki_use_initial_Gamma = kind(initial_Gamma)
  REAL*8 initial_atmosphere_factor
  integer, parameter :: cctki_use_initial_atmosphere_factor = kind(initial_atmosphere_factor)
  REAL*8 initial_k
  integer, parameter :: cctki_use_initial_k = kind(initial_k)
  REAL*8 initial_rho_abs_min
  integer, parameter :: cctki_use_initial_rho_abs_min = kind(initial_rho_abs_min)
  REAL*8 initial_rho_rel_min
  integer, parameter :: cctki_use_initial_rho_rel_min = kind(initial_rho_rel_min)
  REAL*8 kap_dc
  integer, parameter :: cctki_use_kap_dc = kind(kap_dc)
  REAL*8 max_magnetic_to_gas_pressure_ratio
  integer, parameter :: cctki_use_max_magnetic_to_gas_pressure_ratio = kind(max_magnetic_to_gas_pressure_ratio)
  REAL*8 min_tracer
  integer, parameter :: cctki_use_min_tracer = kind(min_tracer)
  REAL*8 mp5_alpha
  integer, parameter :: cctki_use_mp5_alpha = kind(mp5_alpha)
  REAL*8 mp5_eps
  integer, parameter :: cctki_use_mp5_eps = kind(mp5_eps)
  REAL*8 myfloor
  integer, parameter :: cctki_use_myfloor = kind(myfloor)
  REAL*8 ppm_epsilon
  integer, parameter :: cctki_use_ppm_epsilon = kind(ppm_epsilon)
  REAL*8 ppm_epsilon_shock
  integer, parameter :: cctki_use_ppm_epsilon_shock = kind(ppm_epsilon_shock)
  REAL*8 ppm_eta1
  integer, parameter :: cctki_use_ppm_eta1 = kind(ppm_eta1)
  REAL*8 ppm_eta2
  integer, parameter :: cctki_use_ppm_eta2 = kind(ppm_eta2)
  REAL*8 ppm_k0
  integer, parameter :: cctki_use_ppm_k0 = kind(ppm_k0)
  REAL*8 ppm_omega1
  integer, parameter :: cctki_use_ppm_omega1 = kind(ppm_omega1)
  REAL*8 ppm_omega2
  integer, parameter :: cctki_use_ppm_omega2 = kind(ppm_omega2)
  REAL*8 ppm_omega_tracer
  integer, parameter :: cctki_use_ppm_omega_tracer = kind(ppm_omega_tracer)
  REAL*8 ppm_small
  integer, parameter :: cctki_use_ppm_small = kind(ppm_small)
  REAL*8 rho_abs_min
  integer, parameter :: cctki_use_rho_abs_min = kind(rho_abs_min)
  REAL*8 rho_abs_min_after_recovery
  integer, parameter :: cctki_use_rho_abs_min_after_recovery = kind(rho_abs_min_after_recovery)
  REAL*8 rho_rel_min
  integer, parameter :: cctki_use_rho_rel_min = kind(rho_rel_min)
  REAL*8 setcharmax
  integer, parameter :: cctki_use_setcharmax = kind(setcharmax)
  REAL*8 setcharmin
  integer, parameter :: cctki_use_setcharmin = kind(setcharmin)
  REAL*8 sqrtdet_thr
  integer, parameter :: cctki_use_sqrtdet_thr = kind(sqrtdet_thr)
  REAL*8 tau_rel_min
  integer, parameter :: cctki_use_tau_rel_min = kind(tau_rel_min)
  REAL*8 weno_eps
  integer, parameter :: cctki_use_weno_eps = kind(weno_eps)
  integer*8 Avec_gauge
  integer, parameter :: cctki_use_Avec_gauge = kind(Avec_gauge)
  integer*8 EoS_Change_type
  integer, parameter :: cctki_use_EoS_Change_type = kind(EoS_Change_type)
  integer*8 GRHydro_c2p_failed_action
  integer, parameter :: cctki_use_GRHydro_c2p_failed_action = kind(GRHydro_c2p_failed_action)
  integer*8 GRHydro_eos_table
  integer, parameter :: cctki_use_GRHydro_eos_table = kind(GRHydro_eos_table)
  integer*8 GRHydro_eos_type
  integer, parameter :: cctki_use_GRHydro_eos_type = kind(GRHydro_eos_type)
  integer*8 HLLE_type
  integer, parameter :: cctki_use_HLLE_type = kind(HLLE_type)
  integer*8 bound
  integer, parameter :: cctki_use_bound = kind(bound)
  integer*8 comoving_attenuate
  integer, parameter :: cctki_use_comoving_attenuate = kind(comoving_attenuate)
  integer*8 comoving_v_method
  integer, parameter :: cctki_use_comoving_v_method = kind(comoving_v_method)
  integer*8 gradient_method
  integer, parameter :: cctki_use_gradient_method = kind(gradient_method)
  integer*8 left_eigenvectors
  integer, parameter :: cctki_use_left_eigenvectors = kind(left_eigenvectors)
  integer*8 method_type
  integer, parameter :: cctki_use_method_type = kind(method_type)
  integer*8 numerical_viscosity
  integer, parameter :: cctki_use_numerical_viscosity = kind(numerical_viscosity)
  integer*8 particle_interpolator
  integer, parameter :: cctki_use_particle_interpolator = kind(particle_interpolator)
  integer*8 ppm_flatten
  integer, parameter :: cctki_use_ppm_flatten = kind(ppm_flatten)
  integer*8 psidcspeed
  integer, parameter :: cctki_use_psidcspeed = kind(psidcspeed)
  integer*8 recon_method
  integer, parameter :: cctki_use_recon_method = kind(recon_method)
  integer*8 recon_vars
  integer, parameter :: cctki_use_recon_vars = kind(recon_vars)
  integer*8 riemann_solver
  integer, parameter :: cctki_use_riemann_solver = kind(riemann_solver)
  integer*8 tvd_limiter
  integer, parameter :: cctki_use_tvd_limiter = kind(tvd_limiter)
  integer*8 use_evolution_mask
  integer, parameter :: cctki_use_use_evolution_mask = kind(use_evolution_mask)
  INTEGER*4 Check_Rho_Minimum
  integer, parameter :: cctki_use_Check_Rho_Minimum = kind(Check_Rho_Minimum)
  INTEGER*4 EoS_Change
  integer, parameter :: cctki_use_EoS_Change = kind(EoS_Change)
  INTEGER*4 GRHydro_MaxNumConstrainedVars
  integer, parameter :: cctki_use_GRHydro_MaxNumConstrainedVars = kind(GRHydro_MaxNumConstrainedVars)
  INTEGER*4 GRHydro_MaxNumEvolvedVars
  integer, parameter :: cctki_use_GRHydro_MaxNumEvolvedVars = kind(GRHydro_MaxNumEvolvedVars)
  INTEGER*4 GRHydro_MaxNumEvolvedVarsSlow
  integer, parameter :: cctki_use_GRHydro_MaxNumEvolvedVarsSlow = kind(GRHydro_MaxNumEvolvedVarsSlow)
  INTEGER*4 GRHydro_MaxNumSandRVars
  integer, parameter :: cctki_use_GRHydro_MaxNumSandRVars = kind(GRHydro_MaxNumSandRVars)
  INTEGER*4 GRHydro_NaN_verbose
  integer, parameter :: cctki_use_GRHydro_NaN_verbose = kind(GRHydro_NaN_verbose)
  INTEGER*4 GRHydro_c2p_reset_eps_tau_hot_eos
  integer, parameter :: cctki_use_GRHydro_c2p_reset_eps_tau_hot_eos = kind(GRHydro_c2p_reset_eps_tau_hot_eos)
  INTEGER*4 GRHydro_c2p_warn_from_reflevel
  integer, parameter :: cctki_use_GRHydro_c2p_warn_from_reflevel = kind(GRHydro_c2p_warn_from_reflevel)
  INTEGER*4 GRHydro_c2p_warnlevel
  integer, parameter :: cctki_use_GRHydro_c2p_warnlevel = kind(GRHydro_c2p_warnlevel)
  INTEGER*4 GRHydro_countmax
  integer, parameter :: cctki_use_GRHydro_countmax = kind(GRHydro_countmax)
  INTEGER*4 GRHydro_countmin
  integer, parameter :: cctki_use_GRHydro_countmin = kind(GRHydro_countmin)
  INTEGER*4 GRHydro_enable_internal_excision
  integer, parameter :: cctki_use_GRHydro_enable_internal_excision = kind(GRHydro_enable_internal_excision)
  INTEGER*4 GRHydro_eos_hot_eps_fix
  integer, parameter :: cctki_use_GRHydro_eos_hot_eps_fix = kind(GRHydro_eos_hot_eps_fix)
  INTEGER*4 GRHydro_eos_hot_prim2con_warn
  integer, parameter :: cctki_use_GRHydro_eos_hot_prim2con_warn = kind(GRHydro_eos_hot_prim2con_warn)
  INTEGER*4 GRHydro_hydro_excision
  integer, parameter :: cctki_use_GRHydro_hydro_excision = kind(GRHydro_hydro_excision)
  INTEGER*4 GRHydro_oppm_reflevel
  integer, parameter :: cctki_use_GRHydro_oppm_reflevel = kind(GRHydro_oppm_reflevel)
  INTEGER*4 GRHydro_polish
  integer, parameter :: cctki_use_GRHydro_polish = kind(GRHydro_polish)
  INTEGER*4 GRHydro_stencil
  integer, parameter :: cctki_use_GRHydro_stencil = kind(GRHydro_stencil)
  INTEGER*4 WENO_order
  integer, parameter :: cctki_use_WENO_order = kind(WENO_order)
  INTEGER*4 apply_H_viscosity
  integer, parameter :: cctki_use_apply_H_viscosity = kind(apply_H_viscosity)
  INTEGER*4 c2p_reset_pressure
  integer, parameter :: cctki_use_c2p_reset_pressure = kind(c2p_reset_pressure)
  INTEGER*4 c2p_resort_to_bisection
  integer, parameter :: cctki_use_c2p_resort_to_bisection = kind(c2p_resort_to_bisection)
  INTEGER*4 calculate_bcom
  integer, parameter :: cctki_use_calculate_bcom = kind(calculate_bcom)
  INTEGER*4 check_for_trivial_rp
  integer, parameter :: cctki_use_check_for_trivial_rp = kind(check_for_trivial_rp)
  INTEGER*4 clean_divergence
  integer, parameter :: cctki_use_clean_divergence = kind(clean_divergence)
  INTEGER*4 con2prim_oct_hack
  integer, parameter :: cctki_use_con2prim_oct_hack = kind(con2prim_oct_hack)
  INTEGER*4 decouple_normal_Bfield
  integer, parameter :: cctki_use_decouple_normal_Bfield = kind(decouple_normal_Bfield)
  INTEGER*4 eno_order
  integer, parameter :: cctki_use_eno_order = kind(eno_order)
  INTEGER*4 evolve_tracer
  integer, parameter :: cctki_use_evolve_tracer = kind(evolve_tracer)
  INTEGER*4 mp5_adaptive_eps
  integer, parameter :: cctki_use_mp5_adaptive_eps = kind(mp5_adaptive_eps)
  INTEGER*4 number_of_arrays
  integer, parameter :: cctki_use_number_of_arrays = kind(number_of_arrays)
  INTEGER*4 number_of_particles
  integer, parameter :: cctki_use_number_of_particles = kind(number_of_particles)
  INTEGER*4 number_of_tracers
  integer, parameter :: cctki_use_number_of_tracers = kind(number_of_tracers)
  INTEGER*4 particle_interpolation_order
  integer, parameter :: cctki_use_particle_interpolation_order = kind(particle_interpolation_order)
  INTEGER*4 ppm_detect
  integer, parameter :: cctki_use_ppm_detect = kind(ppm_detect)
  INTEGER*4 ppm_mppm
  integer, parameter :: cctki_use_ppm_mppm = kind(ppm_mppm)
  INTEGER*4 ppm_mppm_debug_eigenvalues
  integer, parameter :: cctki_use_ppm_mppm_debug_eigenvalues = kind(ppm_mppm_debug_eigenvalues)
  INTEGER*4 reconstruct_Wv
  integer, parameter :: cctki_use_reconstruct_Wv = kind(reconstruct_Wv)
  INTEGER*4 reconstruct_temper
  integer, parameter :: cctki_use_reconstruct_temper = kind(reconstruct_temper)
  INTEGER*4 set_trivial_rp_grid_function
  integer, parameter :: cctki_use_set_trivial_rp_grid_function = kind(set_trivial_rp_grid_function)
  INTEGER*4 sources_spatial_order
  integer, parameter :: cctki_use_sources_spatial_order = kind(sources_spatial_order)
  INTEGER*4 sync_conserved_only
  integer, parameter :: cctki_use_sync_conserved_only = kind(sync_conserved_only)
  INTEGER*4 track_divB
  integer, parameter :: cctki_use_track_divB = kind(track_divB)
  INTEGER*4 transport_constraints
  integer, parameter :: cctki_use_transport_constraints = kind(transport_constraints)
  INTEGER*4 use_MoL_slow_multirate_sector
  integer, parameter :: cctki_use_use_MoL_slow_multirate_sector = kind(use_MoL_slow_multirate_sector)
  INTEGER*4 use_enhanced_ppm
  integer, parameter :: cctki_use_use_enhanced_ppm = kind(use_enhanced_ppm)
  INTEGER*4 use_min_tracer
  integer, parameter :: cctki_use_use_min_tracer = kind(use_min_tracer)
  INTEGER*4 use_optimized_ppm
  integer, parameter :: cctki_use_use_optimized_ppm = kind(use_optimized_ppm)
  INTEGER*4 use_weighted_fluxes
  integer, parameter :: cctki_use_use_weighted_fluxes = kind(use_weighted_fluxes)
  INTEGER*4 weno_adaptive_epsilon
  integer, parameter :: cctki_use_weno_adaptive_epsilon = kind(weno_adaptive_epsilon)
  INTEGER*4 wk_atmosphere
  integer, parameter :: cctki_use_wk_atmosphere = kind(wk_atmosphere)
  COMMON /GRHydrorest/GRHydro_Y_e_max, GRHydro_Y_e_min, GRHydro_atmo_tolerance, GRHydro_del_ptol, GRHydro_eos_rf_prec, GRHydro_eps_&
  &min, GRHydro_hot_atmo_Y_e, GRHydro_hot_atmo_temp, GRHydro_lorentz_overshoot_cutoff, GRHydro_max_temp, GRHydro_perc_ptol, GRHydro&
  &_rho_central, Tmunu_damping_radius_max, Tmunu_damping_radius_min, atmo_falloff_power, atmo_falloff_radius, atmo_tolerance_power,&
  & atmo_tolerance_radius, c2p_reset_pressure_to_value, comoving_factor, comoving_tanh_factor, comoving_tanh_offset, enhanced_ppm_C&
  &2, enhanced_ppm_C3, initial_Gamma, initial_atmosphere_factor, initial_k, initial_rho_abs_min, initial_rho_rel_min, kap_dc, max_m&
  &agnetic_to_gas_pressure_ratio, min_tracer, mp5_alpha, mp5_eps, myfloor, ppm_epsilon, ppm_epsilon_shock, ppm_eta1, ppm_eta2, ppm_&
  &k0, ppm_omega1, ppm_omega2, ppm_omega_tracer, ppm_small, rho_abs_min, rho_abs_min_after_recovery, rho_rel_min, setcharmax, setch&
  &armin, sqrtdet_thr, tau_rel_min, weno_eps, Avec_gauge, EoS_Change_type, GRHydro_c2p_failed_action, GRHydro_eos_table, GRHydro_eo&
  &s_type, HLLE_type, bound, comoving_attenuate, comoving_v_method, gradient_method, left_eigenvectors, method_type, numerical_visc&
  &osity, particle_interpolator, ppm_flatten, psidcspeed, recon_method, recon_vars, riemann_solver, tvd_limiter, use_evolution_mask&
  &, Check_Rho_Minimum, EoS_Change, GRHydro_MaxNumConstrainedVars, GRHydro_MaxNumEvolvedVars, GRHydro_MaxNumEvolvedVarsSlow, GRHydr&
  &o_MaxNumSandRVars, GRHydro_NaN_verbose, GRHydro_c2p_reset_eps_tau_hot_eos, GRHydro_c2p_warn_from_reflevel, GRHydro_c2p_warnlevel&
  &, GRHydro_countmax, GRHydro_countmin, GRHydro_enable_internal_excision, GRHydro_eos_hot_eps_fix, GRHydro_eos_hot_prim2con_warn, &
  &GRHydro_hydro_excision, GRHydro_oppm_reflevel, GRHydro_polish, GRHydro_stencil, WENO_order, apply_H_viscosity, c2p_reset_pressur&
  &e, c2p_resort_to_bisection, calculate_bcom, check_for_trivial_rp, clean_divergence, con2prim_oct_hack, decouple_normal_Bfield, e&
  &no_order, evolve_tracer, mp5_adaptive_eps, number_of_arrays, number_of_particles, number_of_tracers, particle_interpolation_orde&
  &r, ppm_detect, ppm_mppm, ppm_mppm_debug_eigenvalues, reconstruct_Wv, reconstruct_temper, set_trivial_rp_grid_function, sources_s&
  &patial_order, sync_conserved_only, track_divB, transport_constraints, use_MoL_slow_multirate_sector, use_enhanced_ppm, use_min_t&
  &racer, use_optimized_ppm, use_weighted_fluxes, weno_adaptive_epsilon, wk_atmosphere
  INTEGER*4 constrain_to_1D
  integer, parameter :: cctki_use_constrain_to_1D = kind(constrain_to_1D)
  INTEGER*4 use_cxx_code
  integer, parameter :: cctki_use_use_cxx_code = kind(use_cxx_code)
  INTEGER*4 verbose
  integer, parameter :: cctki_use_verbose = kind(verbose)
  COMMON /GRHydropriv/constrain_to_1D, use_cxx_code, verbose
  integer*8 CCTKH0
  integer, parameter :: cctki_use_CCTKH0 = kind(CCTKH0)
  integer*8 CCTKH1
  integer, parameter :: cctki_use_CCTKH1 = kind(CCTKH1)
  integer*8 CCTKH2
  integer, parameter :: cctki_use_CCTKH2 = kind(CCTKH2)
  integer*8 CCTKH3
  integer, parameter :: cctki_use_CCTKH3 = kind(CCTKH3)
  integer*8 CCTKH4
  integer, parameter :: cctki_use_CCTKH4 = kind(CCTKH4)
  integer*8 CCTKH5
  integer, parameter :: cctki_use_CCTKH5 = kind(CCTKH5)
  integer*8 CCTKH6
  integer, parameter :: cctki_use_CCTKH6 = kind(CCTKH6)
  integer*8 CCTKH7
  integer, parameter :: cctki_use_CCTKH7 = kind(CCTKH7)
  integer*8 initial_shift
  integer, parameter :: cctki_use_initial_shift = kind(initial_shift)
  integer*8 CCTKH8
  integer, parameter :: cctki_use_CCTKH8 = kind(CCTKH8)
  integer*8 CCTKH9
  integer, parameter :: cctki_use_CCTKH9 = kind(CCTKH9)
  integer*8 CCTKH11
  integer, parameter :: cctki_use_CCTKH11 = kind(CCTKH11)
  integer*8 metric_type
  integer, parameter :: cctki_use_metric_type = kind(metric_type)
  integer*8 shift_evolution_method
  integer, parameter :: cctki_use_shift_evolution_method = kind(shift_evolution_method)
  integer*8 CCTKH13
  integer, parameter :: cctki_use_CCTKH13 = kind(CCTKH13)
  INTEGER*4 CCTKH10
  integer, parameter :: cctki_use_CCTKH10 = kind(CCTKH10)
  INTEGER*4 CCTKH12
  integer, parameter :: cctki_use_CCTKH12 = kind(CCTKH12)
  INTEGER*4 CCTKH14
  integer, parameter :: cctki_use_CCTKH14 = kind(CCTKH14)
  COMMON /ADMBASErest/CCTKH0, CCTKH1, CCTKH2, CCTKH3, CCTKH4, CCTKH5, CCTKH6, CCTKH7, initial_shift, CCTKH8, CCTKH9, CCTKH11, metri&
  &c_type, shift_evolution_method, CCTKH13, CCTKH10, CCTKH12, CCTKH14
  REAL*8 CCTKH15
  integer, parameter :: cctki_use_CCTKH15 = kind(CCTKH15)
  REAL*8 CCTKH25
  integer, parameter :: cctki_use_CCTKH25 = kind(CCTKH25)
  REAL*8 CCTKH26
  integer, parameter :: cctki_use_CCTKH26 = kind(CCTKH26)
  REAL*8 CCTKH27 (10)
  integer, parameter :: cctki_use_CCTKH27 = kind(CCTKH27)
  REAL*8 CCTKH28
  integer, parameter :: cctki_use_CCTKH28 = kind(CCTKH28)
  REAL*8 CCTKH29
  integer, parameter :: cctki_use_CCTKH29 = kind(CCTKH29)
  REAL*8 CCTKH30
  integer, parameter :: cctki_use_CCTKH30 = kind(CCTKH30)
  REAL*8 CCTKH31
  integer, parameter :: cctki_use_CCTKH31 = kind(CCTKH31)
  REAL*8 CCTKH32
  integer, parameter :: cctki_use_CCTKH32 = kind(CCTKH32)
  REAL*8 CCTKH33 (9)
  integer, parameter :: cctki_use_CCTKH33 = kind(CCTKH33)
  REAL*8 CCTKH34
  integer, parameter :: cctki_use_CCTKH34 = kind(CCTKH34)
  REAL*8 CCTKH38
  integer, parameter :: cctki_use_CCTKH38 = kind(CCTKH38)
  REAL*8 CCTKH39
  integer, parameter :: cctki_use_CCTKH39 = kind(CCTKH39)
  REAL*8 poly_k
  integer, parameter :: cctki_use_poly_k = kind(poly_k)
  integer*8 CCTKH17
  integer, parameter :: cctki_use_CCTKH17 = kind(CCTKH17)
  integer*8 CCTKH20
  integer, parameter :: cctki_use_CCTKH20 = kind(CCTKH20)
  integer*8 CCTKH24
  integer, parameter :: cctki_use_CCTKH24 = kind(CCTKH24)
  integer*8 CCTKH37
  integer, parameter :: cctki_use_CCTKH37 = kind(CCTKH37)
  INTEGER*4 CCTKH16
  integer, parameter :: cctki_use_CCTKH16 = kind(CCTKH16)
  INTEGER*4 CCTKH18
  integer, parameter :: cctki_use_CCTKH18 = kind(CCTKH18)
  INTEGER*4 CCTKH19
  integer, parameter :: cctki_use_CCTKH19 = kind(CCTKH19)
  INTEGER*4 CCTKH21
  integer, parameter :: cctki_use_CCTKH21 = kind(CCTKH21)
  INTEGER*4 CCTKH22
  integer, parameter :: cctki_use_CCTKH22 = kind(CCTKH22)
  INTEGER*4 CCTKH23
  integer, parameter :: cctki_use_CCTKH23 = kind(CCTKH23)
  INTEGER*4 CCTKH35
  integer, parameter :: cctki_use_CCTKH35 = kind(CCTKH35)
  INTEGER*4 CCTKH36
  integer, parameter :: cctki_use_CCTKH36 = kind(CCTKH36)
  INTEGER*4 CCTKH40
  integer, parameter :: cctki_use_CCTKH40 = kind(CCTKH40)
  INTEGER*4 CCTKH41
  integer, parameter :: cctki_use_CCTKH41 = kind(CCTKH41)
  COMMON /EOS_OMNIrest/CCTKH15, CCTKH25, CCTKH26, CCTKH27, CCTKH28, CCTKH29, CCTKH30, CCTKH31, CCTKH32, CCTKH33, CCTKH34, CCTKH38, &
  &CCTKH39, poly_k, CCTKH17, CCTKH20, CCTKH24, CCTKH37, CCTKH16, CCTKH18, CCTKH19, CCTKH21, CCTKH22, CCTKH23, CCTKH35, CCTKH36, CCT&
  &KH40, CCTKH41
  integer*8 CCTKH42
  integer, parameter :: cctki_use_CCTKH42 = kind(CCTKH42)
  integer*8 Bvec_evolution_method
  integer, parameter :: cctki_use_Bvec_evolution_method = kind(Bvec_evolution_method)
  integer*8 Y_e_evolution_method
  integer, parameter :: cctki_use_Y_e_evolution_method = kind(Y_e_evolution_method)
  integer*8 entropy_evolution_method
  integer, parameter :: cctki_use_entropy_evolution_method = kind(entropy_evolution_method)
  integer*8 evolution_method
  integer, parameter :: cctki_use_evolution_method = kind(evolution_method)
  integer*8 CCTKH43
  integer, parameter :: cctki_use_CCTKH43 = kind(CCTKH43)
  integer*8 CCTKH44
  integer, parameter :: cctki_use_CCTKH44 = kind(CCTKH44)
  integer*8 CCTKH45
  integer, parameter :: cctki_use_CCTKH45 = kind(CCTKH45)
  integer*8 initial_Bvec
  integer, parameter :: cctki_use_initial_Bvec = kind(initial_Bvec)
  integer*8 CCTKH46
  integer, parameter :: cctki_use_CCTKH46 = kind(CCTKH46)
  integer*8 CCTKH47
  integer, parameter :: cctki_use_CCTKH47 = kind(CCTKH47)
  integer*8 CCTKH48
  integer, parameter :: cctki_use_CCTKH48 = kind(CCTKH48)
  integer*8 CCTKH49
  integer, parameter :: cctki_use_CCTKH49 = kind(CCTKH49)
  integer*8 prolongation_type
  integer, parameter :: cctki_use_prolongation_type = kind(prolongation_type)
  integer*8 temperature_evolution_method
  integer, parameter :: cctki_use_temperature_evolution_method = kind(temperature_evolution_method)
  INTEGER*4 hydro_excision
  integer, parameter :: cctki_use_hydro_excision = kind(hydro_excision)
  INTEGER*4 timelevels
  integer, parameter :: cctki_use_timelevels = kind(timelevels)
  COMMON /HYDROBASErest/CCTKH42, Bvec_evolution_method, Y_e_evolution_method, entropy_evolution_method, evolution_method, CCTKH43, &
  &CCTKH44, CCTKH45, initial_Bvec, CCTKH46, CCTKH47, CCTKH48, CCTKH49, prolongation_type, temperature_evolution_method, hydro_excis&
  &ion, timelevels
  INTEGER*4 MoL_Max_Evolved_Array_Size
  integer, parameter :: cctki_use_MoL_Max_Evolved_Array_Size = kind(MoL_Max_Evolved_Array_Size)
  INTEGER*4 CCTKH50
  integer, parameter :: cctki_use_CCTKH50 = kind(CCTKH50)
  INTEGER*4 MoL_Num_ArrayEvolved_Vars
  integer, parameter :: cctki_use_MoL_Num_ArrayEvolved_Vars = kind(MoL_Num_ArrayEvolved_Vars)
  INTEGER*4 CCTKH51
  integer, parameter :: cctki_use_CCTKH51 = kind(CCTKH51)
  INTEGER*4 MoL_Num_Constrained_Vars
  integer, parameter :: cctki_use_MoL_Num_Constrained_Vars = kind(MoL_Num_Constrained_Vars)
  INTEGER*4 MoL_Num_Evolved_Vars
  integer, parameter :: cctki_use_MoL_Num_Evolved_Vars = kind(MoL_Num_Evolved_Vars)
  INTEGER*4 MoL_Num_Evolved_Vars_Slow
  integer, parameter :: cctki_use_MoL_Num_Evolved_Vars_Slow = kind(MoL_Num_Evolved_Vars_Slow)
  INTEGER*4 MoL_Num_SaveAndRestore_Vars
  integer, parameter :: cctki_use_MoL_Num_SaveAndRestore_Vars = kind(MoL_Num_SaveAndRestore_Vars)
  INTEGER*4 CCTKH52
  integer, parameter :: cctki_use_CCTKH52 = kind(CCTKH52)
  COMMON /METHODOFLINESrest/MoL_Max_Evolved_Array_Size, CCTKH50, MoL_Num_ArrayEvolved_Vars, CCTKH51, MoL_Num_Constrained_Vars, MoL_&
  &Num_Evolved_Vars, MoL_Num_Evolved_Vars_Slow, MoL_Num_SaveAndRestore_Vars, CCTKH52
  INTEGER*4 use_mask
  integer, parameter :: cctki_use_use_mask = kind(use_mask)
  COMMON /SPACEMASKrest/use_mask
  

  REAL*8 rhor,velxr,velyr,velzr,epsr
  REAL*8 rhol,velxl,velyl,velzl,epsl
  REAL*8 densl,sxl,syl,szl,taul
  REAL*8 densr,sxr,syr,szr,taur

  REAL*8 lam(5),p(5,5),q(5,5),dw(5)
  REAL*8 rfluxr(5),rfluxl(5)
  REAL*8 gxx,gxy,gxz,gyy,gyz,gzz,u,det
  REAL*8 alp,beta,flux1r,flux2r,flux3r,flux4r,flux5r
  REAL*8 flux1l,flux2l,flux3l,flux4l,flux5l
  REAL*8 flux1,flux2,flux3,flux4,flux5



  REAL*8 du(5),aa(5,5)
  integer i,j,k
  REAL*8 paug(5,6),tmp1,tmp2,sump,summ,f_du(5)
  REAL*8 leivec1l(5),leivec2l(5),leivec3l(5),leivecpl(5),leivecml(5)
  REAL*8 leivec1r(5),leivec2r(5),leivec3r(5),leivecpr(5),leivecmr(5)
  REAL*8 reivec1l(5),reivec2l(5),reivec3l(5),reivecpl(5),reivecml(5)
  REAL*8 reivec1r(5),reivec2r(5),reivec3r(5),reivecpr(5),reivecmr(5)
  REAL*8 lam1l,lam2l,lam3l,lamml,lampl,lamm_nobetal,lamp_nobetal
  REAL*8 lam1r,lam2r,lam3r,lammr,lampr,lamm_nobetar,lamp_nobetar
  REAL*8 lam1,lam2,lam3,lamm,lamp,lamp_tmp

  REAL*8 cs2l,cs2r,one,two
  REAL*8 vlowxr,vlowyr,vlowzr,v2r,wr
  REAL*8 vlowxl,vlowyl,vlowzl,v2l,wl

  REAL*8 pressl,dpdrhol,dpdepsl,enthalpyl,kappal
  REAL*8 pressr,dpdrhor,dpdepsr,enthalpyr,kappar
  REAL*8 axpl,axml,vxpl,vxml,xsil,dltl
  REAL*8 axpr,axmr,vxpr,vxmr,xsir,dltr
  REAL*8 cxx,cxy,cxz,gam,vxa,vxb
  INTEGER*4 handle
  REAL*8 invrhol, betainvalp, invfacl
  REAL*8 invrhor, invfacr
  REAL*8 invulpfac, invulmfac, invurpfac, invurmfac


  INTEGER*4  :: n 
  INTEGER*4  :: keytemp 
  INTEGER*4  :: anyerr 
  INTEGER*4  :: keyerr(1) 
  REAL*8 :: xpress 
  REAL*8 :: xeps 
  REAL*8 :: xtemp
  REAL*8 :: xye 

  n=1;keytemp=0;anyerr=0;keyerr(1)=0
  xpress=0.0d0;xeps=0.0d0;xtemp=0.0d0;xye=0.0d0

  one = 1.0d0
  two = 2.0d0




  invrhol = one / rhol

  call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
       rhol,epsl,xtemp,xye,pressl,keyerr,anyerr)


  call EOS_Omni_DPressByDEps(handle,keytemp,GRHydro_eos_rf_prec,n,&
       rhol,epsl,xtemp,xye,dpdepsl,keyerr,anyerr)
  call EOS_Omni_DPressByDRho(handle,keytemp,GRHydro_eos_rf_prec,n,&
       rhol,epsl,xtemp,xye,dpdrhol,keyerr,anyerr)
  cs2l = (dpdrhol + pressl * dpdepsl * invrhol**2)/ &
       (1.0d0 + epsl + pressl*invrhol)

  enthalpyl = one + epsl + pressl * invrhol

  vlowxl = gxx*velxl + gxy*velyl + gxz*velzl
  vlowyl = gxy*velxl + gyy*velyl + gyz*velzl
  vlowzl = gxz*velxl + gyz*velyl + gzz*velzl
  v2l = vlowxl*velxl + vlowyl*velyl + vlowzl*velzl

  if (v2l .ge. 1.d0-1.0d-10) then
    v2l = 1.d0-1.0d-10
  endif


  wl = one / sqrt(one - v2l)



  betainvalp = beta / alp
  lam1l = velxl - betainvalp
  lam2l = lam1l
  lam3l = lam1l
  lamp_tmp = cs2l*(one-v2l)*(u*(one-v2l*cs2l) - velxl**2*(one-cs2l))
  if (lamp_tmp .le. 1.0d-10) then
    lamp_tmp = 1.0d-10;
  endif
  invfacl = one/(one-v2l*cs2l)
  lamp_nobetal = (velxl*(one-cs2l) + sqrt(lamp_tmp))*invfacl
  lamm_nobetal = (velxl*(one-cs2l) - sqrt(lamp_tmp))*invfacl

  lampl = lamp_nobetal - betainvalp
  lamml = lamm_nobetal - betainvalp




  invrhor = one / rhor

  call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
       rhor,epsr,xtemp,xye,pressr,keyerr,anyerr)


  call EOS_Omni_DPressByDEps(handle,keytemp,GRHydro_eos_rf_prec,n,&
       rhor,epsr,xtemp,xye,dpdepsr,keyerr,anyerr)
  call EOS_Omni_DPressByDRho(handle,keytemp,GRHydro_eos_rf_prec,n,&
       rhor,epsr,xtemp,xye,dpdrhor,keyerr,anyerr)
  cs2r = (dpdrhor + pressr * dpdepsr * invrhor**2)/ &
       (1.0d0 + epsr + pressr*invrhor)

  enthalpyr = one + epsr + pressr * invrhor

  vlowxr = gxx*velxr + gxy*velyr + gxz*velzr
  vlowyr = gxy*velxr + gyy*velyr + gyz*velzr
  vlowzr = gxz*velxr + gyz*velyr + gzz*velzr
  v2r = vlowxr*velxr + vlowyr*velyr + vlowzr*velzr

  if (v2r .ge. 1.d0-1.0d-10) then
    v2r = 1.d0-1.0d-10
  endif

  wr = one / sqrt(one - v2r)



  lam1r = velxr - betainvalp
  lam2r = lam1r
  lam3r = lam1r
  lamp_tmp = cs2r*(one-v2r)*(u*(one-v2r*cs2r) - velxr**2*(one-cs2r))
  if (lamp_tmp .le. 1.0d-10) then
    lamp_tmp = 1.0d-10;
  endif
  invfacr = one/(one-v2r*cs2r)
  lamp_nobetar = (velxr*(one-cs2r) + sqrt(lamp_tmp))*invfacr
  lamm_nobetar = (velxr*(one-cs2r) - sqrt(lamp_tmp))*invfacr

  lampr = lamp_nobetar - betainvalp
  lammr = lamm_nobetar - betainvalp



  lam1 = max(abs(lam1l),abs(lam1r))
  lam2 = lam1
  lam3 = lam1
  lamp = max(abs(lampl),abs(lampr))
  lamm = max(abs(lamml),abs(lammr))







  lam(1) = lamm
  lam(5) = lam1
  lam(3) = lam2
  lam(4) = lam3
  lam(2) = lamp
  




  invulpfac = one / (u - velxl*lamp_nobetal)
  invulmfac = one / (u - velxl*lamm_nobetal)
  axpl = (u - velxl*velxl)*invulpfac
  axml = (u - velxl*velxl)*invulmfac
  vxpl = (velxl - lamp_nobetal)*invulpfac
  vxml = (velxl - lamm_nobetal)*invulmfac



  kappal = dpdepsl / (dpdepsl - rhol * cs2l)




  reivec1l(1) = kappal / (enthalpyl * wl)
  reivec1l(2) = vlowxl
  reivec1l(3) = vlowyl
  reivec1l(4) = vlowzl
  reivec1l(5) = one - reivec1l(1)



  reivec2l(1) = wl * vlowyl
  reivec2l(2) = enthalpyl * (gxy + two * wl * wl * vlowxl * vlowyl)
  reivec2l(3) = enthalpyl * (gyy + two * wl * wl * vlowyl * vlowyl)
  reivec2l(4) = enthalpyl * (gyz + two * wl * wl * vlowyl * vlowzl)
  reivec2l(5) = vlowyl * wl * (two * wl * enthalpyl - one)



  reivec3l(1) = wl * vlowzl
  reivec3l(2) = enthalpyl * (gxz + two * wl * wl * vlowxl * vlowzl)
  reivec3l(3) = enthalpyl * (gyz + two * wl * wl * vlowyl * vlowzl)
  reivec3l(4) = enthalpyl * (gzz + two * wl * wl * vlowzl * vlowzl)
  reivec3l(5) = vlowzl * wl * (two * wl * enthalpyl - one)



  reivecpl(1) = one
  reivecpl(2) = enthalpyl * wl * (vlowxl - vxpl)
  reivecpl(3) = enthalpyl * wl * vlowyl
  reivecpl(4) = enthalpyl * wl * vlowzl
  reivecpl(5) = enthalpyl * wl * axpl - one



  reivecml(1) = one
  reivecml(2) = enthalpyl * wl * (vlowxl - vxml)
  reivecml(3) = enthalpyl * wl * vlowyl
  reivecml(4) = enthalpyl * wl * vlowzl
  reivecml(5) = enthalpyl * wl * axml - one





  invurpfac = one / (u - velxr*lamp_nobetar)
  invurmfac = one / (u - velxr*lamm_nobetar)
  axpr = (u - velxr*velxr)*invurpfac
  axmr = (u - velxr*velxr)*invurmfac
  vxpr = (velxr - lamp_nobetar)*invurpfac
  vxmr = (velxr - lamm_nobetar)*invurmfac



  kappar = dpdepsr / (dpdepsr - rhor * cs2r)



  reivec1r(1) = kappar / (enthalpyr * wr)
  reivec1r(2) = vlowxr
  reivec1r(3) = vlowyr
  reivec1r(4) = vlowzr
  reivec1r(5) = one - reivec1r(1)



  reivec2r(1) = wr * vlowyr
  reivec2r(2) = enthalpyr * (gxy + two * wr * wr * vlowxr * vlowyr)
  reivec2r(3) = enthalpyr * (gyy + two * wr * wr * vlowyr * vlowyr)
  reivec2r(4) = enthalpyr * (gyz + two * wr * wr * vlowyr * vlowzr)
  reivec2r(5) = vlowyr * wr * (two * wr * enthalpyr - one)



  reivec3r(1) = wr * vlowzr
  reivec3r(2) = enthalpyr * (gxz + two * wr * wr * vlowxr * vlowzr)
  reivec3r(3) = enthalpyr * (gyz + two * wr * wr * vlowyr * vlowzr)
  reivec3r(4) = enthalpyr * (gzz + two * wr * wr * vlowzr * vlowzr)
  reivec3r(5) = vlowzr * wr * (two * wr * enthalpyr - one)



  reivecpr(1) = one
  reivecpr(2) = enthalpyr * wr * (vlowxr - vxpr)
  reivecpr(3) = enthalpyr * wr * vlowyr
  reivecpr(4) = enthalpyr * wr * vlowzr
  reivecpr(5) = enthalpyr * wr * axpr - one



  reivecmr(1) = one
  reivecmr(2) = enthalpyr * wr * (vlowxr - vxmr)
  reivecmr(3) = enthalpyr * wr * vlowyr
  reivecmr(4) = enthalpyr * wr * vlowzr
  reivecmr(5) = enthalpyr * wr * axmr - one



  if (ANALYTICAL) then

    cxx = gyy * gzz - gyz * gyz
    cxy = gxz * gyz - gxy * gzz
    cxz = gxy * gyz - gxz * gyy
    gam = gxx * cxx + gxy * cxy + gxz * cxz



    xsil = cxx - gam * velxl * velxl
    dltl = enthalpyl**3 * wl * (kappal - one) * (vxml - vxpl) * xsil



    tmp1 = wl / (kappal - one)

    leivec1l(1) = tmp1 * (enthalpyl - wl)
    leivec1l(2) = tmp1 * wl * velxl
    leivec1l(3) = tmp1 * wl * velyl
    leivec1l(4) = tmp1 * wl * velzl
    leivec1l(5) =-tmp1 * wl



    tmp1 = one / (xsil * enthalpyl)

    leivec2l(1) = (gyz * vlowzl - gzz * vlowyl) * tmp1
    leivec2l(2) = (gzz * vlowyl - gyz * vlowzl) * tmp1 * velxl
    leivec2l(3) = (gzz * (one - velxl * vlowxl) + gxz * vlowzl * velxl ) * tmp1
    leivec2l(4) = (gyz * (velxl * vlowxl - one) - gxz * velxl * vlowyl) * tmp1
    leivec2l(5) = (gyz * vlowzl - gzz * vlowyl) * tmp1



    leivec3l(1) = (gyz * vlowyl - gyy * vlowzl) * tmp1
    leivec3l(2) = (gyy * vlowzl - gyz * vlowyl) * tmp1 * velxl
    leivec3l(3) = (gyz * (velxl * vlowxl - one) - gxy * velxl * vlowzl) * tmp1
    leivec3l(4) = (gyy * (one - velxl * vlowxl) + gxy * velxl * vlowyl) * tmp1
    leivec3l(5) = (gyz * vlowyl - gyy * vlowzl) * tmp1



    tmp1 = enthalpyl * enthalpyl / dltl
    tmp2 = wl * wl * xsil

    leivecpl(1) = - (enthalpyl * wl * vxml * xsil + (one - kappal) * &
      (vxml * (tmp2 - cxx) - gam * velxl) - kappal * tmp2 * vxml) * tmp1
    leivecpl(2) = - (cxx * (one - kappal * axml) + (two * kappal - one) * &
      vxml * (tmp2 * velxl - cxx * velxl)) * tmp1
    leivecpl(3) = - (cxy * (one - kappal * axml) + (two * kappal - one) * &
      vxml * (tmp2 * velyl - cxy * velxl)) * tmp1
    leivecpl(4) = - (cxz * (one - kappal * axml) + (two * kappal - one) * &
      vxml * (tmp2 * velzl - cxz * velxl)) * tmp1
    leivecpl(5) = - ((one - kappal) * (vxml * (tmp2 - cxx) - gam * velxl) - &
      kappal * tmp2 * vxml) * tmp1



    leivecml(1) = (enthalpyl * wl * vxpl * xsil + (one - kappal) * &
      (vxpl * (tmp2 - cxx) - gam * velxl) - kappal * tmp2 * vxpl) * tmp1
    leivecml(2) = (cxx * (one - kappal * axpl) + (two * kappal - one) * &
      vxpl * (tmp2 * velxl - cxx * velxl)) * tmp1
    leivecml(3) = (cxy * (one - kappal * axpl) + (two * kappal - one) * &
      vxpl * (tmp2 * velyl - cxy * velxl)) * tmp1
    leivecml(4) = (cxz * (one - kappal * axpl) + (two * kappal - one) * &
      vxpl * (tmp2 * velzl - cxz * velxl)) * tmp1
    leivecml(5) = ((one - kappal) * (vxpl * (tmp2 - cxx) - gam * &
      velxl) - kappal * tmp2 * vxpl) * tmp1



    xsir = cxx - gam * velxr * velxr
    dltr = enthalpyr**3 * wr * (kappar - one) * (vxmr - vxpr) * xsir



    tmp1 = wr / (kappar - one)

    leivec1r(1) = tmp1 * (enthalpyr - wr)
    leivec1r(2) = tmp1 * wr * velxr
    leivec1r(3) = tmp1 * wr * velyr
    leivec1r(4) = tmp1 * wr * velzr
    leivec1r(5) =-tmp1 * wr



    tmp1 = one / (xsir * enthalpyr)

    leivec2r(1) = (gyz * vlowzr - gzz * vlowyr) * tmp1
    leivec2r(2) = (gzz * vlowyr - gyz * vlowzr) * tmp1 * velxr
    leivec2r(3) = (gzz * (one - velxr * vlowxr) + gxz * vlowzr * velxr) * tmp1
    leivec2r(4) = (gyz * (velxr * vlowxr - one) - gxz * velxr * vlowyr) * tmp1
    leivec2r(5) = (gyz * vlowzr - gzz * vlowyr) * tmp1



    leivec3r(1) = (gyz * vlowyr - gyy * vlowzr) * tmp1
    leivec3r(2) = (gyy * vlowzr - gyz * vlowyr) * tmp1 * velxr
    leivec3r(3) = (gyz * (velxr * vlowxr - one) - gxy * velxr * vlowzr) * tmp1
    leivec3r(4) = (gyy * (one - velxr * vlowxr) + gxy * velxr * vlowyr) * tmp1
    leivec3r(5) = (gyz * vlowyr - gyy * vlowzr) * tmp1



    tmp1 = enthalpyr * enthalpyr / dltr
    tmp2 = wr * wr * xsir

    leivecpr(1) = - (enthalpyr * wr * vxmr * xsir + (one - kappar) * &
      (vxmr * (tmp2 - cxx) - gam * velxr) - kappar * tmp2 * vxmr) * tmp1
    leivecpr(2) = - (cxx * (one - kappar * axmr) + (two * kappar - one) * &
      vxmr * (tmp2 * velxr - cxx * velxr)) * tmp1
    leivecpr(3) = - (cxy * (one - kappar * axmr) + (two * kappar - one) * &
      vxmr * (tmp2 * velyr - cxy * velxr)) * tmp1
    leivecpr(4) = - (cxz * (one - kappar * axmr) + (two * kappar - one) * &
      vxmr * (tmp2 * velzr - cxz * velxr)) * tmp1
    leivecpr(5) = - ((one - kappar) * (vxmr * (tmp2 - cxx) - gam * velxr) - &
      kappar * tmp2 * vxmr) * tmp1



    leivecmr(1) = (enthalpyr * wr * vxpr * xsir + (one - kappar) * &
      (vxpr * (tmp2 - cxx) - gam * velxr) - kappar * tmp2 * vxpr) * tmp1
    leivecmr(2) = (cxx * (one - kappar * axpr) + (two * kappar - one) * &
      vxpr * (tmp2 * velxr - cxx * velxr)) * tmp1
    leivecmr(3) = (cxy * (one - kappar * axpr) + (two * kappar - one) * &
      vxpr * (tmp2 * velyr - cxy * velxr)) * tmp1
    leivecmr(4) = (cxz * (one - kappar * axpr) + (two * kappar - one) * &
      vxpr * (tmp2 * velzr - cxz * velxr)) * tmp1
    leivecmr(5) = ((one - kappar) * (vxpr * (tmp2 - cxx) - gam * &
      velxr) - kappar * tmp2 * vxpr) * tmp1

  endif










  p(:,1) = reivecml(:)
  p(:,2) = reivecpl(:)
  p(:,3) = reivec2l(:)
  p(:,4) = reivec3l(:)
  p(:,5) = reivec1l(:)



  du(1) = densl
  du(2) = sxl
  du(3) = syl
  du(4) = szl
  du(5) = taul

  if (ANALYTICAL) then

    if (FAST) then

      sump = 0.0d0
      summ = 0.0d0

      do i=1,5
         sump = sump + (lamp - lam1) * leivecpl(i) * du(i)
         summ = summ + (lamm - lam1) * leivecml(i) * du(i)
      enddo

      vxa = sump + summ
      vxb =-(sump * vxpl + summ * vxml)

      rfluxl(1) = lam1 * du(1) + vxa
      rfluxl(2) = lam1 * du(2) + enthalpyl * wl * (vlowxl * vxa + vxb)
      rfluxl(3) = lam1 * du(3) + enthalpyl * wl * (vlowyl * vxa)
      rfluxl(4) = lam1 * du(4) + enthalpyl * wl * (vlowzl * vxa)
      rfluxl(5) = lam1 * du(5) + enthalpyl * wl * (velxl  * vxb + vxa) - vxa

    else



      q(1,:) = leivecml(:)
      q(2,:) = leivecpl(:)
      q(3,:) = leivec2l(:)
      q(4,:) = leivec3l(:)
      q(5,:) = leivec1l(:)

      do i=1,5
        dw(i) = 0.0d0
        do j=1,5
          dw(i) = dw(i) + q(i,j) * du(j)
        enddo
      enddo

      do i = 1, 5
        rfluxl(i) = 0.d0
        do j = 1, 5
          rfluxl(i) = rfluxl(i) + p(i,j) * lam(j) * dw(j)
        end do
      end do

    endif

  else



    dw=du
    aa = p

    do i=1,5
      paug(:,i) = p(:,i)
    enddo

















    paug(1,6) = du(1)
    paug(2,6) = du(2)
    paug(3,6) = du(3)
    paug(4,6) = du(4)
    paug(5,6) = du(5)


    do i=1,5

      tmp1 = one / paug(i,i)
      do j=i,6
        paug(i,j) = paug(i,j)*tmp1
      enddo

      do j=i+1,5
        tmp1 = - (paug(j,i))
        do k=i,6
          paug(j,k) = paug(j,k) + tmp1*paug(i,k)
        enddo
      enddo
    enddo

    f_du(5) = paug(5,6)
    do i=4,1,-1
      f_du(i) = paug(i,6)
      do j=i+1,5
        f_du(i) = f_du(i) - paug(i,j)*f_du(j)
      enddo
    enddo

    dw(1) = f_du(1)
    dw(2) = f_du(2)
    dw(3) = f_du(3)
    dw(4) = f_du(4)
    dw(5) = f_du(5)

    do i = 1, 5
      rfluxl(i) = 0.d0
      do j = 1, 5
        rfluxl(i) = rfluxl(i) + p(i,j) * lam(j) * dw(j)
      end do
    end do

  endif

    flux1l = rfluxl(1)
    flux2l = rfluxl(2)
    flux3l = rfluxl(3)
    flux4l = rfluxl(4)
    flux5l = rfluxl(5)











  p(:,1) = reivecmr(:)
  p(:,2) = reivecpr(:)
  p(:,3) = reivec2r(:)
  p(:,4) = reivec3r(:)
  p(:,5) = reivec1r(:)

  du(1) = densr
  du(2) = sxr
  du(3) = syr
  du(4) = szr
  du(5) = taur

  if (ANALYTICAL) then

    if (FAST) then

      sump = 0.0d0
      summ = 0.0d0

      do i=1,5
         sump = sump + (lamp - lam1) * leivecpr(i) * du(i)
         summ = summ + (lamm - lam1) * leivecmr(i) * du(i)
      enddo

      vxa = sump + summ
      vxb =-(sump * vxpr + summ * vxmr)

      rfluxr(1) = lam1 * du(1) + vxa
      rfluxr(2) = lam1 * du(2) + enthalpyr * wr * (vlowxr * vxa + vxb)
      rfluxr(3) = lam1 * du(3) + enthalpyr * wr * (vlowyr * vxa)
      rfluxr(4) = lam1 * du(4) + enthalpyr * wr * (vlowzr * vxa)
      rfluxr(5) = lam1 * du(5) + enthalpyr * wr * (velxr  * vxb + vxa) - vxa

    else



      q(1,:) = leivecmr(:)
      q(2,:) = leivecpr(:)
      q(3,:) = leivec2r(:)
      q(4,:) = leivec3r(:)
      q(5,:) = leivec1r(:)

      do i=1,5
        dw(i) = 0.0d0
        do j=1,5
          dw(i) = dw(i) + q(i,j) * du(j)
        enddo
      enddo

      do i = 1, 5
        rfluxr(i) = 0.d0
        do j = 1, 5
          rfluxr(i) = rfluxr(i) + p(i,j) * lam(j) * dw(j)
        end do
      end do

    endif

  else



    do i=1,5
      dw(i) = du(i)
      do j=1,5
        aa(i,j) = p(i,j)
      end do
    enddo

    do i=1,5
      paug(i,1) = p(i,1)
      paug(i,2) = p(i,2)
      paug(i,3) = p(i,3)
      paug(i,4) = p(i,4)
      paug(i,5) = p(i,5)
    enddo

    paug(1,6) = du(1)
    paug(2,6) = du(2)
    paug(3,6) = du(3)
    paug(4,6) = du(4)
    paug(5,6) = du(5)



    do i=1,5

      tmp1 = one / paug(i,i)
      do j=i,6
        paug(i,j) = paug(i,j)*tmp1
      enddo

      do j=i+1,5
        tmp1 = - (paug(j,i))
        do k=i,6
          paug(j,k) = paug(j,k) + tmp1*paug(i,k)
        enddo
      enddo
    enddo


    f_du(5) = paug(5,6)
    do i=4,1,-1
      f_du(i) = paug(i,6)
      do j=i+1,5
        f_du(i) = f_du(i) - paug(i,j)*f_du(j)
      enddo
    enddo

    dw(1) = f_du(1)
    dw(2) = f_du(2)
    dw(3) = f_du(3)
    dw(4) = f_du(4)
    dw(5) = f_du(5)

    do i = 1, 5
      rfluxr(i) = 0.d0
      do j = 1, 5
        rfluxr(i) = rfluxr(i) + p(i,j) * lam(j) * dw(j)
      end do
    end do

  endif

  flux1r = rfluxr(1)
  flux2r = rfluxr(2)
  flux3r = rfluxr(3)
  flux4r = rfluxr(4)
  flux5r = rfluxr(5)

  flux1 = flux1r - flux1l 
  flux2 = flux2r - flux2l 
  flux3 = flux3r - flux3l 
  flux4 = flux4r - flux4l 
  flux5 = flux5r - flux5l 

  return

end subroutine eigenproblem_marquina

subroutine eigenproblem_marquina_hot(handle,keytemp,rhor,velxr,velyr,&
     velzr,epsr,rhol,velxl,velyl,velzl,epsl,templ,tempr,ye_l,ye_r,gxx,gxy,gxz,&
     gyy,gyz,gzz,u,det,alp,beta,densl,sxl,syl,szl,taul,&
     densr,sxr,syr,szr,taur,flux1,flux2,flux3,flux4,flux5)

  USE GRHydro_Scalars

  implicit none

  REAL*8 GRHydro_Y_e_max
  integer, parameter :: cctki_use_GRHydro_Y_e_max = kind(GRHydro_Y_e_max)
  REAL*8 GRHydro_Y_e_min
  integer, parameter :: cctki_use_GRHydro_Y_e_min = kind(GRHydro_Y_e_min)
  REAL*8 GRHydro_atmo_tolerance
  integer, parameter :: cctki_use_GRHydro_atmo_tolerance = kind(GRHydro_atmo_tolerance)
  REAL*8 GRHydro_del_ptol
  integer, parameter :: cctki_use_GRHydro_del_ptol = kind(GRHydro_del_ptol)
  REAL*8 GRHydro_eos_rf_prec
  integer, parameter :: cctki_use_GRHydro_eos_rf_prec = kind(GRHydro_eos_rf_prec)
  REAL*8 GRHydro_eps_min
  integer, parameter :: cctki_use_GRHydro_eps_min = kind(GRHydro_eps_min)
  REAL*8 GRHydro_hot_atmo_Y_e
  integer, parameter :: cctki_use_GRHydro_hot_atmo_Y_e = kind(GRHydro_hot_atmo_Y_e)
  REAL*8 GRHydro_hot_atmo_temp
  integer, parameter :: cctki_use_GRHydro_hot_atmo_temp = kind(GRHydro_hot_atmo_temp)
  REAL*8 GRHydro_lorentz_overshoot_cutoff
  integer, parameter :: cctki_use_GRHydro_lorentz_overshoot_cutoff = kind(GRHydro_lorentz_overshoot_cutoff)
  REAL*8 GRHydro_max_temp
  integer, parameter :: cctki_use_GRHydro_max_temp = kind(GRHydro_max_temp)
  REAL*8 GRHydro_perc_ptol
  integer, parameter :: cctki_use_GRHydro_perc_ptol = kind(GRHydro_perc_ptol)
  REAL*8 GRHydro_rho_central
  integer, parameter :: cctki_use_GRHydro_rho_central = kind(GRHydro_rho_central)
  REAL*8 Tmunu_damping_radius_max
  integer, parameter :: cctki_use_Tmunu_damping_radius_max = kind(Tmunu_damping_radius_max)
  REAL*8 Tmunu_damping_radius_min
  integer, parameter :: cctki_use_Tmunu_damping_radius_min = kind(Tmunu_damping_radius_min)
  REAL*8 atmo_falloff_power
  integer, parameter :: cctki_use_atmo_falloff_power = kind(atmo_falloff_power)
  REAL*8 atmo_falloff_radius
  integer, parameter :: cctki_use_atmo_falloff_radius = kind(atmo_falloff_radius)
  REAL*8 atmo_tolerance_power
  integer, parameter :: cctki_use_atmo_tolerance_power = kind(atmo_tolerance_power)
  REAL*8 atmo_tolerance_radius
  integer, parameter :: cctki_use_atmo_tolerance_radius = kind(atmo_tolerance_radius)
  REAL*8 c2p_reset_pressure_to_value
  integer, parameter :: cctki_use_c2p_reset_pressure_to_value = kind(c2p_reset_pressure_to_value)
  REAL*8 comoving_factor
  integer, parameter :: cctki_use_comoving_factor = kind(comoving_factor)
  REAL*8 comoving_tanh_factor
  integer, parameter :: cctki_use_comoving_tanh_factor = kind(comoving_tanh_factor)
  REAL*8 comoving_tanh_offset
  integer, parameter :: cctki_use_comoving_tanh_offset = kind(comoving_tanh_offset)
  REAL*8 enhanced_ppm_C2
  integer, parameter :: cctki_use_enhanced_ppm_C2 = kind(enhanced_ppm_C2)
  REAL*8 enhanced_ppm_C3
  integer, parameter :: cctki_use_enhanced_ppm_C3 = kind(enhanced_ppm_C3)
  REAL*8 initial_Gamma
  integer, parameter :: cctki_use_initial_Gamma = kind(initial_Gamma)
  REAL*8 initial_atmosphere_factor
  integer, parameter :: cctki_use_initial_atmosphere_factor = kind(initial_atmosphere_factor)
  REAL*8 initial_k
  integer, parameter :: cctki_use_initial_k = kind(initial_k)
  REAL*8 initial_rho_abs_min
  integer, parameter :: cctki_use_initial_rho_abs_min = kind(initial_rho_abs_min)
  REAL*8 initial_rho_rel_min
  integer, parameter :: cctki_use_initial_rho_rel_min = kind(initial_rho_rel_min)
  REAL*8 kap_dc
  integer, parameter :: cctki_use_kap_dc = kind(kap_dc)
  REAL*8 max_magnetic_to_gas_pressure_ratio
  integer, parameter :: cctki_use_max_magnetic_to_gas_pressure_ratio = kind(max_magnetic_to_gas_pressure_ratio)
  REAL*8 min_tracer
  integer, parameter :: cctki_use_min_tracer = kind(min_tracer)
  REAL*8 mp5_alpha
  integer, parameter :: cctki_use_mp5_alpha = kind(mp5_alpha)
  REAL*8 mp5_eps
  integer, parameter :: cctki_use_mp5_eps = kind(mp5_eps)
  REAL*8 myfloor
  integer, parameter :: cctki_use_myfloor = kind(myfloor)
  REAL*8 ppm_epsilon
  integer, parameter :: cctki_use_ppm_epsilon = kind(ppm_epsilon)
  REAL*8 ppm_epsilon_shock
  integer, parameter :: cctki_use_ppm_epsilon_shock = kind(ppm_epsilon_shock)
  REAL*8 ppm_eta1
  integer, parameter :: cctki_use_ppm_eta1 = kind(ppm_eta1)
  REAL*8 ppm_eta2
  integer, parameter :: cctki_use_ppm_eta2 = kind(ppm_eta2)
  REAL*8 ppm_k0
  integer, parameter :: cctki_use_ppm_k0 = kind(ppm_k0)
  REAL*8 ppm_omega1
  integer, parameter :: cctki_use_ppm_omega1 = kind(ppm_omega1)
  REAL*8 ppm_omega2
  integer, parameter :: cctki_use_ppm_omega2 = kind(ppm_omega2)
  REAL*8 ppm_omega_tracer
  integer, parameter :: cctki_use_ppm_omega_tracer = kind(ppm_omega_tracer)
  REAL*8 ppm_small
  integer, parameter :: cctki_use_ppm_small = kind(ppm_small)
  REAL*8 rho_abs_min
  integer, parameter :: cctki_use_rho_abs_min = kind(rho_abs_min)
  REAL*8 rho_abs_min_after_recovery
  integer, parameter :: cctki_use_rho_abs_min_after_recovery = kind(rho_abs_min_after_recovery)
  REAL*8 rho_rel_min
  integer, parameter :: cctki_use_rho_rel_min = kind(rho_rel_min)
  REAL*8 setcharmax
  integer, parameter :: cctki_use_setcharmax = kind(setcharmax)
  REAL*8 setcharmin
  integer, parameter :: cctki_use_setcharmin = kind(setcharmin)
  REAL*8 sqrtdet_thr
  integer, parameter :: cctki_use_sqrtdet_thr = kind(sqrtdet_thr)
  REAL*8 tau_rel_min
  integer, parameter :: cctki_use_tau_rel_min = kind(tau_rel_min)
  REAL*8 weno_eps
  integer, parameter :: cctki_use_weno_eps = kind(weno_eps)
  integer*8 Avec_gauge
  integer, parameter :: cctki_use_Avec_gauge = kind(Avec_gauge)
  integer*8 EoS_Change_type
  integer, parameter :: cctki_use_EoS_Change_type = kind(EoS_Change_type)
  integer*8 GRHydro_c2p_failed_action
  integer, parameter :: cctki_use_GRHydro_c2p_failed_action = kind(GRHydro_c2p_failed_action)
  integer*8 GRHydro_eos_table
  integer, parameter :: cctki_use_GRHydro_eos_table = kind(GRHydro_eos_table)
  integer*8 GRHydro_eos_type
  integer, parameter :: cctki_use_GRHydro_eos_type = kind(GRHydro_eos_type)
  integer*8 HLLE_type
  integer, parameter :: cctki_use_HLLE_type = kind(HLLE_type)
  integer*8 bound
  integer, parameter :: cctki_use_bound = kind(bound)
  integer*8 comoving_attenuate
  integer, parameter :: cctki_use_comoving_attenuate = kind(comoving_attenuate)
  integer*8 comoving_v_method
  integer, parameter :: cctki_use_comoving_v_method = kind(comoving_v_method)
  integer*8 gradient_method
  integer, parameter :: cctki_use_gradient_method = kind(gradient_method)
  integer*8 left_eigenvectors
  integer, parameter :: cctki_use_left_eigenvectors = kind(left_eigenvectors)
  integer*8 method_type
  integer, parameter :: cctki_use_method_type = kind(method_type)
  integer*8 numerical_viscosity
  integer, parameter :: cctki_use_numerical_viscosity = kind(numerical_viscosity)
  integer*8 particle_interpolator
  integer, parameter :: cctki_use_particle_interpolator = kind(particle_interpolator)
  integer*8 ppm_flatten
  integer, parameter :: cctki_use_ppm_flatten = kind(ppm_flatten)
  integer*8 psidcspeed
  integer, parameter :: cctki_use_psidcspeed = kind(psidcspeed)
  integer*8 recon_method
  integer, parameter :: cctki_use_recon_method = kind(recon_method)
  integer*8 recon_vars
  integer, parameter :: cctki_use_recon_vars = kind(recon_vars)
  integer*8 riemann_solver
  integer, parameter :: cctki_use_riemann_solver = kind(riemann_solver)
  integer*8 tvd_limiter
  integer, parameter :: cctki_use_tvd_limiter = kind(tvd_limiter)
  integer*8 use_evolution_mask
  integer, parameter :: cctki_use_use_evolution_mask = kind(use_evolution_mask)
  INTEGER*4 Check_Rho_Minimum
  integer, parameter :: cctki_use_Check_Rho_Minimum = kind(Check_Rho_Minimum)
  INTEGER*4 EoS_Change
  integer, parameter :: cctki_use_EoS_Change = kind(EoS_Change)
  INTEGER*4 GRHydro_MaxNumConstrainedVars
  integer, parameter :: cctki_use_GRHydro_MaxNumConstrainedVars = kind(GRHydro_MaxNumConstrainedVars)
  INTEGER*4 GRHydro_MaxNumEvolvedVars
  integer, parameter :: cctki_use_GRHydro_MaxNumEvolvedVars = kind(GRHydro_MaxNumEvolvedVars)
  INTEGER*4 GRHydro_MaxNumEvolvedVarsSlow
  integer, parameter :: cctki_use_GRHydro_MaxNumEvolvedVarsSlow = kind(GRHydro_MaxNumEvolvedVarsSlow)
  INTEGER*4 GRHydro_MaxNumSandRVars
  integer, parameter :: cctki_use_GRHydro_MaxNumSandRVars = kind(GRHydro_MaxNumSandRVars)
  INTEGER*4 GRHydro_NaN_verbose
  integer, parameter :: cctki_use_GRHydro_NaN_verbose = kind(GRHydro_NaN_verbose)
  INTEGER*4 GRHydro_c2p_reset_eps_tau_hot_eos
  integer, parameter :: cctki_use_GRHydro_c2p_reset_eps_tau_hot_eos = kind(GRHydro_c2p_reset_eps_tau_hot_eos)
  INTEGER*4 GRHydro_c2p_warn_from_reflevel
  integer, parameter :: cctki_use_GRHydro_c2p_warn_from_reflevel = kind(GRHydro_c2p_warn_from_reflevel)
  INTEGER*4 GRHydro_c2p_warnlevel
  integer, parameter :: cctki_use_GRHydro_c2p_warnlevel = kind(GRHydro_c2p_warnlevel)
  INTEGER*4 GRHydro_countmax
  integer, parameter :: cctki_use_GRHydro_countmax = kind(GRHydro_countmax)
  INTEGER*4 GRHydro_countmin
  integer, parameter :: cctki_use_GRHydro_countmin = kind(GRHydro_countmin)
  INTEGER*4 GRHydro_enable_internal_excision
  integer, parameter :: cctki_use_GRHydro_enable_internal_excision = kind(GRHydro_enable_internal_excision)
  INTEGER*4 GRHydro_eos_hot_eps_fix
  integer, parameter :: cctki_use_GRHydro_eos_hot_eps_fix = kind(GRHydro_eos_hot_eps_fix)
  INTEGER*4 GRHydro_eos_hot_prim2con_warn
  integer, parameter :: cctki_use_GRHydro_eos_hot_prim2con_warn = kind(GRHydro_eos_hot_prim2con_warn)
  INTEGER*4 GRHydro_hydro_excision
  integer, parameter :: cctki_use_GRHydro_hydro_excision = kind(GRHydro_hydro_excision)
  INTEGER*4 GRHydro_oppm_reflevel
  integer, parameter :: cctki_use_GRHydro_oppm_reflevel = kind(GRHydro_oppm_reflevel)
  INTEGER*4 GRHydro_polish
  integer, parameter :: cctki_use_GRHydro_polish = kind(GRHydro_polish)
  INTEGER*4 GRHydro_stencil
  integer, parameter :: cctki_use_GRHydro_stencil = kind(GRHydro_stencil)
  INTEGER*4 WENO_order
  integer, parameter :: cctki_use_WENO_order = kind(WENO_order)
  INTEGER*4 apply_H_viscosity
  integer, parameter :: cctki_use_apply_H_viscosity = kind(apply_H_viscosity)
  INTEGER*4 c2p_reset_pressure
  integer, parameter :: cctki_use_c2p_reset_pressure = kind(c2p_reset_pressure)
  INTEGER*4 c2p_resort_to_bisection
  integer, parameter :: cctki_use_c2p_resort_to_bisection = kind(c2p_resort_to_bisection)
  INTEGER*4 calculate_bcom
  integer, parameter :: cctki_use_calculate_bcom = kind(calculate_bcom)
  INTEGER*4 check_for_trivial_rp
  integer, parameter :: cctki_use_check_for_trivial_rp = kind(check_for_trivial_rp)
  INTEGER*4 clean_divergence
  integer, parameter :: cctki_use_clean_divergence = kind(clean_divergence)
  INTEGER*4 con2prim_oct_hack
  integer, parameter :: cctki_use_con2prim_oct_hack = kind(con2prim_oct_hack)
  INTEGER*4 decouple_normal_Bfield
  integer, parameter :: cctki_use_decouple_normal_Bfield = kind(decouple_normal_Bfield)
  INTEGER*4 eno_order
  integer, parameter :: cctki_use_eno_order = kind(eno_order)
  INTEGER*4 evolve_tracer
  integer, parameter :: cctki_use_evolve_tracer = kind(evolve_tracer)
  INTEGER*4 mp5_adaptive_eps
  integer, parameter :: cctki_use_mp5_adaptive_eps = kind(mp5_adaptive_eps)
  INTEGER*4 number_of_arrays
  integer, parameter :: cctki_use_number_of_arrays = kind(number_of_arrays)
  INTEGER*4 number_of_particles
  integer, parameter :: cctki_use_number_of_particles = kind(number_of_particles)
  INTEGER*4 number_of_tracers
  integer, parameter :: cctki_use_number_of_tracers = kind(number_of_tracers)
  INTEGER*4 particle_interpolation_order
  integer, parameter :: cctki_use_particle_interpolation_order = kind(particle_interpolation_order)
  INTEGER*4 ppm_detect
  integer, parameter :: cctki_use_ppm_detect = kind(ppm_detect)
  INTEGER*4 ppm_mppm
  integer, parameter :: cctki_use_ppm_mppm = kind(ppm_mppm)
  INTEGER*4 ppm_mppm_debug_eigenvalues
  integer, parameter :: cctki_use_ppm_mppm_debug_eigenvalues = kind(ppm_mppm_debug_eigenvalues)
  INTEGER*4 reconstruct_Wv
  integer, parameter :: cctki_use_reconstruct_Wv = kind(reconstruct_Wv)
  INTEGER*4 reconstruct_temper
  integer, parameter :: cctki_use_reconstruct_temper = kind(reconstruct_temper)
  INTEGER*4 set_trivial_rp_grid_function
  integer, parameter :: cctki_use_set_trivial_rp_grid_function = kind(set_trivial_rp_grid_function)
  INTEGER*4 sources_spatial_order
  integer, parameter :: cctki_use_sources_spatial_order = kind(sources_spatial_order)
  INTEGER*4 sync_conserved_only
  integer, parameter :: cctki_use_sync_conserved_only = kind(sync_conserved_only)
  INTEGER*4 track_divB
  integer, parameter :: cctki_use_track_divB = kind(track_divB)
  INTEGER*4 transport_constraints
  integer, parameter :: cctki_use_transport_constraints = kind(transport_constraints)
  INTEGER*4 use_MoL_slow_multirate_sector
  integer, parameter :: cctki_use_use_MoL_slow_multirate_sector = kind(use_MoL_slow_multirate_sector)
  INTEGER*4 use_enhanced_ppm
  integer, parameter :: cctki_use_use_enhanced_ppm = kind(use_enhanced_ppm)
  INTEGER*4 use_min_tracer
  integer, parameter :: cctki_use_use_min_tracer = kind(use_min_tracer)
  INTEGER*4 use_optimized_ppm
  integer, parameter :: cctki_use_use_optimized_ppm = kind(use_optimized_ppm)
  INTEGER*4 use_weighted_fluxes
  integer, parameter :: cctki_use_use_weighted_fluxes = kind(use_weighted_fluxes)
  INTEGER*4 weno_adaptive_epsilon
  integer, parameter :: cctki_use_weno_adaptive_epsilon = kind(weno_adaptive_epsilon)
  INTEGER*4 wk_atmosphere
  integer, parameter :: cctki_use_wk_atmosphere = kind(wk_atmosphere)
  COMMON /GRHydrorest/GRHydro_Y_e_max, GRHydro_Y_e_min, GRHydro_atmo_tolerance, GRHydro_del_ptol, GRHydro_eos_rf_prec, GRHydro_eps_&
  &min, GRHydro_hot_atmo_Y_e, GRHydro_hot_atmo_temp, GRHydro_lorentz_overshoot_cutoff, GRHydro_max_temp, GRHydro_perc_ptol, GRHydro&
  &_rho_central, Tmunu_damping_radius_max, Tmunu_damping_radius_min, atmo_falloff_power, atmo_falloff_radius, atmo_tolerance_power,&
  & atmo_tolerance_radius, c2p_reset_pressure_to_value, comoving_factor, comoving_tanh_factor, comoving_tanh_offset, enhanced_ppm_C&
  &2, enhanced_ppm_C3, initial_Gamma, initial_atmosphere_factor, initial_k, initial_rho_abs_min, initial_rho_rel_min, kap_dc, max_m&
  &agnetic_to_gas_pressure_ratio, min_tracer, mp5_alpha, mp5_eps, myfloor, ppm_epsilon, ppm_epsilon_shock, ppm_eta1, ppm_eta2, ppm_&
  &k0, ppm_omega1, ppm_omega2, ppm_omega_tracer, ppm_small, rho_abs_min, rho_abs_min_after_recovery, rho_rel_min, setcharmax, setch&
  &armin, sqrtdet_thr, tau_rel_min, weno_eps, Avec_gauge, EoS_Change_type, GRHydro_c2p_failed_action, GRHydro_eos_table, GRHydro_eo&
  &s_type, HLLE_type, bound, comoving_attenuate, comoving_v_method, gradient_method, left_eigenvectors, method_type, numerical_visc&
  &osity, particle_interpolator, ppm_flatten, psidcspeed, recon_method, recon_vars, riemann_solver, tvd_limiter, use_evolution_mask&
  &, Check_Rho_Minimum, EoS_Change, GRHydro_MaxNumConstrainedVars, GRHydro_MaxNumEvolvedVars, GRHydro_MaxNumEvolvedVarsSlow, GRHydr&
  &o_MaxNumSandRVars, GRHydro_NaN_verbose, GRHydro_c2p_reset_eps_tau_hot_eos, GRHydro_c2p_warn_from_reflevel, GRHydro_c2p_warnlevel&
  &, GRHydro_countmax, GRHydro_countmin, GRHydro_enable_internal_excision, GRHydro_eos_hot_eps_fix, GRHydro_eos_hot_prim2con_warn, &
  &GRHydro_hydro_excision, GRHydro_oppm_reflevel, GRHydro_polish, GRHydro_stencil, WENO_order, apply_H_viscosity, c2p_reset_pressur&
  &e, c2p_resort_to_bisection, calculate_bcom, check_for_trivial_rp, clean_divergence, con2prim_oct_hack, decouple_normal_Bfield, e&
  &no_order, evolve_tracer, mp5_adaptive_eps, number_of_arrays, number_of_particles, number_of_tracers, particle_interpolation_orde&
  &r, ppm_detect, ppm_mppm, ppm_mppm_debug_eigenvalues, reconstruct_Wv, reconstruct_temper, set_trivial_rp_grid_function, sources_s&
  &patial_order, sync_conserved_only, track_divB, transport_constraints, use_MoL_slow_multirate_sector, use_enhanced_ppm, use_min_t&
  &racer, use_optimized_ppm, use_weighted_fluxes, weno_adaptive_epsilon, wk_atmosphere
  INTEGER*4 constrain_to_1D
  integer, parameter :: cctki_use_constrain_to_1D = kind(constrain_to_1D)
  INTEGER*4 use_cxx_code
  integer, parameter :: cctki_use_use_cxx_code = kind(use_cxx_code)
  INTEGER*4 verbose
  integer, parameter :: cctki_use_verbose = kind(verbose)
  COMMON /GRHydropriv/constrain_to_1D, use_cxx_code, verbose
  integer*8 CCTKH0
  integer, parameter :: cctki_use_CCTKH0 = kind(CCTKH0)
  integer*8 CCTKH1
  integer, parameter :: cctki_use_CCTKH1 = kind(CCTKH1)
  integer*8 CCTKH2
  integer, parameter :: cctki_use_CCTKH2 = kind(CCTKH2)
  integer*8 CCTKH3
  integer, parameter :: cctki_use_CCTKH3 = kind(CCTKH3)
  integer*8 CCTKH4
  integer, parameter :: cctki_use_CCTKH4 = kind(CCTKH4)
  integer*8 CCTKH5
  integer, parameter :: cctki_use_CCTKH5 = kind(CCTKH5)
  integer*8 CCTKH6
  integer, parameter :: cctki_use_CCTKH6 = kind(CCTKH6)
  integer*8 CCTKH7
  integer, parameter :: cctki_use_CCTKH7 = kind(CCTKH7)
  integer*8 initial_shift
  integer, parameter :: cctki_use_initial_shift = kind(initial_shift)
  integer*8 CCTKH8
  integer, parameter :: cctki_use_CCTKH8 = kind(CCTKH8)
  integer*8 CCTKH9
  integer, parameter :: cctki_use_CCTKH9 = kind(CCTKH9)
  integer*8 CCTKH11
  integer, parameter :: cctki_use_CCTKH11 = kind(CCTKH11)
  integer*8 metric_type
  integer, parameter :: cctki_use_metric_type = kind(metric_type)
  integer*8 shift_evolution_method
  integer, parameter :: cctki_use_shift_evolution_method = kind(shift_evolution_method)
  integer*8 CCTKH13
  integer, parameter :: cctki_use_CCTKH13 = kind(CCTKH13)
  INTEGER*4 CCTKH10
  integer, parameter :: cctki_use_CCTKH10 = kind(CCTKH10)
  INTEGER*4 CCTKH12
  integer, parameter :: cctki_use_CCTKH12 = kind(CCTKH12)
  INTEGER*4 CCTKH14
  integer, parameter :: cctki_use_CCTKH14 = kind(CCTKH14)
  COMMON /ADMBASErest/CCTKH0, CCTKH1, CCTKH2, CCTKH3, CCTKH4, CCTKH5, CCTKH6, CCTKH7, initial_shift, CCTKH8, CCTKH9, CCTKH11, metri&
  &c_type, shift_evolution_method, CCTKH13, CCTKH10, CCTKH12, CCTKH14
  REAL*8 CCTKH15
  integer, parameter :: cctki_use_CCTKH15 = kind(CCTKH15)
  REAL*8 CCTKH25
  integer, parameter :: cctki_use_CCTKH25 = kind(CCTKH25)
  REAL*8 CCTKH26
  integer, parameter :: cctki_use_CCTKH26 = kind(CCTKH26)
  REAL*8 CCTKH27 (10)
  integer, parameter :: cctki_use_CCTKH27 = kind(CCTKH27)
  REAL*8 CCTKH28
  integer, parameter :: cctki_use_CCTKH28 = kind(CCTKH28)
  REAL*8 CCTKH29
  integer, parameter :: cctki_use_CCTKH29 = kind(CCTKH29)
  REAL*8 CCTKH30
  integer, parameter :: cctki_use_CCTKH30 = kind(CCTKH30)
  REAL*8 CCTKH31
  integer, parameter :: cctki_use_CCTKH31 = kind(CCTKH31)
  REAL*8 CCTKH32
  integer, parameter :: cctki_use_CCTKH32 = kind(CCTKH32)
  REAL*8 CCTKH33 (9)
  integer, parameter :: cctki_use_CCTKH33 = kind(CCTKH33)
  REAL*8 CCTKH34
  integer, parameter :: cctki_use_CCTKH34 = kind(CCTKH34)
  REAL*8 CCTKH38
  integer, parameter :: cctki_use_CCTKH38 = kind(CCTKH38)
  REAL*8 CCTKH39
  integer, parameter :: cctki_use_CCTKH39 = kind(CCTKH39)
  REAL*8 poly_k
  integer, parameter :: cctki_use_poly_k = kind(poly_k)
  integer*8 CCTKH17
  integer, parameter :: cctki_use_CCTKH17 = kind(CCTKH17)
  integer*8 CCTKH20
  integer, parameter :: cctki_use_CCTKH20 = kind(CCTKH20)
  integer*8 CCTKH24
  integer, parameter :: cctki_use_CCTKH24 = kind(CCTKH24)
  integer*8 CCTKH37
  integer, parameter :: cctki_use_CCTKH37 = kind(CCTKH37)
  INTEGER*4 CCTKH16
  integer, parameter :: cctki_use_CCTKH16 = kind(CCTKH16)
  INTEGER*4 CCTKH18
  integer, parameter :: cctki_use_CCTKH18 = kind(CCTKH18)
  INTEGER*4 CCTKH19
  integer, parameter :: cctki_use_CCTKH19 = kind(CCTKH19)
  INTEGER*4 CCTKH21
  integer, parameter :: cctki_use_CCTKH21 = kind(CCTKH21)
  INTEGER*4 CCTKH22
  integer, parameter :: cctki_use_CCTKH22 = kind(CCTKH22)
  INTEGER*4 CCTKH23
  integer, parameter :: cctki_use_CCTKH23 = kind(CCTKH23)
  INTEGER*4 CCTKH35
  integer, parameter :: cctki_use_CCTKH35 = kind(CCTKH35)
  INTEGER*4 CCTKH36
  integer, parameter :: cctki_use_CCTKH36 = kind(CCTKH36)
  INTEGER*4 CCTKH40
  integer, parameter :: cctki_use_CCTKH40 = kind(CCTKH40)
  INTEGER*4 CCTKH41
  integer, parameter :: cctki_use_CCTKH41 = kind(CCTKH41)
  COMMON /EOS_OMNIrest/CCTKH15, CCTKH25, CCTKH26, CCTKH27, CCTKH28, CCTKH29, CCTKH30, CCTKH31, CCTKH32, CCTKH33, CCTKH34, CCTKH38, &
  &CCTKH39, poly_k, CCTKH17, CCTKH20, CCTKH24, CCTKH37, CCTKH16, CCTKH18, CCTKH19, CCTKH21, CCTKH22, CCTKH23, CCTKH35, CCTKH36, CCT&
  &KH40, CCTKH41
  integer*8 CCTKH42
  integer, parameter :: cctki_use_CCTKH42 = kind(CCTKH42)
  integer*8 Bvec_evolution_method
  integer, parameter :: cctki_use_Bvec_evolution_method = kind(Bvec_evolution_method)
  integer*8 Y_e_evolution_method
  integer, parameter :: cctki_use_Y_e_evolution_method = kind(Y_e_evolution_method)
  integer*8 entropy_evolution_method
  integer, parameter :: cctki_use_entropy_evolution_method = kind(entropy_evolution_method)
  integer*8 evolution_method
  integer, parameter :: cctki_use_evolution_method = kind(evolution_method)
  integer*8 CCTKH43
  integer, parameter :: cctki_use_CCTKH43 = kind(CCTKH43)
  integer*8 CCTKH44
  integer, parameter :: cctki_use_CCTKH44 = kind(CCTKH44)
  integer*8 CCTKH45
  integer, parameter :: cctki_use_CCTKH45 = kind(CCTKH45)
  integer*8 initial_Bvec
  integer, parameter :: cctki_use_initial_Bvec = kind(initial_Bvec)
  integer*8 CCTKH46
  integer, parameter :: cctki_use_CCTKH46 = kind(CCTKH46)
  integer*8 CCTKH47
  integer, parameter :: cctki_use_CCTKH47 = kind(CCTKH47)
  integer*8 CCTKH48
  integer, parameter :: cctki_use_CCTKH48 = kind(CCTKH48)
  integer*8 CCTKH49
  integer, parameter :: cctki_use_CCTKH49 = kind(CCTKH49)
  integer*8 prolongation_type
  integer, parameter :: cctki_use_prolongation_type = kind(prolongation_type)
  integer*8 temperature_evolution_method
  integer, parameter :: cctki_use_temperature_evolution_method = kind(temperature_evolution_method)
  INTEGER*4 hydro_excision
  integer, parameter :: cctki_use_hydro_excision = kind(hydro_excision)
  INTEGER*4 timelevels
  integer, parameter :: cctki_use_timelevels = kind(timelevels)
  COMMON /HYDROBASErest/CCTKH42, Bvec_evolution_method, Y_e_evolution_method, entropy_evolution_method, evolution_method, CCTKH43, &
  &CCTKH44, CCTKH45, initial_Bvec, CCTKH46, CCTKH47, CCTKH48, CCTKH49, prolongation_type, temperature_evolution_method, hydro_excis&
  &ion, timelevels
  INTEGER*4 MoL_Max_Evolved_Array_Size
  integer, parameter :: cctki_use_MoL_Max_Evolved_Array_Size = kind(MoL_Max_Evolved_Array_Size)
  INTEGER*4 CCTKH50
  integer, parameter :: cctki_use_CCTKH50 = kind(CCTKH50)
  INTEGER*4 MoL_Num_ArrayEvolved_Vars
  integer, parameter :: cctki_use_MoL_Num_ArrayEvolved_Vars = kind(MoL_Num_ArrayEvolved_Vars)
  INTEGER*4 CCTKH51
  integer, parameter :: cctki_use_CCTKH51 = kind(CCTKH51)
  INTEGER*4 MoL_Num_Constrained_Vars
  integer, parameter :: cctki_use_MoL_Num_Constrained_Vars = kind(MoL_Num_Constrained_Vars)
  INTEGER*4 MoL_Num_Evolved_Vars
  integer, parameter :: cctki_use_MoL_Num_Evolved_Vars = kind(MoL_Num_Evolved_Vars)
  INTEGER*4 MoL_Num_Evolved_Vars_Slow
  integer, parameter :: cctki_use_MoL_Num_Evolved_Vars_Slow = kind(MoL_Num_Evolved_Vars_Slow)
  INTEGER*4 MoL_Num_SaveAndRestore_Vars
  integer, parameter :: cctki_use_MoL_Num_SaveAndRestore_Vars = kind(MoL_Num_SaveAndRestore_Vars)
  INTEGER*4 CCTKH52
  integer, parameter :: cctki_use_CCTKH52 = kind(CCTKH52)
  COMMON /METHODOFLINESrest/MoL_Max_Evolved_Array_Size, CCTKH50, MoL_Num_ArrayEvolved_Vars, CCTKH51, MoL_Num_Constrained_Vars, MoL_&
  &Num_Evolved_Vars, MoL_Num_Evolved_Vars_Slow, MoL_Num_SaveAndRestore_Vars, CCTKH52
  INTEGER*4 use_mask
  integer, parameter :: cctki_use_use_mask = kind(use_mask)
  COMMON /SPACEMASKrest/use_mask
  

  REAL*8 rhor,velxr,velyr,velzr,epsr
  REAL*8 rhol,velxl,velyl,velzl,epsl
  REAL*8 densl,sxl,syl,szl,taul
  REAL*8 densr,sxr,syr,szr,taur
  REAL*8 templ,tempr,ye_l,ye_r

  REAL*8 lam(5),p(5,5),q(5,5),dw(5)
  REAL*8 rfluxr(5),rfluxl(5)
  REAL*8 gxx,gxy,gxz,gyy,gyz,gzz,u,det
  REAL*8 alp,beta,flux1r,flux2r,flux3r,flux4r,flux5r
  REAL*8 flux1l,flux2l,flux3l,flux4l,flux5l
  REAL*8 flux1,flux2,flux3,flux4,flux5



  REAL*8 du(5),aa(5,5)
  integer i,j,k
  REAL*8 paug(5,6),tmp1,tmp2,sump,summ,f_du(5)
  REAL*8 leivec1l(5),leivec2l(5),leivec3l(5),leivecpl(5),leivecml(5)
  REAL*8 leivec1r(5),leivec2r(5),leivec3r(5),leivecpr(5),leivecmr(5)
  REAL*8 reivec1l(5),reivec2l(5),reivec3l(5),reivecpl(5),reivecml(5)
  REAL*8 reivec1r(5),reivec2r(5),reivec3r(5),reivecpr(5),reivecmr(5)
  REAL*8 lam1l,lam2l,lam3l,lamml,lampl,lamm_nobetal,lamp_nobetal
  REAL*8 lam1r,lam2r,lam3r,lammr,lampr,lamm_nobetar,lamp_nobetar
  REAL*8 lam1,lam2,lam3,lamm,lamp,lamp_tmp

  REAL*8 cs2l,cs2r,one,two
  REAL*8 vlowxr,vlowyr,vlowzr,v2r,wr
  REAL*8 vlowxl,vlowyl,vlowzl,v2l,wl

  REAL*8 pressl,dpdepsl,enthalpyl,kappal
  REAL*8 pressr,dpdepsr,enthalpyr,kappar
  REAL*8 axpl,axml,vxpl,vxml,xsil,dltl
  REAL*8 axpr,axmr,vxpr,vxmr,xsir,dltr
  REAL*8 cxx,cxy,cxz,gam,vxa,vxb
  INTEGER*4 handle
  REAL*8 invrhol, betainvalp, invfacl
  REAL*8 invrhor, invfacr
  REAL*8 invulpfac, invulmfac, invurpfac, invurmfac



  INTEGER*4  :: n 
  INTEGER*4  :: keytemp
  INTEGER*4  :: anyerr
  INTEGER*4  :: keyerr(1)
  REAL*8 :: xpress
  REAL*8 :: xeps
  REAL*8 :: xtemp
  REAL*8 :: xye

  n=1;anyerr=0;keyerr(1)=0
  xpress=0.0d0;xeps=0.0d0;xtemp=0.0d0;xye=0.0d0

  one = 1.0d0
  two = 2.0d0




  invrhol = one / rhol
  call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
       rhol,epsl,templ,ye_l,pressl,keyerr,anyerr)
  call EOS_Omni_cs2(handle,keytemp,GRHydro_eos_rf_prec,n,&
       rhol,epsl,templ,ye_l,cs2l,keyerr,anyerr)
  call EOS_Omni_DPressByDEps(handle,keytemp,GRHydro_eos_rf_prec,n,&
       rhol,epsl,templ,ye_l,dpdepsl,keyerr,anyerr)

  enthalpyl = one + epsl + pressl * invrhol

  vlowxl = gxx*velxl + gxy*velyl + gxz*velzl
  vlowyl = gxy*velxl + gyy*velyl + gyz*velzl
  vlowzl = gxz*velxl + gyz*velyl + gzz*velzl
  v2l = vlowxl*velxl + vlowyl*velyl + vlowzl*velzl

  if (v2l .ge. 1.d0-1.0d-10) then
    v2l = 1.d0-1.0d-10
  endif


  wl = one / sqrt(one - v2l)



  betainvalp = beta / alp
  lam1l = velxl - betainvalp
  lam2l = lam1l
  lam3l = lam1l
  lamp_tmp = cs2l*(one-v2l)*(u*(one-v2l*cs2l) - velxl**2*(one-cs2l))
  if (lamp_tmp .le. 1.0d-10) then
    lamp_tmp = 1.0d-10;
  endif
  invfacl = one/(one-v2l*cs2l)
  lamp_nobetal = (velxl*(one-cs2l) + sqrt(lamp_tmp))*invfacl
  lamm_nobetal = (velxl*(one-cs2l) - sqrt(lamp_tmp))*invfacl

  lampl = lamp_nobetal - betainvalp
  lamml = lamm_nobetal - betainvalp




  invrhor = one / rhor

  call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
       rhor,epsr,tempr,ye_r,pressr,keyerr,anyerr)
  call EOS_Omni_cs2(handle,keytemp,GRHydro_eos_rf_prec,n,&
       rhor,epsr,tempr,ye_r,cs2r,keyerr,anyerr)
  call EOS_Omni_DPressByDEps(handle,keytemp,GRHydro_eos_rf_prec,n,&
       rhor,epsr,tempr,ye_r,dpdepsr,keyerr,anyerr)

  enthalpyr = one + epsr + pressr * invrhor

  vlowxr = gxx*velxr + gxy*velyr + gxz*velzr
  vlowyr = gxy*velxr + gyy*velyr + gyz*velzr
  vlowzr = gxz*velxr + gyz*velyr + gzz*velzr
  v2r = vlowxr*velxr + vlowyr*velyr + vlowzr*velzr

  if (v2r .ge. 1.d0-1.0d-10) then
    v2r = 1.d0-1.0d-10
  endif

  wr = one / sqrt(one - v2r)



  lam1r = velxr - betainvalp
  lam2r = lam1r
  lam3r = lam1r
  lamp_tmp = cs2r*(one-v2r)*(u*(one-v2r*cs2r) - velxr**2*(one-cs2r))
  if (lamp_tmp .le. 1.0d-10) then
    lamp_tmp = 1.0d-10;
  endif
  invfacr = one/(one-v2r*cs2r)
  lamp_nobetar = (velxr*(one-cs2r) + sqrt(lamp_tmp))*invfacr
  lamm_nobetar = (velxr*(one-cs2r) - sqrt(lamp_tmp))*invfacr

  lampr = lamp_nobetar - betainvalp
  lammr = lamm_nobetar - betainvalp



  lam1 = max(abs(lam1l),abs(lam1r))
  lam2 = lam1
  lam3 = lam1
  lamp = max(abs(lampl),abs(lampr))
  lamm = max(abs(lamml),abs(lammr))







  lam(1) = lamm
  lam(5) = lam1
  lam(3) = lam2
  lam(4) = lam3
  lam(2) = lamp
  




  invulpfac = one / (u - velxl*lamp_nobetal)
  invulmfac = one / (u - velxl*lamm_nobetal)
  axpl = (u - velxl*velxl)*invulpfac
  axml = (u - velxl*velxl)*invulmfac
  vxpl = (velxl - lamp_nobetal)*invulpfac
  vxml = (velxl - lamm_nobetal)*invulmfac



  kappal = dpdepsl / (dpdepsl - rhol * cs2l)




  reivec1l(1) = kappal / (enthalpyl * wl)
  reivec1l(2) = vlowxl
  reivec1l(3) = vlowyl
  reivec1l(4) = vlowzl
  reivec1l(5) = one - reivec1l(1)



  reivec2l(1) = wl * vlowyl
  reivec2l(2) = enthalpyl * (gxy + two * wl * wl * vlowxl * vlowyl)
  reivec2l(3) = enthalpyl * (gyy + two * wl * wl * vlowyl * vlowyl)
  reivec2l(4) = enthalpyl * (gyz + two * wl * wl * vlowyl * vlowzl)
  reivec2l(5) = vlowyl * wl * (two * wl * enthalpyl - one)



  reivec3l(1) = wl * vlowzl
  reivec3l(2) = enthalpyl * (gxz + two * wl * wl * vlowxl * vlowzl)
  reivec3l(3) = enthalpyl * (gyz + two * wl * wl * vlowyl * vlowzl)
  reivec3l(4) = enthalpyl * (gzz + two * wl * wl * vlowzl * vlowzl)
  reivec3l(5) = vlowzl * wl * (two * wl * enthalpyl - one)



  reivecpl(1) = one
  reivecpl(2) = enthalpyl * wl * (vlowxl - vxpl)
  reivecpl(3) = enthalpyl * wl * vlowyl
  reivecpl(4) = enthalpyl * wl * vlowzl
  reivecpl(5) = enthalpyl * wl * axpl - one



  reivecml(1) = one
  reivecml(2) = enthalpyl * wl * (vlowxl - vxml)
  reivecml(3) = enthalpyl * wl * vlowyl
  reivecml(4) = enthalpyl * wl * vlowzl
  reivecml(5) = enthalpyl * wl * axml - one





  invurpfac = one / (u - velxr*lamp_nobetar)
  invurmfac = one / (u - velxr*lamm_nobetar)
  axpr = (u - velxr*velxr)*invurpfac
  axmr = (u - velxr*velxr)*invurmfac
  vxpr = (velxr - lamp_nobetar)*invurpfac
  vxmr = (velxr - lamm_nobetar)*invurmfac



  kappar = dpdepsr / (dpdepsr - rhor * cs2r)



  reivec1r(1) = kappar / (enthalpyr * wr)
  reivec1r(2) = vlowxr
  reivec1r(3) = vlowyr
  reivec1r(4) = vlowzr
  reivec1r(5) = one - reivec1r(1)



  reivec2r(1) = wr * vlowyr
  reivec2r(2) = enthalpyr * (gxy + two * wr * wr * vlowxr * vlowyr)
  reivec2r(3) = enthalpyr * (gyy + two * wr * wr * vlowyr * vlowyr)
  reivec2r(4) = enthalpyr * (gyz + two * wr * wr * vlowyr * vlowzr)
  reivec2r(5) = vlowyr * wr * (two * wr * enthalpyr - one)



  reivec3r(1) = wr * vlowzr
  reivec3r(2) = enthalpyr * (gxz + two * wr * wr * vlowxr * vlowzr)
  reivec3r(3) = enthalpyr * (gyz + two * wr * wr * vlowyr * vlowzr)
  reivec3r(4) = enthalpyr * (gzz + two * wr * wr * vlowzr * vlowzr)
  reivec3r(5) = vlowzr * wr * (two * wr * enthalpyr - one)



  reivecpr(1) = one
  reivecpr(2) = enthalpyr * wr * (vlowxr - vxpr)
  reivecpr(3) = enthalpyr * wr * vlowyr
  reivecpr(4) = enthalpyr * wr * vlowzr
  reivecpr(5) = enthalpyr * wr * axpr - one



  reivecmr(1) = one
  reivecmr(2) = enthalpyr * wr * (vlowxr - vxmr)
  reivecmr(3) = enthalpyr * wr * vlowyr
  reivecmr(4) = enthalpyr * wr * vlowzr
  reivecmr(5) = enthalpyr * wr * axmr - one



  if (ANALYTICAL) then

    cxx = gyy * gzz - gyz * gyz
    cxy = gxz * gyz - gxy * gzz
    cxz = gxy * gyz - gxz * gyy
    gam = gxx * cxx + gxy * cxy + gxz * cxz



    xsil = cxx - gam * velxl * velxl
    dltl = enthalpyl**3 * wl * (kappal - one) * (vxml - vxpl) * xsil



    tmp1 = wl / (kappal - one)

    leivec1l(1) = tmp1 * (enthalpyl - wl)
    leivec1l(2) = tmp1 * wl * velxl
    leivec1l(3) = tmp1 * wl * velyl
    leivec1l(4) = tmp1 * wl * velzl
    leivec1l(5) =-tmp1 * wl



    tmp1 = one / (xsil * enthalpyl)

    leivec2l(1) = (gyz * vlowzl - gzz * vlowyl) * tmp1
    leivec2l(2) = (gzz * vlowyl - gyz * vlowzl) * tmp1 * velxl
    leivec2l(3) = (gzz * (one - velxl * vlowxl) + gxz * vlowzl * velxl ) * tmp1
    leivec2l(4) = (gyz * (velxl * vlowxl - one) - gxz * velxl * vlowyl) * tmp1
    leivec2l(5) = (gyz * vlowzl - gzz * vlowyl) * tmp1



    leivec3l(1) = (gyz * vlowyl - gyy * vlowzl) * tmp1
    leivec3l(2) = (gyy * vlowzl - gyz * vlowyl) * tmp1 * velxl
    leivec3l(3) = (gyz * (velxl * vlowxl - one) - gxy * velxl * vlowzl) * tmp1
    leivec3l(4) = (gyy * (one - velxl * vlowxl) + gxy * velxl * vlowyl) * tmp1
    leivec3l(5) = (gyz * vlowyl - gyy * vlowzl) * tmp1



    tmp1 = enthalpyl * enthalpyl / dltl
    tmp2 = wl * wl * xsil

    leivecpl(1) = - (enthalpyl * wl * vxml * xsil + (one - kappal) * &
      (vxml * (tmp2 - cxx) - gam * velxl) - kappal * tmp2 * vxml) * tmp1
    leivecpl(2) = - (cxx * (one - kappal * axml) + (two * kappal - one) * &
      vxml * (tmp2 * velxl - cxx * velxl)) * tmp1
    leivecpl(3) = - (cxy * (one - kappal * axml) + (two * kappal - one) * &
      vxml * (tmp2 * velyl - cxy * velxl)) * tmp1
    leivecpl(4) = - (cxz * (one - kappal * axml) + (two * kappal - one) * &
      vxml * (tmp2 * velzl - cxz * velxl)) * tmp1
    leivecpl(5) = - ((one - kappal) * (vxml * (tmp2 - cxx) - gam * velxl) - &
      kappal * tmp2 * vxml) * tmp1



    leivecml(1) = (enthalpyl * wl * vxpl * xsil + (one - kappal) * &
      (vxpl * (tmp2 - cxx) - gam * velxl) - kappal * tmp2 * vxpl) * tmp1
    leivecml(2) = (cxx * (one - kappal * axpl) + (two * kappal - one) * &
      vxpl * (tmp2 * velxl - cxx * velxl)) * tmp1
    leivecml(3) = (cxy * (one - kappal * axpl) + (two * kappal - one) * &
      vxpl * (tmp2 * velyl - cxy * velxl)) * tmp1
    leivecml(4) = (cxz * (one - kappal * axpl) + (two * kappal - one) * &
      vxpl * (tmp2 * velzl - cxz * velxl)) * tmp1
    leivecml(5) = ((one - kappal) * (vxpl * (tmp2 - cxx) - gam * &
      velxl) - kappal * tmp2 * vxpl) * tmp1



    xsir = cxx - gam * velxr * velxr
    dltr = enthalpyr**3 * wr * (kappar - one) * (vxmr - vxpr) * xsir



    tmp1 = wr / (kappar - one)

    leivec1r(1) = tmp1 * (enthalpyr - wr)
    leivec1r(2) = tmp1 * wr * velxr
    leivec1r(3) = tmp1 * wr * velyr
    leivec1r(4) = tmp1 * wr * velzr
    leivec1r(5) =-tmp1 * wr



    tmp1 = one / (xsir * enthalpyr)

    leivec2r(1) = (gyz * vlowzr - gzz * vlowyr) * tmp1
    leivec2r(2) = (gzz * vlowyr - gyz * vlowzr) * tmp1 * velxr
    leivec2r(3) = (gzz * (one - velxr * vlowxr) + gxz * vlowzr * velxr) * tmp1
    leivec2r(4) = (gyz * (velxr * vlowxr - one) - gxz * velxr * vlowyr) * tmp1
    leivec2r(5) = (gyz * vlowzr - gzz * vlowyr) * tmp1



    leivec3r(1) = (gyz * vlowyr - gyy * vlowzr) * tmp1
    leivec3r(2) = (gyy * vlowzr - gyz * vlowyr) * tmp1 * velxr
    leivec3r(3) = (gyz * (velxr * vlowxr - one) - gxy * velxr * vlowzr) * tmp1
    leivec3r(4) = (gyy * (one - velxr * vlowxr) + gxy * velxr * vlowyr) * tmp1
    leivec3r(5) = (gyz * vlowyr - gyy * vlowzr) * tmp1



    tmp1 = enthalpyr * enthalpyr / dltr
    tmp2 = wr * wr * xsir

    leivecpr(1) = - (enthalpyr * wr * vxmr * xsir + (one - kappar) * &
      (vxmr * (tmp2 - cxx) - gam * velxr) - kappar * tmp2 * vxmr) * tmp1
    leivecpr(2) = - (cxx * (one - kappar * axmr) + (two * kappar - one) * &
      vxmr * (tmp2 * velxr - cxx * velxr)) * tmp1
    leivecpr(3) = - (cxy * (one - kappar * axmr) + (two * kappar - one) * &
      vxmr * (tmp2 * velyr - cxy * velxr)) * tmp1
    leivecpr(4) = - (cxz * (one - kappar * axmr) + (two * kappar - one) * &
      vxmr * (tmp2 * velzr - cxz * velxr)) * tmp1
    leivecpr(5) = - ((one - kappar) * (vxmr * (tmp2 - cxx) - gam * velxr) - &
      kappar * tmp2 * vxmr) * tmp1



    leivecmr(1) = (enthalpyr * wr * vxpr * xsir + (one - kappar) * &
      (vxpr * (tmp2 - cxx) - gam * velxr) - kappar * tmp2 * vxpr) * tmp1
    leivecmr(2) = (cxx * (one - kappar * axpr) + (two * kappar - one) * &
      vxpr * (tmp2 * velxr - cxx * velxr)) * tmp1
    leivecmr(3) = (cxy * (one - kappar * axpr) + (two * kappar - one) * &
      vxpr * (tmp2 * velyr - cxy * velxr)) * tmp1
    leivecmr(4) = (cxz * (one - kappar * axpr) + (two * kappar - one) * &
      vxpr * (tmp2 * velzr - cxz * velxr)) * tmp1
    leivecmr(5) = ((one - kappar) * (vxpr * (tmp2 - cxx) - gam * &
      velxr) - kappar * tmp2 * vxpr) * tmp1

  endif










  p(:,1) = reivecml(:)
  p(:,2) = reivecpl(:)
  p(:,3) = reivec2l(:)
  p(:,4) = reivec3l(:)
  p(:,5) = reivec1l(:)



  du(1) = densl
  du(2) = sxl
  du(3) = syl
  du(4) = szl
  du(5) = taul

  if (ANALYTICAL) then

    if (FAST) then

      sump = 0.0d0
      summ = 0.0d0

      do i=1,5
         sump = sump + (lamp - lam1) * leivecpl(i) * du(i)
         summ = summ + (lamm - lam1) * leivecml(i) * du(i)
      enddo

      vxa = sump + summ
      vxb =-(sump * vxpl + summ * vxml)

      rfluxl(1) = lam1 * du(1) + vxa
      rfluxl(2) = lam1 * du(2) + enthalpyl * wl * (vlowxl * vxa + vxb)
      rfluxl(3) = lam1 * du(3) + enthalpyl * wl * (vlowyl * vxa)
      rfluxl(4) = lam1 * du(4) + enthalpyl * wl * (vlowzl * vxa)
      rfluxl(5) = lam1 * du(5) + enthalpyl * wl * (velxl  * vxb + vxa) - vxa

    else



      q(1,:) = leivecml(:)
      q(2,:) = leivecpl(:)
      q(3,:) = leivec2l(:)
      q(4,:) = leivec3l(:)
      q(5,:) = leivec1l(:)

      do i=1,5
        dw(i) = 0.0d0
        do j=1,5
          dw(i) = dw(i) + q(i,j) * du(j)
        enddo
      enddo

      do i = 1, 5
        rfluxl(i) = 0.d0
        do j = 1, 5
          rfluxl(i) = rfluxl(i) + p(i,j) * lam(j) * dw(j)
        end do
      end do

    endif

  else



    dw=du
    aa = p

    do i=1,5
      paug(:,i) = p(:,i)
    enddo

















    paug(1,6) = du(1)
    paug(2,6) = du(2)
    paug(3,6) = du(3)
    paug(4,6) = du(4)
    paug(5,6) = du(5)


    do i=1,5

      tmp1 = one / paug(i,i)
      do j=i,6
        paug(i,j) = paug(i,j)*tmp1
      enddo

      do j=i+1,5
        tmp1 = - (paug(j,i))
        do k=i,6
          paug(j,k) = paug(j,k) + tmp1*paug(i,k)
        enddo
      enddo
    enddo

    f_du(5) = paug(5,6)
    do i=4,1,-1
      f_du(i) = paug(i,6)
      do j=i+1,5
        f_du(i) = f_du(i) - paug(i,j)*f_du(j)
      enddo
    enddo

    dw(1) = f_du(1)
    dw(2) = f_du(2)
    dw(3) = f_du(3)
    dw(4) = f_du(4)
    dw(5) = f_du(5)

    do i = 1, 5
      rfluxl(i) = 0.d0
      do j = 1, 5
        rfluxl(i) = rfluxl(i) + p(i,j) * lam(j) * dw(j)
      end do
    end do

  endif

    flux1l = rfluxl(1)
    flux2l = rfluxl(2)
    flux3l = rfluxl(3)
    flux4l = rfluxl(4)
    flux5l = rfluxl(5)











  p(:,1) = reivecmr(:)
  p(:,2) = reivecpr(:)
  p(:,3) = reivec2r(:)
  p(:,4) = reivec3r(:)
  p(:,5) = reivec1r(:)

  du(1) = densr
  du(2) = sxr
  du(3) = syr
  du(4) = szr
  du(5) = taur

  if (ANALYTICAL) then

    if (FAST) then

      sump = 0.0d0
      summ = 0.0d0

      do i=1,5
         sump = sump + (lamp - lam1) * leivecpr(i) * du(i)
         summ = summ + (lamm - lam1) * leivecmr(i) * du(i)
      enddo

      vxa = sump + summ
      vxb =-(sump * vxpr + summ * vxmr)

      rfluxr(1) = lam1 * du(1) + vxa
      rfluxr(2) = lam1 * du(2) + enthalpyr * wr * (vlowxr * vxa + vxb)
      rfluxr(3) = lam1 * du(3) + enthalpyr * wr * (vlowyr * vxa)
      rfluxr(4) = lam1 * du(4) + enthalpyr * wr * (vlowzr * vxa)
      rfluxr(5) = lam1 * du(5) + enthalpyr * wr * (velxr  * vxb + vxa) - vxa

    else



      q(1,:) = leivecmr(:)
      q(2,:) = leivecpr(:)
      q(3,:) = leivec2r(:)
      q(4,:) = leivec3r(:)
      q(5,:) = leivec1r(:)

      do i=1,5
        dw(i) = 0.0d0
        do j=1,5
          dw(i) = dw(i) + q(i,j) * du(j)
        enddo
      enddo

      do i = 1, 5
        rfluxr(i) = 0.d0
        do j = 1, 5
          rfluxr(i) = rfluxr(i) + p(i,j) * lam(j) * dw(j)
        end do
      end do

    endif

  else



    do i=1,5
      dw(i) = du(i)
      do j=1,5
        aa(i,j) = p(i,j)
      end do
    enddo

    do i=1,5
      paug(i,1) = p(i,1)
      paug(i,2) = p(i,2)
      paug(i,3) = p(i,3)
      paug(i,4) = p(i,4)
      paug(i,5) = p(i,5)
    enddo

    paug(1,6) = du(1)
    paug(2,6) = du(2)
    paug(3,6) = du(3)
    paug(4,6) = du(4)
    paug(5,6) = du(5)



    do i=1,5

      tmp1 = one / paug(i,i)
      do j=i,6
        paug(i,j) = paug(i,j)*tmp1
      enddo

      do j=i+1,5
        tmp1 = - (paug(j,i))
        do k=i,6
          paug(j,k) = paug(j,k) + tmp1*paug(i,k)
        enddo
      enddo
    enddo


    f_du(5) = paug(5,6)
    do i=4,1,-1
      f_du(i) = paug(i,6)
      do j=i+1,5
        f_du(i) = f_du(i) - paug(i,j)*f_du(j)
      enddo
    enddo

    dw(1) = f_du(1)
    dw(2) = f_du(2)
    dw(3) = f_du(3)
    dw(4) = f_du(4)
    dw(5) = f_du(5)

    do i = 1, 5
      rfluxr(i) = 0.d0
      do j = 1, 5
        rfluxr(i) = rfluxr(i) + p(i,j) * lam(j) * dw(j)
      end do
    end do

  endif

  flux1r = rfluxr(1)
  flux2r = rfluxr(2)
  flux3r = rfluxr(3)
  flux4r = rfluxr(4)
  flux5r = rfluxr(5)

  flux1 = flux1r - flux1l 
  flux2 = flux2r - flux2l 
  flux3 = flux3r - flux3l 
  flux4 = flux4r - flux4l 
  flux5 = flux5r - flux5l 

  return

end subroutine eigenproblem_marquina_hot

