chem/KPP/kpp/kpp-2.1/int/user_contributed/ros2_manual.f90

   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