source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/gillespie.f90 @ 3800

Last change on this file since 3800 was 2696, checked in by kanani, 7 years ago

Merge of branch palm4u into trunk

File size: 2.5 KB
Line 
1MODULE KPP_ROOT_Integrator
2
3  USE KPP_ROOT_Parameters, ONLY : NVAR, NFIX, NREACT 
4  USE KPP_ROOT_Global, ONLY : TIME, RCONST, Volume 
5  USE KPP_ROOT_Stoichiom 
6  USE KPP_ROOT_Stochastic
7  USE KPP_ROOT_Rates
8  IMPLICIT NONE
9
10CONTAINS
11
12!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
13  SUBROUTINE Gillespie(Nevents, T, SCT, NmlcV, NmlcF)
14!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
15!
16!   Gillespie stochastic integration
17!   INPUT:
18!       Nevents = no. of individual reaction events to be simulated
19!       SCT     = stochastic rate constants
20!       T       = time
21!       NmlcV, NmlcF = no. of molecules for variable and fixed species
22!   OUTPUT:
23!       T       = updated time (after Nevents reactions)
24!       NmlcV   = updated no. of molecules for variable species
25!     
26!
27!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
28      IMPLICIT NONE
29      KPP_REAL:: T
30      INTEGER :: Nevents     
31      INTEGER :: NmlcV(NVAR), NmlcF(NFIX)
32      INTEGER :: i, m, issa
33      REAL   :: r1, r2
34      KPP_REAL :: A(NREACT), SCT(NREACT), x
35   
36      DO issa = 1, Nevents
37
38          ! Uniformly distributed random numbers
39          CALL RANDOM_NUMBER(r1)
40          CALL RANDOM_NUMBER(r2)
41          ! Avoid log of zero
42          r2 = MAX(r2,1.e-14)
43         
44          ! Propensity vector
45          CALL  Propensity ( NmlcV, NmlcF, SCT, A )
46          ! Cumulative sum of propensities
47          DO i = 2, NREACT
48            A(i) = A(i-1)+A(i);
49          END DO
50         
51          ! Index of next reaction
52          x = r1*A(NREACT)
53          DO i = 1, NREACT
54            IF (A(i)>=x) THEN
55              m = i;
56              EXIT
57            END IF
58          END DO
59         
60          ! Update time with time to next reaction
61          T = T - LOG(r2)/A(NREACT);
62
63          ! Update state vector
64          CALL MoleculeChange( m, NmlcV )
65       
66     END DO
67   
68  CONTAINS
69
70    SUBROUTINE PropensityTemplate( T, NmlcV, NmlcF, Prop )
71      KPP_REAL, INTENT(IN)  :: T
72      INTEGER, INTENT(IN)   :: NmlcV(NVAR), NmlcF(NFIX)
73      KPP_REAL, INTENT(OUT) :: Prop(NREACT)
74      KPP_REAL :: Tsave, SCT(NREACT)
75    ! Update the stochastic reaction rates, which may be time dependent
76      Tsave = TIME
77      TIME = T
78      CALL Update_RCONST()
79      CALL StochasticRates( RCONST, Volume, SCT )   
80      CALL Propensity ( NmlcV, NmlcF, SCT, Prop )
81      TIME = Tsave
82    END SUBROUTINE PropensityTemplate   
83   
84  END SUBROUTINE Gillespie
85
86END MODULE KPP_ROOT_Integrator
Note: See TracBrowser for help on using the repository browser.