rhs Subroutine

   subroutine rhs(t, v, vdot)
   !! This subroutine computes the *numerical approximation* to the right hand side of:
   !!```
   !!     du(i,j,t)/dt = -1/dx1(i)*( f1(u(x1(i+1/2),x2(j),t)) - f1(u(x1(i-1/2),x2(j),t)) )
   !!                    -1/dx2(j)*( f2(u(x1(i),x2(j+1/2),t)) - f2(u(x1(i),x2(j-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 along each dimension. 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(z,t) values
      real(rk), intent(out) :: vdot(:)
         !! vector(N) with v'(z,t) values

      real(rk) :: vl1(nc(1)), vl2(nc(2)), vr1(nc(1)), vr2(nc(2))
      real(rk) :: fedges1(0:nc(1), 0:nc(2)), fedges2(0:nc(2), 0:nc(1))
      integer :: i, j

      ! Fluxes along x1 and x2 at interior cell boundaries
      !$omp parallel
      !$omp do private(vl1, vr1)
      do j = 1, nc(2)
         call myweno(1)%reconstruct(v(1+(j-1)*nc(1):j*nc(1)), vl1, vr1)
         do i = 1, nc(1) - 1
            fedges1(i, j) = godunov(flux1, vr1(i), vl1(i + 1), &
                                    [gx(1)%right(i), gx(2)%center(j)], t)
         end do
      end do
      !$omp end do nowait
      !$omp do private(vl2, vr2)
      do i = 1, nc(1)
         call myweno(2)%reconstruct(v(i:i+(nc(2)-1)*nc(1):nc(1)), vl2, vr2)
         do j = 1, nc(2) - 1
            fedges2(j, i) = godunov(flux2, vr2(j), vl2(j + 1), &
                                    [gx(1)%center(i), gx(2)%right(j)], t)
         end do
      end do
      !$omp end do
      !$omp end parallel

      ! Apply problem-specific flux constraints at domain boundaries
      fedges1(0, :) = 0
      fedges1(nc(1), :) = 0
      fedges2(0, :) = 0
      fedges2(nc(2), :) = 0

      ! Evaluate du/dt
      do concurrent(i=1:nc(1), j=1:nc(2))
         vdot((j - 1)*nc(1) + i) = &
            - (fedges1(i, j) - fedges1(i - 1, j))/gx(1)%width(i) &
            - (fedges2(j, i) - fedges2(j - 1, i))/gx(2)%width(j)
      end do

   end subroutine rhs