














































 














































































 


































































































































































 

















































































































































































































































































































































































































































 






































































































 










 








































































































































































































































 










                                
                                

                                
                                

                                
                                

                                



















































































































 

subroutine SimplePPM_1dM(handle,poly,nx,dx,&
     rho,velx,vely,velz,Bvcx,Bvcy,Bvcz,psidc,eps,press,&
     rhominus,velxminus,velyminus,velzminus,Bvcxminus,Bvcyminus,Bvczminus,psidcminus,epsminus,&
     rhoplus,velxplus,velyplus,velzplus,Bvcxplus,Bvcyplus,Bvczplus,psidcplus,epsplus,dc_flag,&
     trivial_rp, hydro_excision_mask,&
     gxx, gxy, gxz, gyy, gyz, gzz, beta, alp, w_lorentz, &
     dir, ni, nj, nrx, nry, nrz, ev_l, ev_r, xw)

  USE GRHydro_Scalars
  USE GRHydro_EigenproblemM

  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
  interface
  integer*8 function CCTK_PointerTo (var)
  implicit none
  type(*),dimension(..),TARGET :: var
  end function CCTK_PointerTo
  end interface
  interface
  integer function CCTK_Equals (arg1, arg2)
  implicit none
  integer*8 arg1
  character(*) arg2
  end function CCTK_Equals
  integer function CCTK_MyProc (cctkGH)
  implicit none
  integer*8 cctkGH
  end function CCTK_MyProc
  integer function CCTK_nProcs (cctkGH)
  implicit none
  integer*8 cctkGH
  end function CCTK_nProcs
  integer function CCTK_IsThornActive (name)
  implicit none
  character(*) name
  end function CCTK_IsThornActive
  integer*8 function CCTK_NullPointer ()
  implicit none
  end function CCTK_NullPointer
  end interface
  interface
  INTEGER*4 function Boundary_SelectGroupForBC (GH, faces, boundary_width, table_handle, group_name, bc_name)
  implicit none
  integer*8 GH
  INTEGER*4 faces
  INTEGER*4 boundary_width
  INTEGER*4 table_handle
  character(*) group_name
  character(*) bc_name
  end function Boundary_SelectGroupForBC
  end interface
  interface
  INTEGER*4 function Boundary_SelectVarForBC (GH, faces, boundary_width, table_handle, var_name, bc_name)
  implicit none
  integer*8 GH
  INTEGER*4 faces
  INTEGER*4 boundary_width
  INTEGER*4 table_handle
  character(*) var_name
  character(*) bc_name
  end function Boundary_SelectVarForBC
  end interface
  interface
  subroutine Con2PrimGen (handle, dens, sx, sy, sz, tau, rho, velx, vely, velz, epsilon, pressure, w_lorentz, uxx, uxy, uxz, uyy, u&
  &yz, uzz, det, x, y, z, r, epsnegative, GRHydro_rho_min, pmin, epsmin, GRHydro_reflevel, retval)
  implicit none
  INTEGER*4 handle
  REAL*8 dens
  REAL*8 sx
  REAL*8 sy
  REAL*8 sz
  REAL*8 tau
  REAL*8 rho
  REAL*8 velx
  REAL*8 vely
  REAL*8 velz
  REAL*8 epsilon
  REAL*8 pressure
  REAL*8 w_lorentz
  REAL*8 uxx
  REAL*8 uxy
  REAL*8 uxz
  REAL*8 uyy
  REAL*8 uyz
  REAL*8 uzz
  REAL*8 det
  REAL*8 x
  REAL*8 y
  REAL*8 z
  REAL*8 r
  INTEGER*4 epsnegative
  REAL*8 GRHydro_rho_min
  REAL*8 pmin
  REAL*8 epsmin
  INTEGER*4 GRHydro_reflevel
  REAL*8 retval
  end subroutine Con2PrimGen
  end interface
  interface
  subroutine Con2PrimGenM (handle, keytemp, prec, gamma_eos, dens, sx, sy, sz, tau, Bconsx, Bconsy, Bconsz, y_e, temp, rho, velx, v&
  &ely, velz, epsilon, pressure, Bvecx, Bvecy, Bvecz, bvecsq, w_lorentz, gxx, gxy, gxz, gyy, gyz, gzz, uxx, uxy, uxz, uyy, uyz, uzz&
  &, det, epsnegative, retval)
  implicit none
  INTEGER*4 handle
  INTEGER*4 keytemp
  REAL*8 prec
  REAL*8 gamma_eos
  REAL*8 dens
  REAL*8 sx
  REAL*8 sy
  REAL*8 sz
  REAL*8 tau
  REAL*8 Bconsx
  REAL*8 Bconsy
  REAL*8 Bconsz
  REAL*8 y_e
  REAL*8 temp
  REAL*8 rho
  REAL*8 velx
  REAL*8 vely
  REAL*8 velz
  REAL*8 epsilon
  REAL*8 pressure
  REAL*8 Bvecx
  REAL*8 Bvecy
  REAL*8 Bvecz
  REAL*8 bvecsq
  REAL*8 w_lorentz
  REAL*8 gxx
  REAL*8 gxy
  REAL*8 gxz
  REAL*8 gyy
  REAL*8 gyz
  REAL*8 gzz
  REAL*8 uxx
  REAL*8 uxy
  REAL*8 uxz
  REAL*8 uyy
  REAL*8 uyz
  REAL*8 uzz
  REAL*8 det
  INTEGER*4 epsnegative
  REAL*8 retval
  end subroutine Con2PrimGenM
  end interface
  interface
  subroutine Con2PrimGenMee (handle, keytemp, prec, gamma_eos, dens, sx, sy, sz, tau, Bconsx, Bconsy, Bconsz, entropycons, y_e, tem&
  &p, rho, velx, vely, velz, epsilon, pressure, Bvecx, Bvecy, Bvecz, Bvecsq, w_lorentz, gxx, gxy, gxz, gyy, gyz, gzz, uxx, uxy, uxz&
  &, uyy, uyz, uzz, det, epsnegative, retval)
  implicit none
  INTEGER*4 handle
  INTEGER*4 keytemp
  REAL*8 prec
  REAL*8 gamma_eos
  REAL*8 dens
  REAL*8 sx
  REAL*8 sy
  REAL*8 sz
  REAL*8 tau
  REAL*8 Bconsx
  REAL*8 Bconsy
  REAL*8 Bconsz
  REAL*8 entropycons
  REAL*8 y_e
  REAL*8 temp
  REAL*8 rho
  REAL*8 velx
  REAL*8 vely
  REAL*8 velz
  REAL*8 epsilon
  REAL*8 pressure
  REAL*8 Bvecx
  REAL*8 Bvecy
  REAL*8 Bvecz
  REAL*8 Bvecsq
  REAL*8 w_lorentz
  REAL*8 gxx
  REAL*8 gxy
  REAL*8 gxz
  REAL*8 gyy
  REAL*8 gyz
  REAL*8 gzz
  REAL*8 uxx
  REAL*8 uxy
  REAL*8 uxz
  REAL*8 uyy
  REAL*8 uyz
  REAL*8 uzz
  REAL*8 det
  INTEGER*4 epsnegative
  REAL*8 retval
  end subroutine Con2PrimGenMee
  end interface
  interface
  subroutine Con2PrimPoly (handle, dens, sx, sy, sz, tau, rho, velx, vely, velz, epsilon, press, w_lorentz, uxx, uxy, uxz, uyy, uyz&
  &, uzz, det, x, y, z, r, rho_min, GRHydro_reflevel, GRHydro_C2P_failed)
  implicit none
  INTEGER*4 handle
  REAL*8 dens
  REAL*8 sx
  REAL*8 sy
  REAL*8 sz
  REAL*8 tau
  REAL*8 rho
  REAL*8 velx
  REAL*8 vely
  REAL*8 velz
  REAL*8 epsilon
  REAL*8 press
  REAL*8 w_lorentz
  REAL*8 uxx
  REAL*8 uxy
  REAL*8 uxz
  REAL*8 uyy
  REAL*8 uyz
  REAL*8 uzz
  REAL*8 det
  REAL*8 x
  REAL*8 y
  REAL*8 z
  REAL*8 r
  REAL*8 rho_min
  INTEGER*4 GRHydro_reflevel
  REAL*8 GRHydro_C2P_failed
  end subroutine Con2PrimPoly
  end interface
  interface
  subroutine Con2PrimPolyM (handle, gamma_eos, dens, sx, sy, sz, sc, Bconsx, Bconsy, Bconsz, rho, velx, vely, velz, epsilon, pressu&
  &re, Bvecx, Bvecy, Bvecz, Bvecsq, w_lorentz, gxx, gxy, gxz, gyy, gyz, gzz, uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative, retval)
  implicit none
  INTEGER*4 handle
  REAL*8 gamma_eos
  REAL*8 dens
  REAL*8 sx
  REAL*8 sy
  REAL*8 sz
  REAL*8 sc
  REAL*8 Bconsx
  REAL*8 Bconsy
  REAL*8 Bconsz
  REAL*8 rho
  REAL*8 velx
  REAL*8 vely
  REAL*8 velz
  REAL*8 epsilon
  REAL*8 pressure
  REAL*8 Bvecx
  REAL*8 Bvecy
  REAL*8 Bvecz
  REAL*8 Bvecsq
  REAL*8 w_lorentz
  REAL*8 gxx
  REAL*8 gxy
  REAL*8 gxz
  REAL*8 gyy
  REAL*8 gyz
  REAL*8 gzz
  REAL*8 uxx
  REAL*8 uxy
  REAL*8 uxz
  REAL*8 uyy
  REAL*8 uyz
  REAL*8 uzz
  REAL*8 det
  INTEGER*4 epsnegative
  REAL*8 retval
  end subroutine Con2PrimPolyM
  end interface
  interface
  subroutine EOS_Omni_DPressByDEps (eoskey, havetemp, rf_precision, npoints, rho, eps, temp, ye, dpdepsrho, keyerr, anyerr)
  implicit none
  INTEGER*4 eoskey
  INTEGER*4 havetemp
  REAL*8 rf_precision
  INTEGER*4 npoints
  REAL*8 rho(*)
  REAL*8 eps(*)
  REAL*8 temp(*)
  REAL*8 ye(*)
  REAL*8 dpdepsrho(*)
  INTEGER*4 keyerr(*)
  INTEGER*4 anyerr
  end subroutine EOS_Omni_DPressByDEps
  end interface
  interface
  subroutine EOS_Omni_DPressByDRho (eoskey, havetemp, rf_precision, npoints, rho, eps, temp, ye, dpdrhoeps, keyerr, anyerr)
  implicit none
  INTEGER*4 eoskey
  INTEGER*4 havetemp
  REAL*8 rf_precision
  INTEGER*4 npoints
  REAL*8 rho(*)
  REAL*8 eps(*)
  REAL*8 temp(*)
  REAL*8 ye(*)
  REAL*8 dpdrhoeps(*)
  INTEGER*4 keyerr(*)
  INTEGER*4 anyerr
  end subroutine EOS_Omni_DPressByDRho
  end interface
  interface
  subroutine EOS_Omni_EpsFromPress (eoskey, havetemp, rf_precision, npoints, rho, eps, temp, ye, press, xeps, keyerr, anyerr)
  implicit none
  INTEGER*4 eoskey
  INTEGER*4 havetemp
  REAL*8 rf_precision
  INTEGER*4 npoints
  REAL*8 rho(*)
  REAL*8 eps(*)
  REAL*8 temp(*)
  REAL*8 ye(*)
  REAL*8 press(*)
  REAL*8 xeps(*)
  INTEGER*4 keyerr(*)
  INTEGER*4 anyerr
  end subroutine EOS_Omni_EpsFromPress
  end interface
  interface
  INTEGER*4 function EOS_Omni_GetHandle (name)
  implicit none
  character(*) name
  end function EOS_Omni_GetHandle
  end interface
  interface
  subroutine EOS_Omni_cs2 (eoskey, havetemp, rf_precision, npoints, rho, eps, temp, ye, cs2, keyerr, anyerr)
  implicit none
  INTEGER*4 eoskey
  INTEGER*4 havetemp
  REAL*8 rf_precision
  INTEGER*4 npoints
  REAL*8 rho(*)
  REAL*8 eps(*)
  REAL*8 temp(*)
  REAL*8 ye(*)
  REAL*8 cs2(*)
  INTEGER*4 keyerr(*)
  INTEGER*4 anyerr
  end subroutine EOS_Omni_cs2
  end interface
  interface
  subroutine EOS_Omni_dpderho_dpdrhoe (eoskey, havetemp, rf_precision, npoints, rho, eps, temp, ye, dpderho, dpdrhoe, keyerr, anyer&
  &r)
  implicit none
  INTEGER*4 eoskey
  INTEGER*4 havetemp
  REAL*8 rf_precision
  INTEGER*4 npoints
  REAL*8 rho(*)
  REAL*8 eps(*)
  REAL*8 temp(*)
  REAL*8 ye(*)
  REAL*8 dpderho(*)
  REAL*8 dpdrhoe(*)
  INTEGER*4 keyerr(*)
  INTEGER*4 anyerr
  end subroutine EOS_Omni_dpderho_dpdrhoe
  end interface
  interface
  subroutine EOS_Omni_press (eoskey, havetemp, rf_precision, npoints, rho, eps, temp, ye, press, keyerr, anyerr)
  implicit none
  INTEGER*4 eoskey
  INTEGER*4 havetemp
  REAL*8 rf_precision
  INTEGER*4 npoints
  REAL*8 rho(*)
  REAL*8 eps(*)
  REAL*8 temp(*)
  REAL*8 ye(*)
  REAL*8 press(*)
  INTEGER*4 keyerr(*)
  INTEGER*4 anyerr
  end subroutine EOS_Omni_press
  end interface
  interface
  subroutine EOS_Omni_pressOMP (eoskey, havetemp, rf_precision, npoints, rho, eps, temp, ye, press, keyerr, anyerr)
  implicit none
  INTEGER*4 eoskey
  INTEGER*4 havetemp
  REAL*8 rf_precision
  INTEGER*4 npoints
  REAL*8 rho(*)
  REAL*8 eps(*)
  REAL*8 temp(*)
  REAL*8 ye(*)
  REAL*8 press(*)
  INTEGER*4 keyerr(*)
  INTEGER*4 anyerr
  end subroutine EOS_Omni_pressOMP
  end interface
  interface
  subroutine EOS_Omni_press_cs2 (eoskey, havetemp, rf_precision, npoints, rho, eps, temp, ye, press, cs2, keyerr, anyerr)
  implicit none
  INTEGER*4 eoskey
  INTEGER*4 havetemp
  REAL*8 rf_precision
  INTEGER*4 npoints
  REAL*8 rho(*)
  REAL*8 eps(*)
  REAL*8 temp(*)
  REAL*8 ye(*)
  REAL*8 press(*)
  REAL*8 cs2(*)
  INTEGER*4 keyerr(*)
  INTEGER*4 anyerr
  end subroutine EOS_Omni_press_cs2
  end interface
  interface
  INTEGER*4 function MoLRegisterConstrained (ConstrainedIndex)
  implicit none
  INTEGER*4 ConstrainedIndex
  end function MoLRegisterConstrained
  end interface
  interface
  INTEGER*4 function MoLRegisterConstrainedGroup (ConstrainedIndex)
  implicit none
  INTEGER*4 ConstrainedIndex
  end function MoLRegisterConstrainedGroup
  end interface
  interface
  INTEGER*4 function MoLRegisterEvolved (EvolvedIndex, RHSIndex)
  implicit none
  INTEGER*4 EvolvedIndex
  INTEGER*4 RHSIndex
  end function MoLRegisterEvolved
  end interface
  interface
  INTEGER*4 function MoLRegisterEvolvedGroup (EvolvedIndex, RHSIndex)
  implicit none
  INTEGER*4 EvolvedIndex
  INTEGER*4 RHSIndex
  end function MoLRegisterEvolvedGroup
  end interface
  interface
  INTEGER*4 function MoLRegisterEvolvedGroupSlow (EvolvedIndex, RHSIndexSlow)
  implicit none
  INTEGER*4 EvolvedIndex
  INTEGER*4 RHSIndexSlow
  end function MoLRegisterEvolvedGroupSlow
  end interface
  interface
  INTEGER*4 function MoLRegisterEvolvedSlow (EvolvedIndex, RHSIndexSlow)
  implicit none
  INTEGER*4 EvolvedIndex
  INTEGER*4 RHSIndexSlow
  end function MoLRegisterEvolvedSlow
  end interface
  interface
  INTEGER*4 function MoLRegisterSaveAndRestoreGroup (SandRIndex)
  implicit none
  INTEGER*4 SandRIndex
  end function MoLRegisterSaveAndRestoreGroup
  end interface
  interface
  subroutine Prim2ConGen (handle, gxx, gxy, gxz, gyy, gyz, gzz, det, dens, sx, sy, sz, tau, rho, velx, vely, velz, epsilon, press, &
  &w_lorentz)
  implicit none
  INTEGER*4 handle
  REAL*8 gxx
  REAL*8 gxy
  REAL*8 gxz
  REAL*8 gyy
  REAL*8 gyz
  REAL*8 gzz
  REAL*8 det
  REAL*8 dens
  REAL*8 sx
  REAL*8 sy
  REAL*8 sz
  REAL*8 tau
  REAL*8 rho
  REAL*8 velx
  REAL*8 vely
  REAL*8 velz
  REAL*8 epsilon
  REAL*8 press
  REAL*8 w_lorentz
  end subroutine Prim2ConGen
  end interface
  interface
  subroutine Prim2ConGenM (handle, gxx, gxy, gxz, gyy, gyz, gzz, det, dens, sx, sy, sz, tau, Bconsx, Bconsy, Bconsz, rho, velx, vel&
  &y, velz, epsilon, press, Bvecx, Bvecy, Bvecz, w_lorentz)
  implicit none
  INTEGER*4 handle
  REAL*8 gxx
  REAL*8 gxy
  REAL*8 gxz
  REAL*8 gyy
  REAL*8 gyz
  REAL*8 gzz
  REAL*8 det
  REAL*8 dens
  REAL*8 sx
  REAL*8 sy
  REAL*8 sz
  REAL*8 tau
  REAL*8 Bconsx
  REAL*8 Bconsy
  REAL*8 Bconsz
  REAL*8 rho
  REAL*8 velx
  REAL*8 vely
  REAL*8 velz
  REAL*8 epsilon
  REAL*8 press
  REAL*8 Bvecx
  REAL*8 Bvecy
  REAL*8 Bvecz
  REAL*8 w_lorentz
  end subroutine Prim2ConGenM
  end interface
  interface
  subroutine Prim2ConGenM_hot (handle, GRHydro_reflevel, i, j, k, x, y, z, gxx, gxy, gxz, gyy, gyz, gzz, det, dens, sx, sy, sz, tau&
  &, Bconsx, Bconsy, Bconsz, rho, velx, vely, velz, epsilon, press, Bvecx, Bvecy, Bvecz, w_lorentz, temperature, Y_e)
  implicit none
  INTEGER*4 handle
  INTEGER*4 GRHydro_reflevel
  INTEGER*4 i
  INTEGER*4 j
  INTEGER*4 k
  REAL*8 x
  REAL*8 y
  REAL*8 z
  REAL*8 gxx
  REAL*8 gxy
  REAL*8 gxz
  REAL*8 gyy
  REAL*8 gyz
  REAL*8 gzz
  REAL*8 det
  REAL*8 dens
  REAL*8 sx
  REAL*8 sy
  REAL*8 sz
  REAL*8 tau
  REAL*8 Bconsx
  REAL*8 Bconsy
  REAL*8 Bconsz
  REAL*8 rho
  REAL*8 velx
  REAL*8 vely
  REAL*8 velz
  REAL*8 epsilon
  REAL*8 press
  REAL*8 Bvecx
  REAL*8 Bvecy
  REAL*8 Bvecz
  REAL*8 w_lorentz
  REAL*8 temperature
  REAL*8 Y_e
  end subroutine Prim2ConGenM_hot
  end interface
  interface
  subroutine Prim2ConPoly (handle, gxx, gxy, gxz, gyy, gyz, gzz, det, dens, sx, sy, sz, tau, rho, velx, vely, velz, epsilon, press,&
  & w_lorentz)
  implicit none
  INTEGER*4 handle
  REAL*8 gxx
  REAL*8 gxy
  REAL*8 gxz
  REAL*8 gyy
  REAL*8 gyz
  REAL*8 gzz
  REAL*8 det
  REAL*8 dens
  REAL*8 sx
  REAL*8 sy
  REAL*8 sz
  REAL*8 tau
  REAL*8 rho
  REAL*8 velx
  REAL*8 vely
  REAL*8 velz
  REAL*8 epsilon
  REAL*8 press
  REAL*8 w_lorentz
  end subroutine Prim2ConPoly
  end interface
  interface
  subroutine Prim2ConPolyM (handle, gxx, gxy, gxz, gyy, gyz, gzz, det, dens, sx, sy, sz, tau, Bconsx, Bconsy, Bconsz, rho, velx, ve&
  &ly, velz, epsilon, press, Bvecx, Bvecy, Bvecz, w_lorentz)
  implicit none
  INTEGER*4 handle
  REAL*8 gxx
  REAL*8 gxy
  REAL*8 gxz
  REAL*8 gyy
  REAL*8 gyz
  REAL*8 gzz
  REAL*8 det
  REAL*8 dens
  REAL*8 sx
  REAL*8 sy
  REAL*8 sz
  REAL*8 tau
  REAL*8 Bconsx
  REAL*8 Bconsy
  REAL*8 Bconsz
  REAL*8 rho
  REAL*8 velx
  REAL*8 vely
  REAL*8 velz
  REAL*8 epsilon
  REAL*8 press
  REAL*8 Bvecx
  REAL*8 Bvecy
  REAL*8 Bvecz
  REAL*8 w_lorentz
  end subroutine Prim2ConPolyM
  end interface
  interface
  subroutine SetDriftCorrectPosition (GH_Pointer, x, y, z)
  implicit none
  integer*8 GH_Pointer
  REAL*8 x
  REAL*8 y
  REAL*8 z
  end subroutine SetDriftCorrectPosition
  end interface
  interface
  subroutine SpatialDet (gxx, gxy, gxz, gyy, gyz, gzz, det)
  implicit none
  REAL*8 gxx
  REAL*8 gxy
  REAL*8 gxz
  REAL*8 gyy
  REAL*8 gyz
  REAL*8 gzz
  REAL*8 det
  end subroutine SpatialDet
  end interface
  interface
  subroutine UpperMet (uxx, uxy, uxz, uyy, uyz, uzz, det, gxx, gxy, gxz, gyy, gyz, gzz)
  implicit none
  REAL*8 uxx
  REAL*8 uxy
  REAL*8 uxz
  REAL*8 uyy
  REAL*8 uyz
  REAL*8 uzz
  REAL*8 det
  REAL*8 gxx
  REAL*8 gxy
  REAL*8 gxz
  REAL*8 gyy
  REAL*8 gyz
  REAL*8 gzz
  end subroutine UpperMet
  end interface

  INTEGER*4 :: handle,poly,nx
  REAL*8 :: dx
  REAL*8, dimension(nx) :: rho,velx,vely,velz,eps
  REAL*8, dimension(nx) :: Bvcx,Bvcy,Bvcz,psidc
  REAL*8, dimension(nx) :: rhominus,velxminus,velyminus,velzminus,epsminus
  REAL*8, dimension(nx) :: Bvcxminus,Bvcyminus,Bvczminus,psidcminus
  REAL*8, dimension(nx) :: rhoplus,velxplus,velyplus,velzplus,epsplus
  REAL*8, dimension(nx) :: Bvcxplus,Bvcyplus,Bvczplus,psidcplus

  INTEGER*4 :: i,s,dc_flag
  REAL*8, dimension(nx) :: drho,dvelx,dvely,dvelz,deps
  REAL*8, dimension(nx) :: dBvcx,dBvcy,dBvcz,dpsidc
  REAL*8, dimension(nx) :: dmrho,dmvelx,dmvely,dmvelz,dmeps
  REAL*8, dimension(nx) :: dmBvcx,dmBvcy,dmBvcz,dmpsidc
  REAL*8, dimension(nx) :: press,dpress,d2rho,tilde_flatten
  REAL*8 :: dpress2,dvel,w,flatten,eta,etatilde

  logical, dimension(nx) :: trivial_rp

  INTEGER*4, dimension(nx) :: hydro_excision_mask

  REAL*8, dimension(nx) :: gxx, gxy, gxz, gyy, gyz, gzz, &
                              beta, alp, w_lorentz
  INTEGER*4 :: dir, ni, nj, nrx, nry, nrz
  REAL*8, dimension(nrx, nry, nrz) :: ev_l, ev_r, xw

  REAL*8 :: uxx, uxy, uxz, uyy, uyz, uzz, det
  REAL*8, dimension(5) :: lam
  REAL*8 :: agxx, agxy, agxz, agyy, agyz, agzz
  REAL*8, dimension(nx) :: xwind, l_ev_l, l_ev_r

  logical :: cond

  REAL*8, dimension(nx) :: vx,vy,vz
  
  
  
  
  if (reconstruct_Wv.ne.0) then
    
    
    
    vx = w_lorentz*velx
    vy = w_lorentz*vely
    vz = w_lorentz*velz
  else
    vx = velx
    vy = vely
    vz = velz
  end if




