cinit Subroutine

   pure subroutine cinit(c, cdot, pred_ic, rpar)
   !! This routine computes and loads the vectors of initial values.

      use web_par, only: ax, ay, alpha, beta, dx, dy, ns, np, mx, my, mxns, pi
   
      real(rk), intent(out) :: c(:)
      real(rk), intent(out) :: cdot(:)
      real(rk), intent(in)  :: pred_ic
      real(rk), intent(inout) :: rpar(:)

      integer :: i, ioff, iyoff, jx, jy, npp1
      real(rk) :: argx, argy, fac, t, x, y

      ! Load C
      npp1 = np + 1
      do jy = 1, my
         y = (jy - 1)*dy
         argy = 16*y*y*(ay - y)*(ay - y)
         iyoff = mxns*(jy - 1)
         do jx = 1, mx
            x = (jx - 1)*dx
            argx = 16*x*x*(ax - x)*(ax - x)
            ioff = iyoff + ns*(jx - 1)
            fac = one + alpha*x*y + beta*sin(4*pi*x)*sin(4*pi*y)
            do i = 1, np
               c(ioff + i) = 10.0_rk + i*argx*argy
            end do
            do i = npp1, ns
               c(ioff + i) = pred_ic
            end do
         end do
      end do

      ! Load CDOT
      t = zero
      call f(t, c, cdot, rpar)
      do jy = 1, my
         iyoff = mxns*(jy - 1)
         do jx = 1, mx
            ioff = iyoff + ns*(jx - 1)
            do i = npp1, ns
               cdot(ioff + i) = zero
            end do
         end do
      end do

   end subroutine cinit