alloc_state Subroutine

   pure subroutine alloc_state(self, n, dima, stat)
   !! Initializes [[nf_state]], i.e., (re)allocates all allocatable arrays and sets them to
   !! zero.

      class(nf_state), intent(inout) :: self
         !! State object.
      integer, intent(in)  :: n
         !! Problem dimension.
      integer, intent(in)  :: dima
         !! Dimension of the parameter vector `a`.
      integer, intent(out), optional   :: stat
         !! Error status of the allocation.

      integer :: ierr(6)

      if (n <= 0) then
         error stop "Error: 'n' must be positive in hompack_nf_state%alloc()."
      end if

      if (dima <= 0) then
         error stop "Error: 'dima' must be positive in hompack_nf_state%alloc()."
      end if

      if (present(stat)) stat = 0

      ! Deallocate any previously allocated arrays
      self%n = 0
      if (allocated(self%y)) deallocate (self%y)
      if (allocated(self%yold)) deallocate (self%yold)
      if (allocated(self%yp)) deallocate (self%yp)
      if (allocated(self%ypold)) deallocate (self%ypold)
      if (allocated(self%a)) deallocate (self%a)

      ! Allocate/initialize state arrays
      allocate (self%y(n + 1), source=zero, stat=ierr(1))
      allocate (self%yold(n + 1), source=zero, stat=ierr(2))
      allocate (self%yp(n + 1), source=zero, stat=ierr(3))
      allocate (self%ypold(n + 1), source=zero, stat=ierr(4))
      allocate (self%a(dima), source=zero, stat=ierr(5))

      ! Initialize workspace
      call self%workspace%alloc(n, stat=ierr(6))

      if (any(ierr /= 0)) then
         if (present(stat)) then
            stat = ierr(findloc(ierr /= 0, .true., dim=1))
         else
            error stop "Error: Allocation failed in nf_state%alloc()."
         end if
      end if

      ! Set problem dimension now that allocation is successful
      self%n = n

   end subroutine alloc_state