setup_rbgpre Subroutine

   subroutine setup_rbgpre(mx, my, ns, nsd, nxg, nyg, lid, iwork)
   !! This routine sets the mesh and block-grouping parameters needed to use the routines
   !! [[jac_rbgpre]] and [[psol_rbgpre]], assuming a 2D rectangular problem with uniform
   !! rectangular block-grouping. Additionally, it loads the lengths `lrwp` and `liwp` in array
   !! `iwork`.

      integer, intent(in) :: mx
         !! Number of mesh points in x-direction.
      integer, intent(in) :: my
         !! Number of mesh points in y-direction.
      integer, intent(in) :: ns
         !! Number of PDE variables, the size of each block in the block-diagonal preconditioner
         !! matrix \(P_R\).
      integer, intent(in) :: nsd
         !! Number of differential PDE variables. In the DAE system, the first `nsd` variables
         !! at each spatial point have time derivatives, and the remaining `(ns - nsd)` do not.
      integer, intent(in) :: nxg
         !! Number of groups in x-direction.
      integer, intent(in) :: nyg
         !! Number of groups in y-direction.
      integer, intent(in) :: lid
         !! Flag indicating whether to load the ID array in `iwork`. If `lid > 0`, set the ID
         !! array in `iwork`, indicating which components are differential and which are algebraic.
         !! This value is required if either `info(11) = 1` or `info(16) = 1`, in which case
         !! set `lid = 40` or `lid = 40 + neq`, depending on the value of the constraint option
         !! `info(10)`. Otherwise, set `lid = 0`.
      integer, intent(inout) :: iwork(*)
         !! Integer work array.

      integer :: i, i0, jx, jy, maxm

      ! Load the modules variables.
      srur = sqrt(epsilon(one))
      mp = ns
      mpd = nsd
      mpsq = ns*ns
      meshx = mx
      meshy = my
      mxmp = meshx*mp
      ngx = nxg
      ngy = nyg
      ngrp = ngx*ngy

      maxm = max(meshx, meshy)
      allocate(jgx(maxm + 1), jgy(maxm + 1), jigx(maxm), jigy(maxm), jxr(maxm), jyr(maxm))
   
      ! Call set_grouping for each mesh direction to load grouping arrays.
      call set_grouping(meshx, ngx, jgx, jigx, jxr)
      call set_grouping(meshy, ngy, jgy, jigy, jyr)

      ! Set the sizes of the preconditioning storage space segments in RWORK and IWORK.
      iwork(27) = mpsq*ngrp
      iwork(28) = mp*ngrp

      ! If LID > 0, set the ID array in IWORK.
      if (lid == 0) return
      i0 = lid
      do jy = 1, my
         do jx = 1, mx
            do i = 1, mpd
               iwork(i0 + i) = 1
            end do
            do i = mpd + 1, mp
               iwork(i0 + i) = -1
               i0 = i0 + mp
            end do
         end do
      end do

   end subroutine setup_rbgpre