pure subroutine f(t, c, cprime, rpar)
!! This routine computes the right-hand sides of all the equations and returns them in the
!! array `cprime`. The interaction rates are computed by calls to [[rates]], and these are
!! saved in `rpar` for later use in preconditioning.
use web_par, only: cox, coy, ns, mx, my, mxns
real(rk), intent(in) :: t
real(rk), intent(in) :: c(*)
real(rk), intent(out) :: cprime(*)
real(rk), intent(inout) :: rpar(*)
integer :: i, ic, ici, idxl, idxu, idyl, idyu, iyoff, jx, jy
real(rk) :: dcxli, dcxui, dcyli, dcyui
do jy = 1, my
iyoff = mxns*(jy - 1)
idyu = mxns
if (jy == my) idyu = -mxns
idyl = mxns
if (jy == 1) idyl = -mxns
do jx = 1, mx
ic = iyoff + ns*(jx - 1) + 1
! Get interaction rates at one point (X,Y)
call rates(t, jx, jy, c(ic), rpar(ic))
idxu = ns
if (jx == mx) idxu = -ns
idxl = ns
if (jx == 1) idxl = -ns
do i = 1, ns
ici = ic + i - 1
! Do differencing in Y
dcyli = c(ici) - c(ici - idyl)
dcyui = c(ici + idyu) - c(ici)
! Do differencing in X
dcxli = c(ici) - c(ici - idxl)
dcxui = c(ici + idxu) - c(ici)
! Collect terms and load CPRIME elements
cprime(ici) = coy(i)*(dcyui - dcyli) + cox(i)*(dcxui - dcxli) + rpar(ici)
end do
end do
end do
end subroutine f