pure subroutine res(t, u, uprime, cj, delta, ires, rpar, ipar)
!! User-supplied residuals subroutine.
!! 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) :: uprime(*)
real(rk), intent(in) :: cj
real(rk), intent(out) :: delta(*)
integer, intent(inout) :: ires
real(rk), intent(in) :: rpar(*)
integer, intent(in) :: 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) = uprime(i) - (temx + temy - 4*u(i))*coeff
end do
end do
end subroutine res