[2696] | 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 |
---|