rhs Subroutine

   pure subroutine rhs(t, v, vdot)
   !! This subroutine computes the *numerical approximation* to the right hand side of:
   !!```
   !!                du(i,t)/dt = -1/dx(i)*( f(u(x(i+1/2),t)) - f(u(x(i-1/2),t)) )
   !!```
   !!   There are two main steps. First, we use the WENO scheme to obtain the reconstructed
   !! values of 'u' at the left and right cell boundaries. Note that, in general, because of
   !! discontinuities, \( u_{i+1/2}^+ \neq u_{(i+1)+1/2}^- \). Second, we use a suitable flux
   !! method (e.g., Godunov, Lax-Friedrichs) to compute the flux from the reconstructed
   !! 'u' values.
      real(rk), intent(in) :: t
         !! time variable
      real(rk), intent(in) :: v(:)
         !! vector(N) with v(x,t) values
      real(rk), intent(out) :: vdot(:)
         !! vector(N) with v'(x,t) values

      real(rk) :: fedges(0:nc), vl(nc), vr(nc)
      integer :: i

      ! Get reconstructed values at cell boundaries
      call myweno%reconstruct(v, vl, vr)

      ! Fluxes at interior cell boundaries
      ! One can use the Lax-Friedrichs or the Godunov method
      do concurrent(i=1:nc - 1)
         !fedges(i) = lax_friedrichs(flux, vr(i), vl(i+1), gx%r(i), t, alpha = 1._rk)
         fedges(i) = godunov(flux, vr(i), vl(i + 1), [gx%right(i)], t)
      end do

      ! Apply problem-specific flux constraints at domain boundaries
      fedges(0) = fedges(1)
      fedges(nc) = fedges(nc - 1)

      ! Evaluate du/dt
      vdot = -(fedges(1:) - fedges(:nc - 1))/gx%width

   end subroutine rhs