solve Subroutine

   pure 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_exact), 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) pattern solution.
   real(R8P)                                               :: dum, alfa, beta !< Dummies coefficients.
   real(R8P)                                               :: p_2, p_3        !< Pessure of state 2 and 3.
   real(R8P)                                               :: dp2, dp3        !< Derivate of pessure (dp/du) of state 2 and 3.

   call pattern%initialize(eos_left=eos_left, state_left=state_left, eos_right=eos_right, state_right=state_right, normal=normal)

   ! initiale u23 speed
   if (pattern%p_1 < pattern%p_4) then
     dum  = 0.5_R8P * pattern%eos_4%gm1() / pattern%eos_4%g() ! (gamma - 1) / (gamma * 2)
   else
     dum  = 0.5_R8P * pattern%eos_1%gm1() / pattern%eos_1%g() ! (gamma - 1) / (gamma * 2)
   endif
   alfa = (pattern%p_1 / pattern%p_4) ** dum
   beta = alfa * pattern%eos_1%delta() / pattern%a_1 + pattern%eos_4%delta()/ pattern%a_4
   pattern%u23 = (alfa - 1.0_R8P) / beta +               &
                 0.5_R8P * (pattern%u_1 + pattern%u_4) + &
                 0.5_R8P * (pattern%u_1 - pattern%u_4) * &
                 (alfa * pattern%eos_1%delta() / pattern%a_1 - pattern%eos_4%delta()/ pattern%a_4) / beta

   Newton: do
      call pattern%compute_states23_from_u23(p_2=p_2, p_3=p_3)
      ! evaluate the Newton-Rapson convergence
      if (abs(1.0_R8P - (p_2 / p_3)) >= self%tolerance) then
         dp2 = -1._R8P * pattern%eos_1%g() * p_2 / pattern%a_2
         dp3 =  1._R8P * pattern%eos_4%g() * p_3 / pattern%a_3
         pattern%u23 = pattern%u23 - ((p_2 - p_3) / (dp2-dp3))
      else
        pattern%p23 = p_2 ! p_2 ~= p_3
        exit Newton
      endif
   enddo Newton

   call pattern%compute_fluxes(normal=normal, fluxes=fluxes)
   endsubroutine solve