trivial_rp = .true.




  do i = 2, nx - 1
    drho(i) = 0.5d0 * (rho(i+1) - rho(i-1))
    dvelx(i) = 0.5d0 * (vx(i+1) - vx(i-1))
    dvely(i) = 0.5d0 * (vy(i+1) - vy(i-1))
    dvelz(i) = 0.5d0 * (vz(i+1) - vz(i-1))
    dBvcx(i) = 0.5d0 * (Bvcx(i+1) - Bvcx(i-1))
    dBvcy(i) = 0.5d0 * (Bvcy(i+1) - Bvcy(i-1))
    dBvcz(i) = 0.5d0 * (Bvcz(i+1) - Bvcz(i-1))
    if(dc_flag.ne.0)dpsidc(i) = 0.5d0 * (psidc(i+1) - psidc(i-1))

    dpress(i) = press(i+1) - press(i-1)
    d2rho(i) = (rho(i+1) - 2.d0 * rho(i) + rho(i-1))
    
    
  end do
  if (poly .eq. 0) then
    do i = 2, nx - 1
      deps(i) = 0.5d0 * (eps(i+1) - eps(i-1))
    end do
  end if



  do i = 2, nx - 1
    if ( (rho(i+1) - rho(i)) * (rho(i) - rho(i-1)) > 0.d0 ) then
  dmrho(i) = sign(1.d0,  drho(i)) *                                              min(abs( drho(i)), 2.d0 * abs(rho(i) - rho(i-1)), &
  &                                          2.d0 * abs(rho(i+1) - rho(i)))
  else
  dmrho(i) = 0.d0
  end if
    if ( (vx(i+1) - vx(i)) * (vx(i) - vx(i-1)) > 0.d0 ) then
  dmvelx(i) = sign(1.d0,  dvelx(i)) *                                              min(abs( dvelx(i)), 2.d0 * abs(vx(i) - vx(i-1)),&
  &                                           2.d0 * abs(vx(i+1) - vx(i)))
  else
  dmvelx(i) = 0.d0
  end if
    if ( (vy(i+1) - vy(i)) * (vy(i) - vy(i-1)) > 0.d0 ) then
  dmvely(i) = sign(1.d0,  dvely(i)) *                                              min(abs( dvely(i)), 2.d0 * abs(vy(i) - vy(i-1)),&
  &                                           2.d0 * abs(vy(i+1) - vy(i)))
  else
  dmvely(i) = 0.d0
  end if
    if ( (vz(i+1) - vz(i)) * (vz(i) - vz(i-1)) > 0.d0 ) then
  dmvelz(i) = sign(1.d0,  dvelz(i)) *                                              min(abs( dvelz(i)), 2.d0 * abs(vz(i) - vz(i-1)),&
  &                                           2.d0 * abs(vz(i+1) - vz(i)))
  else
  dmvelz(i) = 0.d0
  end if
    if ( (Bvcx(i+1) - Bvcx(i)) * (Bvcx(i) - Bvcx(i-1)) > 0.d0 ) then
  dmBvcx(i) = sign(1.d0,  dBvcx(i)) *                                              min(abs( dBvcx(i)), 2.d0 * abs(Bvcx(i) - Bvcx(i-&
  &1)),                                           2.d0 * abs(Bvcx(i+1) - Bvcx(i)))
  else
  dmBvcx(i) = 0.d0
  end if
    if ( (Bvcy(i+1) - Bvcy(i)) * (Bvcy(i) - Bvcy(i-1)) > 0.d0 ) then
  dmBvcy(i) = sign(1.d0,  dBvcy(i)) *                                              min(abs( dBvcy(i)), 2.d0 * abs(Bvcy(i) - Bvcy(i-&
  &1)),                                           2.d0 * abs(Bvcy(i+1) - Bvcy(i)))
  else
  dmBvcy(i) = 0.d0
  end if
    if ( (Bvcz(i+1) - Bvcz(i)) * (Bvcz(i) - Bvcz(i-1)) > 0.d0 ) then
  dmBvcz(i) = sign(1.d0,  dBvcz(i)) *                                              min(abs( dBvcz(i)), 2.d0 * abs(Bvcz(i) - Bvcz(i-&
  &1)),                                           2.d0 * abs(Bvcz(i+1) - Bvcz(i)))
  else
  dmBvcz(i) = 0.d0
  end if
    if(dc_flag.ne.0) then
       if ( (psidc(i+1) - psidc(i)) * (psidc(i) - psidc(i-1)) > 0.d0 ) then
  dmpsidc(i) = sign(1.d0, dpsidc(i)) *                                              min(abs(dpsidc(i)), 2.d0 * abs(psidc(i) - psidc&
  &(i-1)),                                           2.d0 * abs(psidc(i+1) - psidc(i)))
  else
  dmpsidc(i) = 0.d0
  end if
    endif

  end do
  if (poly .eq. 0) then
    do i = 2, nx - 1
      if ( (eps(i+1) - eps(i)) * (eps(i) - eps(i-1)) > 0.d0 ) then
  dmeps(i) = sign(1.d0,  deps(i)) *                                              min(abs( deps(i)), 2.d0 * abs(eps(i) - eps(i-1)), &
  &                                          2.d0 * abs(eps(i+1) - eps(i)))
  else
  dmeps(i) = 0.d0
  end if
    end do
  end if



  do i = 2, nx-2
    rhoplus(i) = 0.5d0 * (rho(i) + rho(i+1)) + &
         (dmrho(i) - dmrho(i+1)) / 6.d0
    rhominus(i+1) = rhoplus(i)
    velxplus(i) = 0.5d0 * (vx(i) + vx(i+1)) + &
         (dmvelx(i) - dmvelx(i+1)) / 6.d0
    velxminus(i+1) = velxplus(i)
    velyplus(i) = 0.5d0 * (vy(i) + vy(i+1)) + &
         (dmvely(i) - dmvely(i+1)) / 6.d0
    velyminus(i+1) = velyplus(i)
    velzplus(i) = 0.5d0 * (vz(i) + vz(i+1)) + &
         (dmvelz(i) - dmvelz(i+1)) / 6.d0
    velzminus(i+1) = velzplus(i)

    Bvcxplus(i) = 0.5d0 * (Bvcx(i) + Bvcx(i+1)) + &
         (dmBvcx(i) - dmBvcx(i+1)) / 6.d0
    Bvcxminus(i+1) = Bvcxplus(i)
    Bvcyplus(i) = 0.5d0 * (Bvcy(i) + Bvcy(i+1)) + &
         (dmBvcy(i) - dmBvcy(i+1)) / 6.d0
    Bvcyminus(i+1) = Bvcyplus(i)
    Bvczplus(i) = 0.5d0 * (Bvcz(i) + Bvcz(i+1)) + &
         (dmBvcz(i) - dmBvcz(i+1)) / 6.d0
    Bvczminus(i+1) = Bvczplus(i)
    if(dc_flag.ne.0) then
       psidcplus(i) = 0.5d0 * (psidc(i) + psidc(i+1)) + &
            (dmpsidc(i) - dmpsidc(i+1)) / 6.d0
       psidcminus(i+1) = psidcplus(i)
    endif
  end do
  if (poly .eq. 0) then
    do i = 2, nx-2
      epsplus(i) = 0.5d0 * (eps(i) + eps(i+1)) + &
           (dmeps(i) - dmeps(i+1)) / 6.d0
      epsminus(i+1) = epsplus(i)
    end do
  end if








  if (ppm_detect .ne. 0) then

    do i = 3, nx - 2
      if ( (d2rho(i+1)*d2rho(i-1) < 0.d0).and.(abs(rho(i+1)-rho(i-1)) - &
           ppm_epsilon_shock * min(abs(rho(i+1)), abs(rho(i-1))) > 0.d0) ) then
        etatilde = (rho(i-2) - rho(i+2) + 4.d0 * drho(i)) / (drho(i) * 12.d0)
      else
        etatilde = 0.d0
      end if
      eta = max(0.d0, min(1.d0, ppm_eta1 * (etatilde - ppm_eta2)))
      if (ppm_k0 * abs(drho(i)) * min(press(i-1),press(i+1)) < &
           abs(dpress(i)) * min(rho(i-1), rho(i+1))) then
        eta = 0.d0
      end if
      if (eta > 0.d0) then
        trivial_rp(i-1) = .false.
        trivial_rp(i) = .false.
      end if
      rhominus(i) = rhominus(i) * (1.d0 - eta) + &
           (rho(i-1) + 0.5d0 * dmrho(i-1)) * eta
      rhoplus(i) = rhoplus(i) * (1.d0 - eta) + &
           (rho(i+1) - 0.5d0 * dmrho(i+1)) * eta
    end do

  end if

  



  if (ppm_mppm .gt. 0) then
    l_ev_l=0.d0
    l_ev_r=0.d0
    xwind=0.d0
    do i=3, nx - 3
      agxx = 0.5d0*( gxx(i) + gxx(i+1) )
      agxy = 0.5d0*( gxy(i) + gxy(i+1) )
      agxz = 0.5d0*( gxz(i) + gxz(i+1) )
      agyy = 0.5d0*( gyy(i) + gyy(i+1) )
      agyz = 0.5d0*( gyz(i) + gyz(i+1) )
      agzz = 0.5d0*( gzz(i) + gzz(i+1) )
      det = (-( agxz)**2*(                               agyy) + 2.0d0*( agxy)*( agxz)*( agyz) - (agxx)*( agyz)**2 - ( agxy)**2*( a&
  &gzz)    + (agxx)*(                               agyy)*( agzz))

      call UpperMetric (uxx, uxy, uxz, uyy, uyz, uzz, &
                        det, agxx, agxy, agxz, agyy, agyz, agzz)
      call eigenvaluesM(handle,&
                 (0.5d0 * (rho(i) + rho(i+1))), (0.5d0 * (velx(i) + velx(i+1))), (0.5d0 * (vely(i) + vely(i+1))), (0.5d0 * (velz(i)&
  & + velz(i+1))), &
                 (0.5d0 * (eps(i) + eps(i+1))), (0.5d0 * (press(i) + press(i+1))), (0.5d0 * (w_lorentz(i) + w_lorentz(i+1))),&
                 (0.5d0 * (Bvcx(i) + Bvcx(i+1))), (0.5d0 * (Bvcy(i) + Bvcy(i+1))), (0.5d0 * (Bvcz(i) + Bvcz(i+1))), lam, &
                 agxx, agxy, agxz, agyy, agyz, agzz, &
                 uxx, (0.5d0 * (alp(i) + alp(i+1))), (0.5d0 * (beta(i) + beta(i+1))))
      l_ev_l(i)=lam(1)
      l_ev_r(i)=lam(5)
      xwind(i) = (lam(1) + lam(5)) / (abs(lam(1)) + abs(lam(5)))
      xwind(i) = min(1.d0, max(-1.d0, xwind(i)))

      if (xwind(i) .lt. 0.0d0) then
         rhoplus(i)   =       abs(xwind(i))  * (13.d0*rho(i+1)-5.d0*rho(i+2)+rho(i+3)+3.d0*rho(i  ))/12.d0 +                       &
  &                    (1.d0-abs(xwind(i))) *  rhoplus(i)
         rhominus(i+1)=       abs(xwind(i))  * (13.d0*rho(i+1)-5.d0*rho(i+2)+rho(i+3)+3.d0*rho(i  ))/12.d0 +                       &
  &                    (1.d0-abs(xwind(i))) *  rhominus(i+1)
         velxplus(i)   =       abs(xwind(i))  * (13.d0*velx(i+1)-5.d0*velx(i+2)+velx(i+3)+3.d0*velx(i  ))/12.d0 +                  &
  &                         (1.d0-abs(xwind(i))) *  velxplus(i)
         velxminus(i+1)=       abs(xwind(i))  * (13.d0*velx(i+1)-5.d0*velx(i+2)+velx(i+3)+3.d0*velx(i  ))/12.d0 +                  &
  &                         (1.d0-abs(xwind(i))) *  velxminus(i+1)
         velyplus(i)   =       abs(xwind(i))  * (13.d0*vely(i+1)-5.d0*vely(i+2)+vely(i+3)+3.d0*vely(i  ))/12.d0 +                  &
  &                         (1.d0-abs(xwind(i))) *  velyplus(i)
         velyminus(i+1)=       abs(xwind(i))  * (13.d0*vely(i+1)-5.d0*vely(i+2)+vely(i+3)+3.d0*vely(i  ))/12.d0 +                  &
  &                         (1.d0-abs(xwind(i))) *  velyminus(i+1)
         velzplus(i)   =       abs(xwind(i))  * (13.d0*velz(i+1)-5.d0*velz(i+2)+velz(i+3)+3.d0*velz(i  ))/12.d0 +                  &
  &                         (1.d0-abs(xwind(i))) *  velzplus(i)
         velzminus(i+1)=       abs(xwind(i))  * (13.d0*velz(i+1)-5.d0*velz(i+2)+velz(i+3)+3.d0*velz(i  ))/12.d0 +                  &
  &                         (1.d0-abs(xwind(i))) *  velzminus(i+1)
         Bvcxplus(i)   =       abs(xwind(i))  * (13.d0*Bvcx(i+1)-5.d0*Bvcx(i+2)+Bvcx(i+3)+3.d0*Bvcx(i  ))/12.d0 +                  &
  &                         (1.d0-abs(xwind(i))) *  Bvcxplus(i)
         Bvcxminus(i+1)=       abs(xwind(i))  * (13.d0*Bvcx(i+1)-5.d0*Bvcx(i+2)+Bvcx(i+3)+3.d0*Bvcx(i  ))/12.d0 +                  &
  &                         (1.d0-abs(xwind(i))) *  Bvcxminus(i+1)
         Bvcyplus(i)   =       abs(xwind(i))  * (13.d0*Bvcy(i+1)-5.d0*Bvcy(i+2)+Bvcy(i+3)+3.d0*Bvcy(i  ))/12.d0 +                  &
  &                         (1.d0-abs(xwind(i))) *  Bvcyplus(i)
         Bvcyminus(i+1)=       abs(xwind(i))  * (13.d0*Bvcy(i+1)-5.d0*Bvcy(i+2)+Bvcy(i+3)+3.d0*Bvcy(i  ))/12.d0 +                  &
  &                         (1.d0-abs(xwind(i))) *  Bvcyminus(i+1)
         Bvczplus(i)   =       abs(xwind(i))  * (13.d0*Bvcz(i+1)-5.d0*Bvcz(i+2)+Bvcz(i+3)+3.d0*Bvcz(i  ))/12.d0 +                  &
  &                         (1.d0-abs(xwind(i))) *  Bvczplus(i)
         Bvczminus(i+1)=       abs(xwind(i))  * (13.d0*Bvcz(i+1)-5.d0*Bvcz(i+2)+Bvcz(i+3)+3.d0*Bvcz(i  ))/12.d0 +                  &
  &                         (1.d0-abs(xwind(i))) *  Bvczminus(i+1)
        if(dc_flag.ne.0) then
            psidcplus(i)   =       abs(xwind(i))  * (13.d0*psidc(i+1)-5.d0*psidc(i+2)+psidc(i+3)+3.d0*psidc(i  ))/12.d0 +          &
  &                                 (1.d0-abs(xwind(i))) *  psidcplus(i)
            psidcminus(i+1)=       abs(xwind(i))  * (13.d0*psidc(i+1)-5.d0*psidc(i+2)+psidc(i+3)+3.d0*psidc(i  ))/12.d0 +          &
  &                                 (1.d0-abs(xwind(i))) *  psidcminus(i+1)
        endif
        if (poly .eq. 0) then
           epsplus(i)   =       abs(xwind(i))  * (13.d0*eps(i+1)-5.d0*eps(i+2)+eps(i+3)+3.d0*eps(i  ))/12.d0 +                     &
  &                      (1.d0-abs(xwind(i))) *  epsplus(i)
           epsminus(i+1)=       abs(xwind(i))  * (13.d0*eps(i+1)-5.d0*eps(i+2)+eps(i+3)+3.d0*eps(i  ))/12.d0 +                     &
  &                      (1.d0-abs(xwind(i))) *  epsminus(i+1)
        end if
      else
         rhoplus(i)   =       abs(xwind(i))  * (13.d0*rho(i  )-5.d0*rho(i-1)+rho(i-2)+3.d0*rho(i+1))/12.d0 +                       &
  &                    (1.d0-abs(xwind(i))) *  rhoplus(i)
         rhominus(i+1)=       abs(xwind(i))  * (13.d0*rho(i  )-5.d0*rho(i-1)+rho(i-2)+3.d0*rho(i+1))/12.d0 +                       &
  &                    (1.d0-abs(xwind(i))) *  rhominus(i+1)
         velxplus(i)   =       abs(xwind(i))  * (13.d0*velx(i  )-5.d0*velx(i-1)+velx(i-2)+3.d0*velx(i+1))/12.d0 +                  &
  &                         (1.d0-abs(xwind(i))) *  velxplus(i)
         velxminus(i+1)=       abs(xwind(i))  * (13.d0*velx(i  )-5.d0*velx(i-1)+velx(i-2)+3.d0*velx(i+1))/12.d0 +                  &
  &                         (1.d0-abs(xwind(i))) *  velxminus(i+1)
         velyplus(i)   =       abs(xwind(i))  * (13.d0*vely(i  )-5.d0*vely(i-1)+vely(i-2)+3.d0*vely(i+1))/12.d0 +                  &
  &                         (1.d0-abs(xwind(i))) *  velyplus(i)
         velyminus(i+1)=       abs(xwind(i))  * (13.d0*vely(i  )-5.d0*vely(i-1)+vely(i-2)+3.d0*vely(i+1))/12.d0 +                  &
  &                         (1.d0-abs(xwind(i))) *  velyminus(i+1)
         velzplus(i)   =       abs(xwind(i))  * (13.d0*velz(i  )-5.d0*velz(i-1)+velz(i-2)+3.d0*velz(i+1))/12.d0 +                  &
  &                         (1.d0-abs(xwind(i))) *  velzplus(i)
         velzminus(i+1)=       abs(xwind(i))  * (13.d0*velz(i  )-5.d0*velz(i-1)+velz(i-2)+3.d0*velz(i+1))/12.d0 +                  &
  &                         (1.d0-abs(xwind(i))) *  velzminus(i+1)
         Bvcxplus(i)   =       abs(xwind(i))  * (13.d0*Bvcx(i  )-5.d0*Bvcx(i-1)+Bvcx(i-2)+3.d0*Bvcx(i+1))/12.d0 +                  &
  &                         (1.d0-abs(xwind(i))) *  Bvcxplus(i)
         Bvcxminus(i+1)=       abs(xwind(i))  * (13.d0*Bvcx(i  )-5.d0*Bvcx(i-1)+Bvcx(i-2)+3.d0*Bvcx(i+1))/12.d0 +                  &
  &                         (1.d0-abs(xwind(i))) *  Bvcxminus(i+1)
         Bvcyplus(i)   =       abs(xwind(i))  * (13.d0*Bvcy(i  )-5.d0*Bvcy(i-1)+Bvcy(i-2)+3.d0*Bvcy(i+1))/12.d0 +                  &
  &                         (1.d0-abs(xwind(i))) *  Bvcyplus(i)
         Bvcyminus(i+1)=       abs(xwind(i))  * (13.d0*Bvcy(i  )-5.d0*Bvcy(i-1)+Bvcy(i-2)+3.d0*Bvcy(i+1))/12.d0 +                  &
  &                         (1.d0-abs(xwind(i))) *  Bvcyminus(i+1)
         Bvczplus(i)   =       abs(xwind(i))  * (13.d0*Bvcz(i  )-5.d0*Bvcz(i-1)+Bvcz(i-2)+3.d0*Bvcz(i+1))/12.d0 +                  &
  &                         (1.d0-abs(xwind(i))) *  Bvczplus(i)
         Bvczminus(i+1)=       abs(xwind(i))  * (13.d0*Bvcz(i  )-5.d0*Bvcz(i-1)+Bvcz(i-2)+3.d0*Bvcz(i+1))/12.d0 +                  &
  &                         (1.d0-abs(xwind(i))) *  Bvczminus(i+1)
        if(dc_flag.ne.0) then
            psidcplus(i)   =       abs(xwind(i))  * (13.d0*psidc(i  )-5.d0*psidc(i-1)+psidc(i-2)+3.d0*psidc(i+1))/12.d0 +          &
  &                                 (1.d0-abs(xwind(i))) *  psidcplus(i)
            psidcminus(i+1)=       abs(xwind(i))  * (13.d0*psidc(i  )-5.d0*psidc(i-1)+psidc(i-2)+3.d0*psidc(i+1))/12.d0 +          &
  &                                 (1.d0-abs(xwind(i))) *  psidcminus(i+1)
        endif
        if (poly .eq. 0) then
           epsplus(i)   =       abs(xwind(i))  * (13.d0*eps(i  )-5.d0*eps(i-1)+eps(i-2)+3.d0*eps(i+1))/12.d0 +                     &
  &                      (1.d0-abs(xwind(i))) *  epsplus(i)
           epsminus(i+1)=       abs(xwind(i))  * (13.d0*eps(i  )-5.d0*eps(i-1)+eps(i-2)+3.d0*eps(i+1))/12.d0 +                     &
  &                      (1.d0-abs(xwind(i))) *  epsminus(i+1)
        end if
      end if
      if (rho(i+1) .gt. rho(i)) then
  rhoplus(i)=min(rho(i+1),max(rho(i), rhoplus(i)))
  else
  rhoplus(i)=min(rho(i),max(rho(i+1), rhoplus(i)))
  endif
      if (rho(i+1) .gt. rho(i)) then
  rhominus(i+1)=min(rho(i+1),max(rho(i), rhominus(i+1)))
  else
  rhominus(i+1)=min(rho(i),max(rho(i+1), rhominus(i+1)))
  endif
      if (velx(i+1) .gt. velx(i)) then
  velxplus(i)=min(velx(i+1),max(velx(i), velxplus(i)))
  else
  velxplus(i)=min(velx(i),max(velx(i+1), velxplus(i)))
  endif
      if (velx(i+1) .gt. velx(i)) then
  velxminus(i+1)=min(velx(i+1),max(velx(i), velxminus(i+1)))
  else
  velxminus(i+1)=min(velx(i),max(velx(i+1), velxminus(i+1)))
  endif
      if (vely(i+1) .gt. vely(i)) then
  velyplus(i)=min(vely(i+1),max(vely(i), velyplus(i)))
  else
  velyplus(i)=min(vely(i),max(vely(i+1), velyplus(i)))
  endif
      if (vely(i+1) .gt. vely(i)) then
  velyminus(i+1)=min(vely(i+1),max(vely(i), velyminus(i+1)))
  else
  velyminus(i+1)=min(vely(i),max(vely(i+1), velyminus(i+1)))
  endif
      if (velz(i+1) .gt. velz(i)) then
  velzplus(i)=min(velz(i+1),max(velz(i), velzplus(i)))
  else
  velzplus(i)=min(velz(i),max(velz(i+1), velzplus(i)))
  endif
      if (velz(i+1) .gt. velz(i)) then
  velzminus(i+1)=min(velz(i+1),max(velz(i), velzminus(i+1)))
  else
  velzminus(i+1)=min(velz(i),max(velz(i+1), velzminus(i+1)))
  endif
      if (Bvcx(i+1) .gt. Bvcx(i)) then
  Bvcxplus(i)=min(Bvcx(i+1),max(Bvcx(i), Bvcxplus(i)))
  else
  Bvcxplus(i)=min(Bvcx(i),max(Bvcx(i+1), Bvcxplus(i)))
  endif
      if (Bvcx(i+1) .gt. Bvcx(i)) then
  Bvcxminus(i+1)=min(Bvcx(i+1),max(Bvcx(i), Bvcxminus(i+1)))
  else
  Bvcxminus(i+1)=min(Bvcx(i),max(Bvcx(i+1), Bvcxminus(i+1)))
  endif
      if (Bvcy(i+1) .gt. Bvcy(i)) then
  Bvcyplus(i)=min(Bvcy(i+1),max(Bvcy(i), Bvcyplus(i)))
  else
  Bvcyplus(i)=min(Bvcy(i),max(Bvcy(i+1), Bvcyplus(i)))
  endif
      if (Bvcy(i+1) .gt. Bvcy(i)) then
  Bvcyminus(i+1)=min(Bvcy(i+1),max(Bvcy(i), Bvcyminus(i+1)))
  else
  Bvcyminus(i+1)=min(Bvcy(i),max(Bvcy(i+1), Bvcyminus(i+1)))
  endif
      if (Bvcz(i+1) .gt. Bvcz(i)) then
  Bvczplus(i)=min(Bvcz(i+1),max(Bvcz(i), Bvczplus(i)))
  else
  Bvczplus(i)=min(Bvcz(i),max(Bvcz(i+1), Bvczplus(i)))
  endif
      if (Bvcz(i+1) .gt. Bvcz(i)) then
  Bvczminus(i+1)=min(Bvcz(i+1),max(Bvcz(i), Bvczminus(i+1)))
  else
  Bvczminus(i+1)=min(Bvcz(i),max(Bvcz(i+1), Bvczminus(i+1)))
  endif
      if(dc_flag.ne.0) then
         if (psidc(i+1) .gt. psidc(i)) then
  psidcplus(i)=min(psidc(i+1),max(psidc(i), psidcplus(i)))
  else
  psidcplus(i)=min(psidc(i),max(psidc(i+1), psidcplus(i)))
  endif
         if (psidc(i+1) .gt. psidc(i)) then
  psidcminus(i+1)=min(psidc(i+1),max(psidc(i), psidcminus(i+1)))
  else
  psidcminus(i+1)=min(psidc(i),max(psidc(i+1), psidcminus(i+1)))
  endif
      endif
      if (poly .eq. 0) then
        if (eps(i+1) .gt. eps(i)) then
  epsplus(i)=min(eps(i+1),max(eps(i), epsplus(i)))
  else
  epsplus(i)=min(eps(i),max(eps(i+1), epsplus(i)))
  endif
        if (eps(i+1) .gt. eps(i)) then
  epsminus(i+1)=min(eps(i+1),max(eps(i), epsminus(i+1)))
  else
  epsminus(i+1)=min(eps(i),max(eps(i+1), epsminus(i+1)))
  endif
      end if



    end do
    
    if (ppm_mppm_debug_eigenvalues .gt. 0) then
      if (dir .eq. 1) then
        ev_l(:,ni,nj) = l_ev_l
        ev_r(:,ni,nj) = l_ev_r
        xw(:,ni,nj) = xwind
      else if (dir .eq. 2) then
        ev_l(ni,:,nj) = l_ev_l
        ev_r(ni,:,nj) = l_ev_r
        xw(ni,:,nj) = xwind
      else if (dir .eq. 3) then
        ev_l(ni,nj,:) = l_ev_l
        ev_r(ni,nj,:) = l_ev_r
        xw(ni,nj,:) = xwind
      else
        write (*,*) "flux direction not 1 to 3 ?"
      end if
    end if
  end if



  do i = 3, nx - 2
    dpress2 = press(i+2) - press(i-2)
    dvel = vx(i+1) - vx(i-1)
    if ( (abs(dpress(i)) >  ppm_epsilon * min(press(i-1),press(i+1))) .and. &
         (dvel < 0.d0) ) then
      w = 1.d0
    else
      w = 0.d0
    end if
    if (abs(dpress2) < ppm_small) then
      tilde_flatten(i) = 1.d0
    else
      tilde_flatten(i) = max(0.d0, 1.d0 - w * max(0.d0, ppm_omega2 * &
           (dpress(i) / dpress2 - ppm_omega1)))
    end if
  end do



  if (PPM3) then 
    do i = 3, nx - 2
      flatten = tilde_flatten(i)
      if (abs(1.d0 - flatten) > 0.d0) then
        trivial_rp(i-1) = .false.
        trivial_rp(i) = .false.
      end if
      rhoplus(i) = flatten * rhoplus(i) + (1.d0 - flatten) * rho(i)
      rhominus(i) = flatten * rhominus(i) + (1.d0 - flatten) * rho(i)
      velxplus(i) = flatten * velxplus(i) + (1.d0 - flatten) * vx(i)
      velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * vx(i)
      velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vy(i)
      velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vy(i)
      velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * vz(i)
      velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * vz(i)
      Bvcxplus(i) = flatten * Bvcxplus(i) + (1.d0 - flatten) * Bvcx(i)
      Bvcxminus(i) = flatten * Bvcxminus(i) + (1.d0 - flatten) * Bvcx(i)
      Bvcyplus(i) = flatten * Bvcyplus(i) + (1.d0 - flatten) * Bvcy(i)
      Bvcyminus(i) = flatten * Bvcyminus(i) + (1.d0 - flatten) * Bvcy(i)
      Bvczplus(i) = flatten * Bvczplus(i) + (1.d0 - flatten) * Bvcz(i)
      Bvczminus(i) = flatten * Bvczminus(i) + (1.d0 - flatten) * Bvcz(i)
      if(dc_flag.ne.0) then
         psidcplus(i) = flatten * psidcplus(i) + (1.d0 - flatten) * psidc(i)
         psidcminus(i) = flatten * psidcminus(i) + (1.d0 - flatten) * psidc(i)
      endif
      if (poly .eq. 0) then
        epsplus(i) = flatten * epsplus(i) + (1.d0 - flatten) * eps(i)
        epsminus(i) = flatten * epsminus(i) + (1.d0 - flatten) * eps(i)
      end if
    end do
  else  
    do i = 4, nx - 3
      s=int(sign(1.d0, -dpress(i)))
      flatten = max(tilde_flatten(i), tilde_flatten(i+s))  
      if (abs(1.d0 - flatten) > 0.d0) then
        trivial_rp(i-1) = .false.
        trivial_rp(i) = .false.
      end if
      rhoplus(i) = flatten * rhoplus(i) + (1.d0 - flatten) * rho(i)
      rhominus(i) = flatten * rhominus(i) + (1.d0 - flatten) * rho(i)
      velxplus(i) = flatten * velxplus(i) + (1.d0 - flatten) * vx(i)
      velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * vx(i)
      velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vy(i)
      velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vy(i)
      velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * vz(i)
      velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * vz(i)
      Bvcxplus(i) = flatten * Bvcxplus(i) + (1.d0 - flatten) * Bvcx(i)
      Bvcxminus(i) = flatten * Bvcxminus(i) + (1.d0 - flatten) * Bvcx(i)
      Bvcyplus(i) = flatten * Bvcyplus(i) + (1.d0 - flatten) * Bvcy(i)
      Bvcyminus(i) = flatten * Bvcyminus(i) + (1.d0 - flatten) * Bvcy(i)
      Bvczplus(i) = flatten * Bvczplus(i) + (1.d0 - flatten) * Bvcz(i)
      Bvczminus(i) = flatten * Bvczminus(i) + (1.d0 - flatten) * Bvcz(i)
      if(dc_flag.ne.0) then
         psidcplus(i) = flatten * psidcplus(i) + (1.d0 - flatten) * psidc(i)
         psidcminus(i) = flatten * psidcminus(i) + (1.d0 - flatten) * psidc(i)
      endif
      if (poly .eq. 0) then
        epsplus(i) = flatten * epsplus(i) + (1.d0 - flatten) * eps(i)
        epsminus(i) = flatten * epsminus(i) + (1.d0 - flatten) * eps(i)
      end if
    end do
  end if




