mpf Subroutine

   pure subroutine mpf(u, c, kwee, k25, k25p, print_every, tout, root)
   !! If `root` is not zero, then it returns value of time when `u=root` in `tout`. Else, runs
   !! until `tout` and returns the value in `u`. If `print_every` is non-zero then the solution
   !! is printed every `print_every` time units or every `h` (which ever is greater).
   !!
   !! This routine is not meant to be precise, it is only intended to be good enough for
   !! providing a working example of ODRPACK95 with bounds. 4th order Runge Kutta and linear
   !! interpolation are used for numerical integration and root finding, respectively.
   !!
   !! `u`: MPF
   !! `c`: Total Cyclin
   !! `kwee`, `k25`, `k25p`: Model parameters corresponding to `beta(1:3)`

      real(kind=wp), intent(out) :: u
      real(kind=wp), intent(in) :: c, kwee, k25, k25p, print_every, root
      real(kind=wp), intent(inout) :: tout
      real(kind=wp), parameter :: h = 1.0E-1_wp
      real(kind=wp) :: last_print, last_u, last_t, t
      real(kind=wp) :: k1, k2, k3, k4

      u = zero
      t = zero

      last_print = zero
      if (print_every > zero) then
         !write (*, *) t, u
      end if

      do while (t < tout)
         last_t = t
         last_u = u
         k1 = h*dudt(u, c, kwee, k25, k25p)
         k2 = h*dudt(u + k1/2, c, kwee, k25, k25p)
         k3 = h*dudt(u + k2/2, c, kwee, k25, k25p)
         k4 = h*dudt(u + k3, c, kwee, k25, k25p)
         u = u + (k1 + 2*k2 + 2*k3 + k4)/6
         t = t + h
         if (t >= print_every + last_print .and. print_every > zero) &
            then
            !write (*, *) t, m
            last_print = t
         end if
         if (root > zero) then
            if (last_u <= root .and. root < u) then
               tout = (t - last_t)/(u - last_u)*(root - last_u) + last_t
               return
            end if
         end if
      end do

   contains

      ! Equation from Zwolak et al. 2001.
      real(kind=wp) pure function dudt(u_, c_, kwee_, k25_, k25p_) result(res)
         real(kind=wp), intent(in) :: u_, c_, kwee_, k25_, k25p_
         res = kwee_*u_ + (k25_ + k25p_*u_**2)*(c_ - u_)
      end function dudt

   end subroutine mpf