source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/user_contributed/ros2_manual.f90 @ 2696

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

Merge of branch palm4u into trunk

File size: 3.4 KB
Line 
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
17MODULE 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
29CONTAINS
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
119END MODULE KPP_ROOT_Integrator
Note: See TracBrowser for help on using the repository browser.