compute_fluxes Subroutine

   elemental subroutine compute_fluxes(self, normal, fluxes)
   !< Compute fluxes at initial discontinuity interface.
   class(riemann_pattern_compressible_object), intent(in)    :: self   !< Riemann (states) pattern solution.
   type(vector),                               intent(in)    :: normal !< Normal (versor) of face where fluxes are given.
   class(conservative_object),                 intent(inout) :: fluxes !< Fluxes at initial discontinuity interface.
   real(R8P)                                                 :: a      !< Speed of sound at interface.
   real(R8P)                                                 :: p      !< Pressure at interface.
   real(R8P)                                                 :: r      !< Desnity at interface.

   ! call fluxes%destroy
   associate(s1=>self%s_1, s2=>self%s_2, u23=>self%u23, s3=>self%s_3, s4=>self%s_4, &
             eos1=>self%eos_1, eos4=>self%eos_4,                                    &
             p1=>self%p_1, r1=>self%r_1, u1=>self%u_1, a1=>self%a_1,                &
             p4=>self%p_4, r4=>self%r_4, u4=>self%u_4, a4=>self%a_4,                &
             p23=>self%p23, r2=>self%r_2, r3=>self%r_3)
      select type(fluxes)
      class is(conservative_compressible)
         select case(minloc([-s1, s1 * s2, s2 * u23, u23 * s3, s3 * s4, s4], dim=1))
         case(1) ! left supersonic
            call fluxes%compute_fluxes_from_primitive(eos=eos1, p=p1, r=r1, u=u1, normal=normal)
         case(2) ! left transonic
            a = (a1 + u1 * eos1%delta()) / (1._R8P + eos1%delta())
            p = p1* (a / a1) ** eos1%eta()
            r = eos1%density(pressure=p, speed_of_sound=a)
            call fluxes%compute_fluxes_from_primitive(eos=eos1, p=p, r=r, u=a, normal=normal)
         case(3) ! left subsonic
            call fluxes%compute_fluxes_from_primitive(eos=eos1, p=p23, r=r2, u=u23, normal=normal)
         case(4) ! right subsonic
            call fluxes%compute_fluxes_from_primitive(eos=eos4, p=p23, r=r3, u=u23, normal=normal)
         case(5) ! right transonic
            a = (a4 - u4 * eos4%delta()) / (1._R8P + eos4%delta())
            p = p4 * (a / a4) ** eos4%eta()
            r = eos4%density(pressure=p, speed_of_sound=a)
            call fluxes%compute_fluxes_from_primitive(eos=eos4, p=p, r=r, u=-a, normal=normal)
         case(6) ! right supersonic
            call fluxes%compute_fluxes_from_primitive(eos=eos4, p=p4, r=r4, u=u4, normal=normal)
         endselect
      endselect
   endassociate
   endsubroutine compute_fluxes