MODULE KPP_ROOT_Integrator USE KPP_ROOT_Parameters, ONLY : NVAR, NFIX, NREACT USE KPP_ROOT_Global, ONLY : TIME, RCONST, Volume USE KPP_ROOT_Stoichiom USE KPP_ROOT_Stochastic USE KPP_ROOT_Rates IMPLICIT NONE CONTAINS !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE Gillespie(Nevents, T, SCT, NmlcV, NmlcF) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! ! Gillespie stochastic integration ! INPUT: ! Nevents = no. of individual reaction events to be simulated ! SCT = stochastic rate constants ! T = time ! NmlcV, NmlcF = no. of molecules for variable and fixed species ! OUTPUT: ! T = updated time (after Nevents reactions) ! NmlcV = updated no. of molecules for variable species ! ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IMPLICIT NONE KPP_REAL:: T INTEGER :: Nevents INTEGER :: NmlcV(NVAR), NmlcF(NFIX) INTEGER :: i, m, issa REAL :: r1, r2 KPP_REAL :: A(NREACT), SCT(NREACT), x DO issa = 1, Nevents ! Uniformly distributed random numbers CALL RANDOM_NUMBER(r1) CALL RANDOM_NUMBER(r2) ! Avoid log of zero r2 = MAX(r2,1.e-14) ! Propensity vector CALL Propensity ( NmlcV, NmlcF, SCT, A ) ! Cumulative sum of propensities DO i = 2, NREACT A(i) = A(i-1)+A(i); END DO ! Index of next reaction x = r1*A(NREACT) DO i = 1, NREACT IF (A(i)>=x) THEN m = i; EXIT END IF END DO ! Update time with time to next reaction T = T - LOG(r2)/A(NREACT); ! Update state vector CALL MoleculeChange( m, NmlcV ) END DO CONTAINS SUBROUTINE PropensityTemplate( T, NmlcV, NmlcF, Prop ) KPP_REAL, INTENT(IN) :: T INTEGER, INTENT(IN) :: NmlcV(NVAR), NmlcF(NFIX) KPP_REAL, INTENT(OUT) :: Prop(NREACT) KPP_REAL :: Tsave, SCT(NREACT) ! Update the stochastic reaction rates, which may be time dependent Tsave = TIME TIME = T CALL Update_RCONST() CALL StochasticRates( RCONST, Volume, SCT ) CALL Propensity ( NmlcV, NmlcF, SCT, Prop ) TIME = Tsave END SUBROUTINE PropensityTemplate END SUBROUTINE Gillespie END MODULE KPP_ROOT_Integrator