do i = GRHydro_stencil, nx - GRHydro_stencil + 1

    if (.not.( (rhoplus(i).eq.rho(i)) .and. (rho(i).eq.rhominus(i)) )             .and. ((rhoplus(i)-rho(i))*(rho(i)-rhominus(i)) .&
  &le. 0.d0)) then
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  rhominus(i) = rho(i)
  rhoplus(i) = rho(i)
  else if (6.d0 * (rhoplus(i) - rhominus(i)) * (rho(i) - 0.5d0 *                      (rhoplus(i) + rhominus(i))) >                &
  &                         (rhoplus(i) - rhominus(i))**2) then
  rhominus(i) = 3.d0 * rho(i) - 2.d0 * rhoplus(i)
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  else if (6.d0 * (rhoplus(i) - rhominus(i)) * (rho(i) - 0.5d0 *                      (rhoplus(i) + rhominus(i))) <                &
  &                        -(rhoplus(i) - rhominus(i))**2) then
  rhoplus(i) = 3.d0 * rho(i) - 2.d0 * rhominus(i)
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  end if
  if (.not.( (rhoplus(i).eq.rho(i)) .and. (rho(i).eq.rhominus(i)) ) ) then
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  end if
    if (.not.( (velxplus(i).eq.vx(i)) .and. (vx(i).eq.velxminus(i)) )             .and. ((velxplus(i)-vx(i))*(vx(i)-velxminus(i)) .&
  &le. 0.d0)) then
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  velxminus(i) = vx(i)
  velxplus(i) = vx(i)
  else if (6.d0 * (velxplus(i) - velxminus(i)) * (vx(i) - 0.5d0 *                      (velxplus(i) + velxminus(i))) >             &
  &                            (velxplus(i) - velxminus(i))**2) then
  velxminus(i) = 3.d0 * vx(i) - 2.d0 * velxplus(i)
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  else if (6.d0 * (velxplus(i) - velxminus(i)) * (vx(i) - 0.5d0 *                      (velxplus(i) + velxminus(i))) <             &
  &                           -(velxplus(i) - velxminus(i))**2) then
  velxplus(i) = 3.d0 * vx(i) - 2.d0 * velxminus(i)
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  end if
  if (.not.( (velxplus(i).eq.vx(i)) .and. (vx(i).eq.velxminus(i)) ) ) then
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  end if
    if (.not.( (velyplus(i).eq.vy(i)) .and. (vy(i).eq.velyminus(i)) )             .and. ((velyplus(i)-vy(i))*(vy(i)-velyminus(i)) .&
  &le. 0.d0)) then
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  velyminus(i) = vy(i)
  velyplus(i) = vy(i)
  else if (6.d0 * (velyplus(i) - velyminus(i)) * (vy(i) - 0.5d0 *                      (velyplus(i) + velyminus(i))) >             &
  &                            (velyplus(i) - velyminus(i))**2) then
  velyminus(i) = 3.d0 * vy(i) - 2.d0 * velyplus(i)
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  else if (6.d0 * (velyplus(i) - velyminus(i)) * (vy(i) - 0.5d0 *                      (velyplus(i) + velyminus(i))) <             &
  &                           -(velyplus(i) - velyminus(i))**2) then
  velyplus(i) = 3.d0 * vy(i) - 2.d0 * velyminus(i)
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  end if
  if (.not.( (velyplus(i).eq.vy(i)) .and. (vy(i).eq.velyminus(i)) ) ) then
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  end if
    if (.not.( (velzplus(i).eq.vz(i)) .and. (vz(i).eq.velzminus(i)) )             .and. ((velzplus(i)-vz(i))*(vz(i)-velzminus(i)) .&
  &le. 0.d0)) then
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  velzminus(i) = vz(i)
  velzplus(i) = vz(i)
  else if (6.d0 * (velzplus(i) - velzminus(i)) * (vz(i) - 0.5d0 *                      (velzplus(i) + velzminus(i))) >             &
  &                            (velzplus(i) - velzminus(i))**2) then
  velzminus(i) = 3.d0 * vz(i) - 2.d0 * velzplus(i)
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  else if (6.d0 * (velzplus(i) - velzminus(i)) * (vz(i) - 0.5d0 *                      (velzplus(i) + velzminus(i))) <             &
  &                           -(velzplus(i) - velzminus(i))**2) then
  velzplus(i) = 3.d0 * vz(i) - 2.d0 * velzminus(i)
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  end if
  if (.not.( (velzplus(i).eq.vz(i)) .and. (vz(i).eq.velzminus(i)) ) ) then
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  end if
    if (.not.( (Bvcxplus(i).eq.Bvcx(i)) .and. (Bvcx(i).eq.Bvcxminus(i)) )             .and. ((Bvcxplus(i)-Bvcx(i))*(Bvcx(i)-Bvcxmin&
  &us(i)) .le. 0.d0)) then
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  Bvcxminus(i) = Bvcx(i)
  Bvcxplus(i) = Bvcx(i)
  else if (6.d0 * (Bvcxplus(i) - Bvcxminus(i)) * (Bvcx(i) - 0.5d0 *                      (Bvcxplus(i) + Bvcxminus(i))) >           &
  &                              (Bvcxplus(i) - Bvcxminus(i))**2) then
  Bvcxminus(i) = 3.d0 * Bvcx(i) - 2.d0 * Bvcxplus(i)
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  else if (6.d0 * (Bvcxplus(i) - Bvcxminus(i)) * (Bvcx(i) - 0.5d0 *                      (Bvcxplus(i) + Bvcxminus(i))) <           &
  &                             -(Bvcxplus(i) - Bvcxminus(i))**2) then
  Bvcxplus(i) = 3.d0 * Bvcx(i) - 2.d0 * Bvcxminus(i)
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  end if
  if (.not.( (Bvcxplus(i).eq.Bvcx(i)) .and. (Bvcx(i).eq.Bvcxminus(i)) ) ) then
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  end if
    if (.not.( (Bvcyplus(i).eq.Bvcy(i)) .and. (Bvcy(i).eq.Bvcyminus(i)) )             .and. ((Bvcyplus(i)-Bvcy(i))*(Bvcy(i)-Bvcymin&
  &us(i)) .le. 0.d0)) then
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  Bvcyminus(i) = Bvcy(i)
  Bvcyplus(i) = Bvcy(i)
  else if (6.d0 * (Bvcyplus(i) - Bvcyminus(i)) * (Bvcy(i) - 0.5d0 *                      (Bvcyplus(i) + Bvcyminus(i))) >           &
  &                              (Bvcyplus(i) - Bvcyminus(i))**2) then
  Bvcyminus(i) = 3.d0 * Bvcy(i) - 2.d0 * Bvcyplus(i)
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  else if (6.d0 * (Bvcyplus(i) - Bvcyminus(i)) * (Bvcy(i) - 0.5d0 *                      (Bvcyplus(i) + Bvcyminus(i))) <           &
  &                             -(Bvcyplus(i) - Bvcyminus(i))**2) then
  Bvcyplus(i) = 3.d0 * Bvcy(i) - 2.d0 * Bvcyminus(i)
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  end if
  if (.not.( (Bvcyplus(i).eq.Bvcy(i)) .and. (Bvcy(i).eq.Bvcyminus(i)) ) ) then
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  end if
    if (.not.( (Bvczplus(i).eq.Bvcz(i)) .and. (Bvcz(i).eq.Bvczminus(i)) )             .and. ((Bvczplus(i)-Bvcz(i))*(Bvcz(i)-Bvczmin&
  &us(i)) .le. 0.d0)) then
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  Bvczminus(i) = Bvcz(i)
  Bvczplus(i) = Bvcz(i)
  else if (6.d0 * (Bvczplus(i) - Bvczminus(i)) * (Bvcz(i) - 0.5d0 *                      (Bvczplus(i) + Bvczminus(i))) >           &
  &                              (Bvczplus(i) - Bvczminus(i))**2) then
  Bvczminus(i) = 3.d0 * Bvcz(i) - 2.d0 * Bvczplus(i)
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  else if (6.d0 * (Bvczplus(i) - Bvczminus(i)) * (Bvcz(i) - 0.5d0 *                      (Bvczplus(i) + Bvczminus(i))) <           &
  &                             -(Bvczplus(i) - Bvczminus(i))**2) then
  Bvczplus(i) = 3.d0 * Bvcz(i) - 2.d0 * Bvczminus(i)
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  end if
  if (.not.( (Bvczplus(i).eq.Bvcz(i)) .and. (Bvcz(i).eq.Bvczminus(i)) ) ) then
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  end if
    if(dc_flag.ne.0) then
       if (.not.( (psidcplus(i).eq.psidc(i)) .and. (psidc(i).eq.psidcminus(i)) )             .and. ((psidcplus(i)-psidc(i))*(psidc(&
  &i)-psidcminus(i)) .le. 0.d0)) then
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  psidcminus(i) = psidc(i)
  psidcplus(i) = psidc(i)
  else if (6.d0 * (psidcplus(i) - psidcminus(i)) * (psidc(i) - 0.5d0 *                      (psidcplus(i) + psidcminus(i))) >      &
  &                                   (psidcplus(i) - psidcminus(i))**2) then
  psidcminus(i) = 3.d0 * psidc(i) - 2.d0 * psidcplus(i)
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  else if (6.d0 * (psidcplus(i) - psidcminus(i)) * (psidc(i) - 0.5d0 *                      (psidcplus(i) + psidcminus(i))) <      &
  &                                  -(psidcplus(i) - psidcminus(i))**2) then
  psidcplus(i) = 3.d0 * psidc(i) - 2.d0 * psidcminus(i)
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  end if
  if (.not.( (psidcplus(i).eq.psidc(i)) .and. (psidc(i).eq.psidcminus(i)) ) ) then
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  end if
    endif
  end do
  if (poly .eq. 0) then
    do i = GRHydro_stencil, nx - GRHydro_stencil + 1
      if (.not.( (epsplus(i).eq.eps(i)) .and. (eps(i).eq.epsminus(i)) )             .and. ((epsplus(i)-eps(i))*(eps(i)-epsminus(i))&
  & .le. 0.d0)) then
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  epsminus(i) = eps(i)
  epsplus(i) = eps(i)
  else if (6.d0 * (epsplus(i) - epsminus(i)) * (eps(i) - 0.5d0 *                      (epsplus(i) + epsminus(i))) >                &
  &                         (epsplus(i) - epsminus(i))**2) then
  epsminus(i) = 3.d0 * eps(i) - 2.d0 * epsplus(i)
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  else if (6.d0 * (epsplus(i) - epsminus(i)) * (eps(i) - 0.5d0 *                      (epsplus(i) + epsminus(i))) <                &
  &                        -(epsplus(i) - epsminus(i))**2) then
  epsplus(i) = 3.d0 * eps(i) - 2.d0 * epsminus(i)
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  end if
  if (.not.( (epsplus(i).eq.eps(i)) .and. (eps(i).eq.epsminus(i)) ) ) then
  trivial_rp(i-1) = .false.
  trivial_rp(i) = .false.
  end if
    end do
  end if

  if (check_for_trivial_rp .eq. 0) then
    trivial_rp = .false.
  end if

  
  do i = 1, nx
    if (GRHydro_enable_internal_excision /= 0 .and. &
        (hydro_excision_mask(i) .ne. 0)) then
      if (i .gt. 1) then
        trivial_rp(i-1)=.true.
      end if
      trivial_rp(i)=.true.
    else
      
      
      
      cond = .false.
      if (i .gt. 1 .and. GRHydro_enable_internal_excision /= 0) then
        cond = hydro_excision_mask(i-1) .ne. 0
      end if
      if (cond) then
        rhominus(i)=rho(i)
        rhoplus(i)=rho(i)
        velxminus(i)=vx(i)
        velxplus(i)=vx(i)
        velyminus(i)=vy(i)
        velyplus(i)=vy(i)
        velzminus(i)=vz(i)
        velzplus(i)=vz(i)
        Bvcxminus(i)=Bvcx(i)  
        Bvcxplus(i)=Bvcx(i)
        Bvcyminus(i)=Bvcy(i)
        Bvcyplus(i)=Bvcy(i)
        Bvczminus(i)=Bvcz(i)
        Bvczplus(i)=Bvcz(i)
        rhominus(i-1)=rho(i)
        rhoplus(i-1)=rho(i)
        velxminus(i-1)=vx(i)
        velxplus(i-1)=vx(i)
        velyminus(i-1)=vy(i)
        velyplus(i-1)=vy(i)
        velzminus(i-1)=vz(i)
        velzplus(i-1)=vz(i)
        Bvcxminus(i-1)=Bvcx(i)
        Bvcxplus(i-1)=Bvcx(i)
        Bvcyminus(i-1)=Bvcy(i)
        Bvcyplus(i-1)=Bvcy(i)
        Bvczminus(i-1)=Bvcz(i)
        Bvczplus(i-1)=Bvcz(i)
        if (dc_flag.ne.0) then
          psidcminus(i)=psidc(i)
          psidcplus(i)=psidc(i)
          psidcminus(i-1)=psidc(i)
          psidcplus(i-1)=psidc(i)
        end if
       if (poly .eq. 0) then
          epsminus(i)=eps(i)
          epsplus(i)=eps(i)
          epsminus(i-1)=eps(i)
          epsplus(i-1)=eps(i)
        end if
      else
        cond = .false.
        if ((i.gt.2) .and. (i.lt.nx) .and. GRHydro_enable_internal_excision /= 0) then
          cond = (ppm_mppm .eq. 0) .and. (hydro_excision_mask(i-2) .ne. 0)
        end if
        if (cond) then
          call PPM_TVD(rho(i-1), rho(i), rho(i+1), rhominus(i), rhoplus(i))
          call PPM_TVD(vx(i-1), vx(i), vx(i+1), velxminus(i), velxplus(i))
          call PPM_TVD(vy(i-1), vy(i), vy(i+1), velyminus(i), velyplus(i))
          call PPM_TVD(vz(i-1), vz(i), vz(i+1), velzminus(i), velzplus(i))
          call PPM_TVD(Bvcx(i-1), Bvcx(i), Bvcx(i+1), Bvcxminus(i), Bvcxplus(i))
          call PPM_TVD(Bvcy(i-1), Bvcy(i), Bvcy(i+1), Bvcyminus(i), Bvcyplus(i))
          call PPM_TVD(Bvcz(i-1), Bvcz(i), Bvcz(i+1), Bvczminus(i), Bvczplus(i))
          if(dc_flag.ne.0) then
             call PPM_TVD(psidc(i-1), psidc(i), psidc(i+1), psidcminus(i), psidcplus(i))
          endif
          if (poly .eq. 0) then
            call PPM_TVD(eps(i-1), eps(i), eps(i+1), epsminus(i), epsplus(i))
          end if
        end if
      end if
      cond = .false.
      if (i .lt. nx .and. GRHydro_enable_internal_excision /= 0) then
        cond = hydro_excision_mask(i+1) .ne. 0
      end if
      if (cond) then
        rhominus(i)=rho(i)
        rhoplus(i)=rho(i)
        velxminus(i)=vx(i)
        velxplus(i)=vx(i)
        velyminus(i)=vy(i)
        velyplus(i)=vy(i)
        velzminus(i)=vz(i)
        velzplus(i)=vz(i)
        Bvcxminus(i)=Bvcx(i)
        Bvcxplus(i)=Bvcx(i)
        Bvcyminus(i)=Bvcy(i)
        Bvcyplus(i)=Bvcy(i)
        Bvczminus(i)=Bvcz(i)
        Bvczplus(i)=Bvcz(i)
        rhominus(i+1)=rho(i)
        rhoplus(i+1)=rho(i)
        velxminus(i+1)=vx(i)
        velxplus(i+1)=vx(i)
        velyminus(i+1)=vy(i)
        velyplus(i+1)=vy(i)
        velzminus(i+1)=vz(i)
        velzplus(i+1)=vz(i)
        Bvcxminus(i+1)=Bvcx(i)
        Bvcxplus(i+1)=Bvcx(i)
        Bvcyminus(i+1)=Bvcy(i)
        Bvcyplus(i+1)=Bvcy(i)
        Bvczminus(i+1)=Bvcz(i)
        Bvczplus(i+1)=Bvcz(i)
        if (dc_flag.ne.0) then
          psidcminus(i)=psidc(i)
          psidcplus(i)=psidc(i)
          psidcminus(i+1)=psidc(i)
          psidcplus(i+1)=psidc(i)
        endif
        if (poly .eq. 0) then
          epsminus(i)=eps(i)
          epsplus(i)=eps(i)
          epsminus(i+1)=eps(i)
          epsplus(i+1)=eps(i)
        endif
      else
        cond = .false.
        if ((i.lt.nx-1) .and. (i.gt.1) .and. GRHydro_enable_internal_excision /= 0) then
          cond = (ppm_mppm .eq. 0) .and. (hydro_excision_mask(i+2) .ne. 0)
        end if
        if (cond) then
          call PPM_TVD(rho(i-1), rho(i), rho(i+1), rhominus(i), rhoplus(i))
          call PPM_TVD(vx(i-1), vx(i), vx(i+1), velxminus(i), velxplus(i))
          call PPM_TVD(vy(i-1), vy(i), vy(i+1), velyminus(i), velyplus(i))
          call PPM_TVD(vz(i-1), vz(i), vz(i+1), velzminus(i), velzplus(i))
          call PPM_TVD(Bvcx(i-1), Bvcx(i), Bvcx(i+1), Bvcxminus(i), Bvcxplus(i))
          call PPM_TVD(Bvcy(i-1), Bvcy(i), Bvcy(i+1), Bvcyminus(i), Bvcyplus(i))
          call PPM_TVD(Bvcz(i-1), Bvcz(i), Bvcz(i+1), Bvczminus(i), Bvczplus(i))
          if(dc_flag.ne.0) then
             call PPM_TVD(Bvcz(i-1), Bvcz(i), Bvcz(i+1), Bvczminus(i), Bvczplus(i))
          endif
          if (poly .eq. 0) then
            call PPM_TVD(eps(i-1), eps(i), eps(i+1), epsminus(i), epsplus(i))
          end if
        end if
      end if
    end if
  end do

  
  if (reconstruct_Wv.ne.0) then
    do i = grhydro_stencil, nx - grhydro_stencil + 1
      
      
      
      
      
      agxx = 0.5d0*( gxx(i) + gxx(i-1) )
      agxy = 0.5d0*( gxy(i) + gxy(i-1) )
      agxz = 0.5d0*( gxz(i) + gxz(i-1) )
      agyy = 0.5d0*( gyy(i) + gyy(i-1) )
      agyz = 0.5d0*( gyz(i) + gyz(i-1) )
      agzz = 0.5d0*( gzz(i) + gzz(i-1) )
      w = sqrt( 1.d0 + agxx*velxminus(i)*velxminus(i) + &
                       agyy*velyminus(i)*velyminus(i) &
              + agzz*velzminus(i)*velzminus(i) + &
                2.d0*agxy*velxminus(i)*velyminus(i) &
              + 2.d0*agxz*velxminus(i)*velzminus(i) + &
                2.d0*agyz*velyminus(i)*velzminus(i) )
      velxminus(i) = velxminus(i)/w
      velyminus(i) = velyminus(i)/w
      velzminus(i) = velzminus(i)/w
      
      agxx = 0.5d0*( gxx(i) + gxx(i+1) )
      agxy = 0.5d0*( gxy(i) + gxy(i+1) )
      agxz = 0.5d0*( gxz(i) + gxz(i+1) )
      agyy = 0.5d0*( gyy(i) + gyy(i+1) )
      agyz = 0.5d0*( gyz(i) + gyz(i+1) )
      agzz = 0.5d0*( gzz(i) + gzz(i+1) )
      w = sqrt( 1.d0 + agxx*velxplus(i)*velxplus(i) + &
                agyy*velyplus(i)*velyplus(i) &
              + agzz*velzplus(i)*velzplus(i) + &
               2.d0*agxy*velxplus(i)*velyplus(i) &
              + 2.d0*agxz*velxplus(i)*velzplus(i) + &
               2.d0*agyz*velyplus(i)*velzplus(i) )
      velxplus(i) = velxplus(i)/w
      velyplus(i) = velyplus(i)/w
      velzplus(i) = velzplus(i)/w
    end do
  end if
  
  
return

end subroutine SimplePPM_1dM


