solve Subroutine

   subroutine solve(self, eos_left, state_left, eos_right, state_right, normal, fluxes)
   !< Solve Riemann Problem.
   !<
   !< Approximate Riemann Solver based on (local) Lax-Friedrichs (known also as Rusanov) algorithm.
   class(riemann_solver_compressible_llf), intent(in)    :: self         !< Solver.
   class(eos_object),                      intent(in)    :: eos_left     !< Equation of state for left state.
   class(conservative_object),             intent(in)    :: state_left   !< Left Riemann state.
   class(eos_object),                      intent(in)    :: eos_right    !< Equation of state for right state.
   class(conservative_object),             intent(in)    :: state_right  !< Right Riemann state.
   type(vector),                           intent(in)    :: normal       !< Normal (versor) of face where fluxes are given.
   class(conservative_object),             intent(inout) :: fluxes       !< Fluxes of the Riemann Problem solution.
   type(riemann_pattern_compressible_pvl)                :: pattern      !< Riemann (states) PVL pattern solution.
   type(conservative_compressible), pointer              :: state_left_  !< Left Riemann state, local variable.
   type(conservative_compressible), pointer              :: state_right_ !< Right Riemann state, local variable.
   type(conservative_compressible)                       :: fluxes_left  !< Fluxes of left state.
   type(conservative_compressible)                       :: fluxes_right !< Fluxes of right state.
   real(R8P)                                             :: lmax         !< Maximum wave speed estimation.

   call pattern%initialize(eos_left=eos_left, state_left=state_left, eos_right=eos_right, state_right=state_right, normal=normal)
   call pattern%compute_waves
   lmax = max(abs(pattern%s_1), abs(pattern%s_4))
   call state_left%compute_fluxes(eos=eos_left, normal=normal, fluxes=fluxes_left)
   call state_right%compute_fluxes(eos=eos_right, normal=normal, fluxes=fluxes_right)
   state_left_ => conservative_compressible_pointer(to=state_left)
   state_right_ => conservative_compressible_pointer(to=state_right)
   select type(fluxes)
   type is(conservative_compressible)
#ifdef __GFORTRAN__
      fluxes = 0.5_R8P * (fluxes_left + fluxes_right - (lmax * (state_right_ - state_left_)))
#else
      error stop 'error: Intel fortran still does not support abstract math!'
#endif
   endselect
   endsubroutine solve