godunov Function

   pure real(rk) function godunov(f, vm, vp, x, t) result(res)
    !!   Monotone Godunov flux. It is less dissipative than the Lax-Friedrichs method, but
    !! computationally more demanding because of the if constructs.
    !!   Source: Equation 2.70, page 21.
    !!
    !! @note
    !!   See note about *elemental* in 'lax_friedrichs'.
      procedure(flux) :: f
        !! flux function, \( f(v, x, t) \)
      real(rk), intent(in) :: vm
        !! left (minus) reconstruction, \( v_{i+1/2}^- \)
      real(rk), intent(in) :: vp
        !! right (plus) reconstruction, \( v_{i+1/2}^+ = v_{(i+1)-1/2}^- \)
      real(rk), intent(in) :: x(:)
        !! \(x\) at flux interface, \( x_{i+1/2} \)
      real(rk), intent(in) :: t
        !! time, \( t \)

      real(rk) :: fm, fp

      fm = f(vm, x, t)
      fp = f(vp, x, t)

      if (vm <= vp) then
         res = min(fm, fp)
      else
         res = max(fm, fp)
      end if

   end function godunov