1 | ! Rosenbrock integrator with manual time step control |
---|
2 | ! Solve: d/dt c = f(t,c) |
---|
3 | |
---|
4 | ! written by Edwin Spee, CWI, Amsterdam. Last update: July 28, 1997 |
---|
5 | ! email: Edwin.Spee@cwi.nl (http://edwin-spee.mypage.org/) |
---|
6 | ! adapted to KPP-2.1 by Rolf Sander, Max-Planck Institute, Mainz, Germany, 2005 |
---|
7 | ! |
---|
8 | ! Integration method for Ros2: |
---|
9 | ! C_{n+1} = C_n + 3/2 dt k_1 + 1/2 dt k_2 |
---|
10 | ! k_1 = S f(t_n, C_n) |
---|
11 | ! k_2 = S [f(t_{n+1},C_n + dt k_1) - 2 k_1] |
---|
12 | ! |
---|
13 | ! where g = 1.0 + sqrt(0.5_dp), |
---|
14 | ! S = (I - g dt J ) ^ {-1} |
---|
15 | ! with J the Jacobian |
---|
16 | |
---|
17 | MODULE KPP_ROOT_Integrator |
---|
18 | |
---|
19 | USE KPP_ROOT_Precision, ONLY: dp |
---|
20 | |
---|
21 | IMPLICIT NONE |
---|
22 | PUBLIC |
---|
23 | SAVE |
---|
24 | |
---|
25 | ! description of the error numbers IERR |
---|
26 | CHARACTER(LEN=50), PARAMETER, DIMENSION(1) :: IERR_NAMES = (/ & |
---|
27 | 'dummy value ' /) |
---|
28 | |
---|
29 | CONTAINS |
---|
30 | |
---|
31 | ! ************************************************************************** |
---|
32 | |
---|
33 | SUBROUTINE INTEGRATE( TIN, TOUT, & |
---|
34 | ICNTRL_U, RCNTRL_U, ISTATUS_U, RSTATUS_U, IERR_U ) |
---|
35 | |
---|
36 | IMPLICIT NONE |
---|
37 | |
---|
38 | KPP_REAL, INTENT(IN) :: TIN ! TIN - Start Time |
---|
39 | KPP_REAL, INTENT(IN) :: TOUT ! TOUT - End Time |
---|
40 | ! Optional input parameters and statistics |
---|
41 | INTEGER, INTENT(IN), OPTIONAL :: ICNTRL_U(20) |
---|
42 | KPP_REAL, INTENT(IN), OPTIONAL :: RCNTRL_U(20) |
---|
43 | INTEGER, INTENT(OUT), OPTIONAL :: ISTATUS_U(20) |
---|
44 | KPP_REAL, INTENT(OUT), OPTIONAL :: RSTATUS_U(20) |
---|
45 | INTEGER, INTENT(OUT), OPTIONAL :: IERR_U |
---|
46 | |
---|
47 | CALL ROS2(TIN, TOUT) |
---|
48 | |
---|
49 | ! if optional parameters are given for output |
---|
50 | ! use them to store information in them |
---|
51 | IF (PRESENT(ISTATUS_U)) ISTATUS_U(:) = 0 |
---|
52 | IF (PRESENT(RSTATUS_U)) THEN |
---|
53 | RSTATUS_U(:) = 0._dp |
---|
54 | RSTATUS_U(1)=TOUT ! put final time into RSTATUS_U |
---|
55 | ENDIF |
---|
56 | IF (PRESENT(IERR_U)) IERR_U = 1 ! dummy value |
---|
57 | |
---|
58 | END SUBROUTINE INTEGRATE |
---|
59 | |
---|
60 | ! ************************************************************************** |
---|
61 | |
---|
62 | SUBROUTINE ROS2(Tstart,Tend) |
---|
63 | |
---|
64 | USE KPP_ROOT_Jacobian, ONLY: Jac_SP |
---|
65 | USE KPP_ROOT_Global, ONLY: VAR, & ! VARiable species |
---|
66 | FIX, & ! FIXed species |
---|
67 | RCONST ! rate coefficients |
---|
68 | USE KPP_ROOT_JacobianSP, ONLY: LU_DIAG |
---|
69 | USE KPP_ROOT_Function, ONLY: Fun |
---|
70 | USE KPP_ROOT_LinearAlgebra, ONLY: KppDecomp, KppSolve |
---|
71 | USE KPP_ROOT_Parameters, ONLY: NVAR, LU_NONZERO |
---|
72 | |
---|
73 | IMPLICIT NONE |
---|
74 | |
---|
75 | KPP_REAL, INTENT(IN) :: Tstart, Tend |
---|
76 | |
---|
77 | KPP_REAL :: dt ! time step |
---|
78 | KPP_REAL :: k1(NVAR), k2(NVAR), w1(NVAR), g, jvs(LU_NONZERO) |
---|
79 | INTEGER ising, i |
---|
80 | |
---|
81 | dt = Tend - Tstart |
---|
82 | g = 1.0 + SQRT(0.5_dp) |
---|
83 | CALL JAC_sp(VAR, FIX, RCONST, jvs) |
---|
84 | jvs(1:LU_NONZERO) = -g*dt*jvs(1:LU_NONZERO) |
---|
85 | |
---|
86 | ! Rolf von Kuhlmann: |
---|
87 | ! optionally cut this out and replace it by directly addressed statements |
---|
88 | DO i=1,NVAR |
---|
89 | jvs(LU_DIAG(i)) = jvs(LU_DIAG(i)) + 1.0_dp |
---|
90 | END DO |
---|
91 | |
---|
92 | CALL KppDecomp (jvs, ising) |
---|
93 | |
---|
94 | IF (ising /= 0) THEN |
---|
95 | PRINT *,'ising <> 0, dt=',dt |
---|
96 | STOP |
---|
97 | END IF |
---|
98 | |
---|
99 | CALL Fun(VAR, FIX, RCONST, k1 ) |
---|
100 | |
---|
101 | CALL KppSolve (jvs,k1) |
---|
102 | |
---|
103 | DO i = 1,NVAR |
---|
104 | w1(i) = MAX(0.0_dp, VAR(i) + dt * k1(i) ) |
---|
105 | END DO |
---|
106 | |
---|
107 | CALL Fun(w1, FIX, RCONST, k2 ) |
---|
108 | |
---|
109 | k2(1:NVAR) = k2(1:NVAR) - 2.0_dp*k1(1:NVAR) |
---|
110 | CALL KppSolve (jvs,k2) |
---|
111 | DO i = 1,NVAR |
---|
112 | VAR(i) = MAX( 0.0_dp, VAR(i)+1.5_dp*dt*k1(i)+0.5_dp*dt*k2(i) ) |
---|
113 | END DO |
---|
114 | |
---|
115 | END SUBROUTINE ROS2 |
---|
116 | |
---|
117 | ! ************************************************************************** |
---|
118 | |
---|
119 | END MODULE KPP_ROOT_Integrator |
---|