jac_rbdpre Subroutine

   subroutine jac_rbdpre(t, u, r0, rblock, r1, rewt, cj, bd, ipbd, ierr)
   !! This routine generates and preprocesses a block-diagonal preconditioner matrix, based on
   !! the part of the Jacobian corresponding to the reaction terms \(R\) of the problem. It 
   !! generates a matrix of the form \(c_J I_d - \partial R/ \partial u\). It calls [[dgefa]]
   !! from LINPACK to do the LU decomposition of each diagonal block. The computation of the
   !! diagonal blocks uses the mesh information in the module variables. One block per spatial
   !! point is computed. The Jacobian elements are generated by difference quotients. 
   !! This routine calls a user routine `rblock` to compute the block `(jx,jy)` of \(R\).

      use linpack, only: dgefa

      real(rk), intent(in) :: t
         !! Current independent variable.
      real(rk), intent(inout) :: u(*)
         !! Current dependent variables.
      real(rk), intent(in) :: r0(*)
         !! Current value of the vector \(R(t,u)\).
      interface
         subroutine rblock(t, jx, jy, uxy, rxy)
            !! User routine that computes a single block of \(R\).
            import :: rk
            real(rk), intent(in) :: t
               !! Current independent variable.
            integer, intent(in) :: jx
               !! Spatial index in x-direction.
            integer, intent(in) :: jy
               !! Spatial index in y-direction.
            real(rk), intent(in) :: uxy(*)
               !! Block of `ns` dependent variables at spatial point `(jx,jy)`.
            real(rk), intent(out) :: rxy(*)
               !! Block `(jx,jy)` of \(R\).
         end subroutine rblock
      end interface
      real(rk), intent(inout) :: r1(*)
         !! Real work space available to this subroutine.
      real(rk), intent(in) :: rewt(*)
         !! Reciprocal error weights.
      real(rk), intent(in) :: cj
         !! Scalar used in forming the system Jacobian.
      real(rk), intent(out) :: bd(*)
         !! LU factors of the diagonal blocks.
      integer, intent(out) :: ipbd(*)
         !! Pivots for the LU factorizations.
      integer, intent(inout) :: ierr
         !! Error flag.
         !! `0`: no error occurred;
         !! `k > 0`: a zero pivot was found at the `k`-th stage in one of the LU factorizations.
      
      real(rk) :: del, dfac, fac, uj
      integer :: i, ibd, idiag, iip, j, j0, js, jx, jy

      ! Make MP calls to RBLOCK to approximate each diagonal block of dR/du.
      dfac = 1e-2_rk
      ibd = 0
      j0 = 0
      do jy = 1, meshy
         do jx = 1, meshx
            ! If R0 has not been set previously as an array of length NEQ, it can
            ! be set here, as an array of length MP, with the call:
            !    CALL RBLOCK (T, JX, JY, U(J0+1), R0)
            ! In this case, change R0(J0+I) below to R0(I).
            do js = 1, mp
               j = j0 + js
               uj = u(j)
               del = max(srur*abs(uj), dfac/rewt(j))
               u(j) = u(j) + del
               fac = -one/del
               call rblock(t, jx, jy, u(j0 + 1), r1)
               do i = 1, mp
                  bd(ibd + i) = (r1(i) - r0(j0 + i))*fac
               end do
               u(j) = uj
               ibd = ibd + mp
            end do
            j0 = j0 + mp
         end do
      end do

      ! Add matrix CJ*Id, and do LU decomposition on blocks.
      ibd = 1
      iip = 1
      do j = 1, meshx*meshy
         idiag = ibd
         do i = 1, mp
            if (i <= mpd) bd(idiag) = bd(idiag) + cj
            idiag = idiag + (mp + 1)
         end do
         call dgefa(bd(ibd), mp, mp, ipbd(iip), ierr)
         if (ierr /= 0) exit
         ibd = ibd + mpsq
         iip = iip + mp
      end do

   end subroutine jac_rbdpre