mpf Subroutine

   pure subroutine mpf(m, c, kwee, k25, k25p, print_every, tout, root)
   !! If ROOT is not zero then returns value of time when M==ROOT in TOUT.  Else,
   !! runs until TOUT and returns value in M.  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.
   !!
   !! M - MPF
   !! C - Total Cyclin
   !! KWEE, K25, K25P - Model parameters (BETA(1:3))

      real(kind=wp), intent(out) :: m
      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_m, last_t, t
      real(kind=wp) :: k1, k2, k3, k4

      m = zero
      t = zero

      last_print = zero
      if (print_every .gt. zero) then
         !write (*, *) t, m
      end if

      do while (t .lt. tout)
         last_t = t
         last_m = m
         k1 = h*dmdt(m, c, kwee, k25, k25p)
         k2 = h*dmdt(m + k1/2, c, kwee, k25, k25p)
         k3 = h*dmdt(m + k2/2, c, kwee, k25, k25p)
         k4 = h*dmdt(m + k3, c, kwee, k25, k25p)
         m = m + (k1 + 2*k2 + 2*k3 + k4)/6
         t = t + h
         if (t .ge. print_every + last_print .and. print_every .gt. zero) &
            then
            !write (*, *) t, m
            last_print = t
         end if
         if (root .gt. zero) then
            if (last_m .le. root .and. root .lt. m) then
               tout = (t - last_t)/(m - last_m)*(root - last_m) + last_t
               return
            end if
         end if
      end do

   contains

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

   end subroutine mpf