res Subroutine

   pure subroutine res(t, u, udot, cj, delta, ires, rpar, ipar)
   !! User-defined residuals routine.
   !! It computes the residuals for the 2D discretized heat equation, with zero boundary values.

      real(rk), intent(in) :: t
      real(rk), intent(in) :: u(*)
      real(rk), intent(in) :: udot(*)
      real(rk), intent(in) :: cj
      real(rk), intent(out) :: delta(*)
      integer, intent(inout) :: ires
      real(rk), intent(inout) :: rpar(*)
      integer, intent(inout) :: ipar(*)

      integer :: i, ioff, j, k, m, m2, neq
      real(rk) :: coeff, temx, temy

      ! Set problem constants using IPAR and RPAR.
      neq = ipar(33)
      m = ipar(34)
      coeff = rpar(4)
      m2 = m + 2

      ! Load U into DELTA, in order to set boundary values.
      delta(1:neq) = u(1:neq)

      ! Loop over interior points, and load residual values.
      do k = 1, m
         ioff = m2*k
         do j = 1, m
            i = ioff + j + 1
            temx = u(i - 1) + u(i + 1)
            temy = u(i - m2) + u(i + m2)
            delta(i) = udot(i) - (temx + temy - 4*u(i))*coeff
         end do
      end do

   end subroutine res