Juan 8/12/2016: add management of LEN_HREC in MNH & SURFEX
[MNH-git_open_source-lfs.git] / src / MNH / mode_RBK90_Integrator.f90
1 !MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
2 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
3 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
4 !MNH_LIC for details. version 1.
5 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
6
7 ! Numerical Integrator (Time-Stepping) File
8
9 ! Generated by KPP-2.2 symbolic chemistry Kinetics PreProcessor
10 !       (http://www.cs.vt.edu/~asandu/Software/KPP)
11 ! KPP is distributed under GPL, the general public licence
12 !       (http://www.gnu.org/copyleft/gpl.html)
13 ! (C) 1995-1997, V. Damian & A. Sandu, CGRER, Univ. Iowa
14 ! (C) 1997-2005, A. Sandu, Michigan Tech, Virginia Tech
15 !     With important contributions from:
16 !        M. Damian, Villanova University, USA
17 !        R. Sander, Max-Planck Institute for Chemistry, Mainz, Germany
18
19 ! File                 : RBK90_Integrator.f90
20 ! Time                 : Mon Apr 16 10:11:55 2007
21 ! Working directory    : /home/pinjp/chimie_num/kpp/kpp-2.2.1.December2006/my-test-NumRec
22 ! Equation file        : RBK90.kpp
23 ! Output root filename : RBK90
24
25 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
26
27
28 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
29
30 ! INTEGRATE - Integrator routine
31 !   Arguments :
32 !      TIN       - Start Time for Integration
33 !      TOUT      - End Time for Integration
34
35 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
36
37 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
38 !  Rosenbrock - Implementation of several Rosenbrock methods:             !
39 !               * Ros1                                                    !
40 !               * Ros2                                                    !
41 !               * Ros3                                                    !
42 !               * Ros4                                                    !
43 !               * Rodas3                                                  !
44 !               * Rodas4                                                  !
45 !  By default the code employs the KPP sparse linear algebra routines     !
46 !  Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK)        !
47 !                                                                         !
48 !    (C)  Adrian Sandu, August 2004                                       !
49 !    Virginia Polytechnic Institute and State University                  !
50 !    Contact: sandu@cs.vt.edu                                             !
51 !    Revised by Philipp Miehe and Adrian Sandu, May 2006                  !                               !
52 !    This implementation is part of KPP - the Kinetic PreProcessor        !
53 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
54
55 MODULE MODE_RBK90_Integrator
56
57   USE MODD_RBK90_JacobianSP_n, ONLY: LU_DIM_SPECIES
58   USE MODD_RBK90_Parameters_n, ONLY: NVAR
59   USE MODD_RBK90_Global_n,     ONLY: STEPMIN
60   IMPLICIT NONE
61   PUBLIC
62   SAVE
63   
64 !~~~>  Statistics on the work performed by the Rosenbrock method
65   INTEGER, PARAMETER :: Nfun=1, Njac=2, Nstp=3, Nacc=4, &
66                         Nrej=5, Ndec=6, Nsol=7, Nsng=8, &
67                         Ntexit=1, Nhexit=2, Nhnew = 3
68
69 CONTAINS
70
71 SUBROUTINE CH_ROSENBROCK( TIN, TOUT, SPECIES, NSPECIES, KMI, KVECNPT, &
72          ATOL, RTOL, ICNTRL_U, RCNTRL_U, ISTATUS_U, RSTATUS_U, IERR_U )
73
74    IMPLICIT NONE
75
76    REAL, INTENT(IN) :: TIN  ! Start Time
77    REAL, INTENT(IN) :: TOUT ! End Time
78 !
79 ! NSPECIES - Number of Variable species
80 INTEGER, INTENT(IN) :: NSPECIES
81 !
82 INTEGER, INTENT(IN) :: KVECNPT
83 INTEGER, INTENT(IN) :: KMI      ! model number
84 !
85 ! SPECIES - Concentrations of variable species (global)
86 REAL, INTENT(INOUT) :: SPECIES(KVECNPT,NSPECIES)
87 ! ATOL - Absolute tolerance
88 REAL,INTENT(IN)  :: ATOL(NSPECIES)
89 ! RTOL - Relative tolerance
90 REAL, INTENT(IN) :: RTOL(NSPECIES)
91 !
92    ! Optional input parameters and statistics
93    INTEGER,       INTENT(IN),  OPTIONAL :: ICNTRL_U(20)
94    REAL, INTENT(IN),  OPTIONAL :: RCNTRL_U(20)
95    INTEGER,       INTENT(OUT), OPTIONAL :: ISTATUS_U(20)
96    REAL, INTENT(OUT), OPTIONAL :: RSTATUS_U(20)
97    INTEGER,       INTENT(OUT), OPTIONAL :: IERR_U
98
99    REAL :: RCNTRL(20), RSTATUS(20)
100    INTEGER       :: ICNTRL(20), ISTATUS(20), IERR
101
102    INTEGER, SAVE :: Ntotal = 0
103    INTEGER       :: IDIM_SOLVER ! KVECNPT times the chemical system size
104    INTEGER       :: JL, JLSHIFT   ! Loop indexes
105    INTEGER       :: ISPECIES
106    REAL, DIMENSION(NVAR) :: ZCONC ! Vectorized SPECIES
107    REAL, DIMENSION(NVAR) :: ZATOL ! Vectorized ATOL
108    REAL, DIMENSION(NVAR) :: ZRTOL ! Vectorized RTOL
109
110    ICNTRL(:)  = 0
111    RCNTRL(:)  = 0.0
112    ISTATUS(:) = 0
113    RSTATUS(:) = 0.0
114
115 !JPP    !~~~> fine-tune the integrator:
116 !JPP   ICNTRL(1) = 0    ! 0 - non-autonomous, 1 - autonomous
117 !JPP   ICNTRL(2) = 0    ! 0 - vector tolerances, 1 - scalars
118
119    ! If optional parameters are given, and if they are >0, 
120    ! then they overwrite default settings. 
121    IF (PRESENT(ICNTRL_U)) THEN
122      WHERE(ICNTRL_U(:) > 0) ICNTRL(:) = ICNTRL_U(:)
123    END IF
124    IF (PRESENT(RCNTRL_U)) THEN
125      WHERE(RCNTRL_U(:) > 0) RCNTRL(:) = RCNTRL_U(:)
126    END IF
127
128    JLSHIFT = 0
129    DO JL = 1,KVECNPT
130      ISPECIES = LU_DIM_SPECIES(JL)
131      ZCONC(JLSHIFT+1:JLSHIFT+ISPECIES) = SPECIES(JL,1:ISPECIES)
132      ZATOL(JLSHIFT+1:JLSHIFT+ISPECIES) = ATOL(1:ISPECIES)
133      ZRTOL(JLSHIFT+1:JLSHIFT+ISPECIES) = RTOL(1:ISPECIES)
134      JLSHIFT = JLSHIFT + ISPECIES
135    END DO
136    IDIM_SOLVER = NVAR ! Note: NVAR is the final value of JLSHIFT
137
138    CALL Rosenbrock(IDIM_SOLVER,ZCONC,TIN,TOUT, &
139         KMI,KVECNPT,                           &
140         ZATOL,ZRTOL,                           &
141         RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR)
142
143    JLSHIFT = 0
144    DO JL = 1,KVECNPT
145      ISPECIES = LU_DIM_SPECIES(JL)
146      SPECIES(JL,1:ISPECIES) = ZCONC(JLSHIFT+1:JLSHIFT+ISPECIES)
147      JLSHIFT = JLSHIFT + ISPECIES
148    END DO
149
150    !~~~> Debug option: show no of steps
151    ! Ntotal = Ntotal + ISTATUS(Nstp)
152    ! PRINT*,'NSTEPS=',ISTATUS(Nstp),' (',Ntotal,')','  O3=', SPECIES(ind_O3)
153
154    STEPMIN = RSTATUS(Nhexit)
155    ! if optional parameters are given for output they 
156    ! are updated with the return information
157    IF (PRESENT(ISTATUS_U)) ISTATUS_U(:) = ISTATUS(:)
158    IF (PRESENT(RSTATUS_U)) RSTATUS_U(:) = RSTATUS(:)
159    IF (PRESENT(IERR_U))    IERR_U       = IERR
160
161 END SUBROUTINE CH_ROSENBROCK
162
163 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
164 SUBROUTINE Rosenbrock(N,Y,Tstart,Tend, &
165            KMI,KVECNPT,                &
166            AbsTol,RelTol,              &
167            RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR)
168 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
169 !
170 !    Solves the system y'=F(t,y) using a Rosenbrock method defined by:
171 !
172 !     G = 1/(H*gamma(1)) - Jac(t0,Y0)
173 !     T_i = t0 + Alpha(i)*H
174 !     Y_i = Y0 + \sum_{j=1}^{i-1} A(i,j)*K_j
175 !     G * K_i = Fun( T_i, Y_i ) + \sum_{j=1}^S C(i,j)/H * K_j +
176 !         gamma(i)*dF/dT(t0, Y0)
177 !     Y1 = Y0 + \sum_{j=1}^S M(j)*K_j
178 !
179 !    For details on Rosenbrock methods and their implementation consult:
180 !      E. Hairer and G. Wanner
181 !      "Solving ODEs II. Stiff and differential-algebraic problems".
182 !      Springer series in computational mathematics, Springer-Verlag, 1996.
183 !    The codes contained in the book inspired this implementation.
184 !
185 !    (C)  Adrian Sandu, August 2004
186 !    Virginia Polytechnic Institute and State University
187 !    Contact: sandu@cs.vt.edu
188 !    Revised by Philipp Miehe and Adrian Sandu, May 2006                  
189 !    This implementation is part of KPP - the Kinetic PreProcessor
190 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
191 !
192 !~~~>   INPUT ARGUMENTS:
193 !
194 !-     Y(N)    = vector of initial conditions (at T=Tstart)
195 !-    [Tstart,Tend]  = time range of integration
196 !     (if Tstart>Tend the integration is performed backwards in time)
197 !-    RelTol, AbsTol = user precribed accuracy
198 !- SUBROUTINE  Fun( T, Y, Ydot ) = ODE function,
199 !                       returns Ydot = Y' = F(T,Y)
200 !- SUBROUTINE  Jac( T, Y, Jcb ) = Jacobian of the ODE function,
201 !                       returns Jcb = dFun/dY
202 !-    ICNTRL(1:20)    = integer inputs parameters
203 !-    RCNTRL(1:20)    = real inputs parameters
204 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
205 !
206 !~~~>     OUTPUT ARGUMENTS:
207 !
208 !-    Y(N)    -> vector of final states (at T->Tend)
209 !-    ISTATUS(1:20)   -> integer output parameters
210 !-    RSTATUS(1:20)   -> real output parameters
211 !-    IERR            -> job status upon return
212 !                        success (positive value) or
213 !                        failure (negative value)
214 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
215 !
216 !~~~>     INPUT PARAMETERS:
217 !
218 !    Note: For input parameters equal to zero the default values of the
219 !       corresponding variables are used.
220 !
221 !    ICNTRL(1) = 1: F = F(y)   Independent of T (AUTONOMOUS)
222 !              = 0: F = F(t,y) Depends on T (NON-AUTONOMOUS)
223 !
224 !    ICNTRL(2) = 0: AbsTol, RelTol are N-dimensional vectors
225 !              = 1: AbsTol, RelTol are scalars
226 !
227 !    ICNTRL(3)  -> selection of a particular Rosenbrock method
228 !        = 0 :    Ros1
229 !        = 1 :    Ros2
230 !        = 2 :    Ros3
231 !        = 3 :    Ros4
232 !        = 4 :    Rodas3
233 !        = 5 :    Rodas4
234 !
235 !    ICNTRL(4)  -> maximum number of integration steps
236 !        For ICNTRL(4)=0) the default value of 100000 is used
237 !
238 !    RCNTRL(1)  -> Hmin, lower bound for the integration step size
239 !          It is strongly recommended to keep Hmin = ZERO
240 !    RCNTRL(2)  -> Hmax, upper bound for the integration step size
241 !    RCNTRL(3)  -> Hstart, starting value for the integration step size
242 !
243 !    RCNTRL(4)  -> FacMin, lower bound on step decrease factor (default=0.2)
244 !    RCNTRL(5)  -> FacMax, upper bound on step increase factor (default=6)
245 !    RCNTRL(6)  -> FacRej, step decrease factor after multiple rejections
246 !                          (default=0.1)
247 !    RCNTRL(7)  -> FacSafe, by which the new step is slightly smaller
248 !         than the predicted value  (default=0.9)
249 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
250 !
251 !
252 !    OUTPUT ARGUMENTS:
253 !    -----------------
254 !
255 !    T           -> T value for which the solution has been computed
256 !                     (after successful return T=Tend).
257 !
258 !    Y(N)        -> Numerical solution at T
259 !
260 !    IDID        -> Reports on successfulness upon return:
261 !                    = 1 for success
262 !                    < 0 for error (value equals error code)
263 !
264 !    ISTATUS(1)  -> No. of function calls
265 !    ISTATUS(2)  -> No. of jacobian calls
266 !    ISTATUS(3)  -> No. of steps
267 !    ISTATUS(4)  -> No. of accepted steps
268 !    ISTATUS(5)  -> No. of rejected steps (except at very beginning)
269 !    ISTATUS(6)  -> No. of LU decompositions
270 !    ISTATUS(7)  -> No. of forward/backward substitutions
271 !    ISTATUS(8)  -> No. of singular matrix decompositions
272 !
273 !    RSTATUS(1)  -> Texit, the time corresponding to the
274 !                     computed Y upon return
275 !    RSTATUS(2)  -> Hexit, last accepted step before exit
276 !    RSTATUS(3)  -> Hnew, last predicted step (not yet taken)
277 !                   For multiple restarts, use Hnew as Hstart 
278 !                     in the subsequent run
279 !
280 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
281
282   USE MODD_RBK90_Parameters_n
283   USE MODE_RBK90_LinearAlgebra
284   IMPLICIT NONE
285
286 !~~~>  Arguments
287    INTEGER,       INTENT(IN)    :: N
288    REAL, INTENT(INOUT) :: Y(N)
289    REAL, INTENT(IN)    :: Tstart,Tend
290 !
291 INTEGER, INTENT(IN) :: KVECNPT
292 INTEGER, INTENT(IN) :: KMI      ! model number
293 !
294    REAL, INTENT(IN)    :: AbsTol(N),RelTol(N)
295    INTEGER,       INTENT(IN)    :: ICNTRL(20)
296    REAL, INTENT(IN)    :: RCNTRL(20)
297    INTEGER,       INTENT(INOUT) :: ISTATUS(20)
298    REAL, INTENT(INOUT) :: RSTATUS(20)
299    INTEGER, INTENT(OUT)   :: IERR
300 !~~~>  Parameters of the Rosenbrock method, up to 6 stages
301    INTEGER ::  ros_S, rosMethod
302    INTEGER, PARAMETER :: RS1=0, RS2=1, RS3=2, RS4=3, RD3=4, RD4=5
303    REAL :: ros_A(15), ros_C(15), ros_M(6), ros_E(6), &
304                     ros_Alpha(6), ros_Gamma(6), ros_ELO
305    LOGICAL :: ros_NewF(6)
306    CHARACTER(LEN=LEN_HREC) :: ros_Name
307 !~~~>  Local variables
308    REAL :: Roundoff, FacMin, FacMax, FacRej, FacSafe
309    REAL :: Hmin, Hmax, Hstart
310    REAL :: Texit
311    INTEGER       :: i, UplimTol, Max_no_steps
312    LOGICAL       :: Autonomous, VectorTol
313 !~~~>   Parameters
314    REAL, PARAMETER :: ZERO = 0.0, ONE  = 1.0
315    REAL, PARAMETER :: DeltaMin = 1.0E-5
316
317 !~~~>  Initialize statistics
318    ISTATUS(1:8) = 0
319    RSTATUS(1:3) = ZERO
320
321 !~~~>  Autonomous or time dependent ODE. Default is time dependent.
322    Autonomous = .NOT.(ICNTRL(1) == 0)
323
324 !~~~>  For Scalar tolerances (ICNTRL(2).NE.0)  the code uses AbsTol(1) and RelTol(1)
325 !   For Vector tolerances (ICNTRL(2) == 0) the code uses AbsTol(1:N) and RelTol(1:N)
326    IF (ICNTRL(2) == 0) THEN
327       VectorTol = .TRUE.
328       UplimTol  = N
329    ELSE
330       VectorTol = .FALSE.
331       UplimTol  = 1
332    END IF
333
334 !~~~>   Initialize the particular Rosenbrock method selected
335    SELECT CASE (ICNTRL(3))
336      CASE (0)
337        CALL Ros1
338      CASE (1)
339        CALL Ros2
340      CASE (2)
341        CALL Ros3
342      CASE (3)
343        CALL Ros4
344      CASE (4)
345        CALL Rodas3
346      CASE (5)
347        CALL Rodas4
348      CASE DEFAULT
349        PRINT * , 'Unknown Rosenbrock method: ICNTRL(3)=',ICNTRL(3) 
350        CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR)
351        RETURN
352    END SELECT
353
354 !~~~>   The maximum number of steps admitted
355    IF (ICNTRL(4) == 0) THEN
356       Max_no_steps = 200000
357    ELSEIF (ICNTRL(4) > 0) THEN
358       Max_no_steps=ICNTRL(4)
359    ELSE
360       PRINT * ,'User-selected max no. of steps: ICNTRL(4)=',ICNTRL(4)
361       CALL ros_ErrorMsg(-1,Tstart,ZERO,IERR)
362       RETURN
363    END IF
364
365 !~~~>  Unit roundoff (1+Roundoff>1)
366    Roundoff = WLAMCH('E')
367
368 !~~~>  Lower bound on the step size: (positive value)
369    IF (RCNTRL(1) == ZERO) THEN
370       Hmin = ZERO
371    ELSEIF (RCNTRL(1) > ZERO) THEN
372       Hmin = RCNTRL(1)
373    ELSE
374       PRINT * , 'User-selected Hmin: RCNTRL(1)=', RCNTRL(1)
375       CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
376       RETURN
377    END IF
378 !~~~>  Upper bound on the step size: (positive value)
379    IF (RCNTRL(2) == ZERO) THEN
380       Hmax = ABS(Tend-Tstart)
381    ELSEIF (RCNTRL(2) > ZERO) THEN
382       Hmax = MIN(ABS(RCNTRL(2)),ABS(Tend-Tstart))
383    ELSE
384       PRINT * , 'User-selected Hmax: RCNTRL(2)=', RCNTRL(2)
385       CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
386       RETURN
387    END IF
388 !~~~>  Starting step size: (positive value)
389    IF (RCNTRL(3) == ZERO) THEN
390       Hstart = MAX(Hmin,DeltaMin)
391    ELSEIF (RCNTRL(3) > ZERO) THEN
392       Hstart = MIN(ABS(RCNTRL(3)),ABS(Tend-Tstart))
393    ELSE
394       PRINT * , 'User-selected Hstart: RCNTRL(3)=', RCNTRL(3)
395       CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
396       RETURN
397    END IF
398 !~~~>  Step size can be changed s.t.  FacMin < Hnew/Hold < FacMax
399    IF (RCNTRL(4) == ZERO) THEN
400       FacMin = 0.2
401    ELSEIF (RCNTRL(4) > ZERO) THEN
402       FacMin = RCNTRL(4)
403    ELSE
404       PRINT * , 'User-selected FacMin: RCNTRL(4)=', RCNTRL(4)
405       CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
406       RETURN
407    END IF
408    IF (RCNTRL(5) == ZERO) THEN
409       FacMax = 6.0
410    ELSEIF (RCNTRL(5) > ZERO) THEN
411       FacMax = RCNTRL(5)
412    ELSE
413       PRINT * , 'User-selected FacMax: RCNTRL(5)=', RCNTRL(5)
414       CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
415       RETURN
416    END IF
417 !~~~>   FacRej: Factor to decrease step after 2 succesive rejections
418    IF (RCNTRL(6) == ZERO) THEN
419       FacRej = 0.1
420    ELSEIF (RCNTRL(6) > ZERO) THEN
421       FacRej = RCNTRL(6)
422    ELSE
423       PRINT * , 'User-selected FacRej: RCNTRL(6)=', RCNTRL(6)
424       CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
425       RETURN
426    END IF
427 !~~~>   FacSafe: Safety Factor in the computation of new step size
428    IF (RCNTRL(7) == ZERO) THEN
429       FacSafe = 0.9
430    ELSEIF (RCNTRL(7) > ZERO) THEN
431       FacSafe = RCNTRL(7)
432    ELSE
433       PRINT * , 'User-selected FacSafe: RCNTRL(7)=', RCNTRL(7)
434       CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
435       RETURN
436    END IF
437 !~~~>  Check if tolerances are reasonable
438     DO i=1,UplimTol
439       IF ( (AbsTol(i) <= ZERO) .OR. (RelTol(i) <= 10.0*Roundoff) &
440          .OR. (RelTol(i) >= 1.0) ) THEN
441         PRINT * , ' AbsTol(',i,') = ',AbsTol(i)
442         PRINT * , ' RelTol(',i,') = ',RelTol(i)
443         CALL ros_ErrorMsg(-5,Tstart,ZERO,IERR)
444         RETURN
445       END IF
446     END DO
447
448
449 !~~~>  CALL Rosenbrock method
450    CALL ros_Integrator(Y, Tstart, Tend, Texit,   &
451         KMI,KVECNPT,                             &
452         AbsTol, RelTol,                          &
453 !  Integration parameters
454         Autonomous, VectorTol, Max_no_steps,     &
455         Roundoff, Hmin, Hmax, Hstart,            &
456         FacMin, FacMax, FacRej, FacSafe,         &
457 !  Error indicator
458         IERR)
459
460 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
461 CONTAINS !  SUBROUTINES internal to Rosenbrock
462 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
463    
464 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
465  SUBROUTINE ros_ErrorMsg(Code,T,H,IERR)
466 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
467 !    Handles all error messages
468 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
469   
470    REAL, INTENT(IN) :: T, H
471    INTEGER, INTENT(IN)  :: Code
472    INTEGER, INTENT(OUT) :: IERR
473    
474    IERR = Code
475    PRINT * , &
476      'Forced exit from Rosenbrock due to the following error:' 
477      
478    SELECT CASE (Code)
479     CASE (-1)    
480       PRINT * , '--> Improper value for maximal no of steps'
481     CASE (-2)    
482       PRINT * , '--> Selected Rosenbrock method not implemented'
483     CASE (-3)    
484       PRINT * , '--> Hmin/Hmax/Hstart must be positive'
485     CASE (-4)    
486       PRINT * , '--> FacMin/FacMax/FacRej must be positive'
487     CASE (-5) 
488       PRINT * , '--> Improper tolerance values'
489     CASE (-6) 
490       PRINT * , '--> No of steps exceeds maximum bound'
491     CASE (-7) 
492       PRINT * , '--> Step size too small: T + 10*H = T', &
493             ' or H < Roundoff'
494     CASE (-8)    
495       PRINT * , '--> Matrix is repeatedly singular'
496     CASE DEFAULT
497       PRINT *, 'Unknown Error code: ', Code
498    END SELECT
499    
500    PRINT *, "T=", T, "T + 10*H=", T+10.*H,"and H=", H
501      
502  END SUBROUTINE ros_ErrorMsg
503
504 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
505  SUBROUTINE ros_Integrator (Y, Tstart, Tend, T,  &
506         KMI,KVECNPT,                             &
507         AbsTol, RelTol,                          &
508 !~~~> Integration parameters
509         Autonomous, VectorTol, Max_no_steps,     &
510         Roundoff, Hmin, Hmax, Hstart,            &
511         FacMin, FacMax, FacRej, FacSafe,         &
512 !~~~> Error indicator
513         IERR )
514 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
515 !   Template for the implementation of a generic Rosenbrock method
516 !      defined by ros_S (no of stages)
517 !      and its coefficients ros_{A,C,M,E,Alpha,Gamma}
518 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
519
520   IMPLICIT NONE
521
522 !~~~> Input: the initial condition at Tstart; Output: the solution at T
523    REAL, INTENT(INOUT) :: Y(N)
524 !~~~> Input: integration interval
525    REAL, INTENT(IN) :: Tstart,Tend
526 !~~~> Output: time at which the solution is returned (T=Tend if success)
527    REAL, INTENT(OUT) ::  T
528 !
529 INTEGER, INTENT(IN) :: KMI      ! model number
530 INTEGER, INTENT(IN) :: KVECNPT
531 !
532 !~~~> Input: tolerances
533    REAL, INTENT(IN) ::  AbsTol(N), RelTol(N)
534 !~~~> Input: integration parameters
535    LOGICAL, INTENT(IN) :: Autonomous, VectorTol
536    REAL, INTENT(IN) :: Hstart, Hmin, Hmax
537    INTEGER, INTENT(IN) :: Max_no_steps
538    REAL, INTENT(IN) :: Roundoff, FacMin, FacMax, FacRej, FacSafe
539 !~~~> Output: Error indicator
540    INTEGER, INTENT(OUT) :: IERR
541 ! ~~~~ Local variables
542    REAL :: Ynew(N), Fcn0(N), Fcn(N)
543    REAL :: K(N*ros_S), dFdT(N)
544 !#ifdef FULL_ALGEBRA    
545 !##   REAL :: Jac0(N,N), Ghimj(N,N)
546 !#else
547    REAL :: Jac0(LU_NONZERO), Ghimj(LU_NONZERO)
548 !#endif
549    REAL :: H, Hnew, HC, HG, Fac, Tau
550    REAL :: Err, Yerr(N)
551    INTEGER :: Pivot(N), Direction, ioffset, j, istage
552    LOGICAL :: RejectLastH, RejectMoreH, Singular
553 !~~~>  Local parameters
554    REAL, PARAMETER :: ZERO = 0.0, ONE  = 1.0
555    REAL, PARAMETER :: DeltaMin = 1.0E-5
556 !~~~>  Locally called functions
557 !    REAL WLAMCH
558 !    EXTERNAL WLAMCH
559 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
560
561
562 !~~~>  Initial preparations
563    T = Tstart
564    RSTATUS(Nhexit) = ZERO
565    H = MIN( MAX(ABS(Hmin),ABS(Hstart)) , ABS(Hmax) )
566    IF (ABS(H) <= 10.0*Roundoff) H = DeltaMin
567
568    IF (Tend  >=  Tstart) THEN
569      Direction = +1
570    ELSE
571      Direction = -1
572    END IF
573    H = Direction*H
574
575    RejectLastH=.FALSE.
576    RejectMoreH=.FALSE.
577
578 !~~~> Time loop begins below
579
580 TimeLoop: DO WHILE ( (Direction > 0).AND.((T-Tend)+Roundoff <= ZERO) &
581        .OR. (Direction < 0).AND.((Tend-T)+Roundoff <= ZERO) )
582
583    IF ( ISTATUS(Nstp) > Max_no_steps ) THEN  ! Too many steps
584       CALL ros_ErrorMsg(-6,T,H,IERR)
585       RETURN
586    END IF
587    IF ( ((T+0.1*H) == T).OR.(H <= Roundoff) ) THEN  ! Step size too small
588       CALL ros_ErrorMsg(-7,T,H,IERR)
589       RETURN
590    END IF
591
592 !~~~>  Limit H if necessary to avoid going beyond Tend
593    H = MIN(H,ABS(Tend-T))
594
595 !~~~>   Compute the function at current time
596    CALL FunTemplate(N,T,Y,Fcn0,KMI,KVECNPT)
597    ISTATUS(Nfun) = ISTATUS(Nfun) + 1
598
599 !~~~>  Compute the function derivative with respect to T
600    IF (.NOT.Autonomous) THEN
601       CALL ros_FunTimeDerivative ( T, Roundoff, Y, &
602                 Fcn0, dFdT,KMI,KVECNPT )
603    END IF
604
605 !~~~>   Compute the Jacobian at current time
606    CALL JacTemplate(N,T,Y,Jac0,KMI,KVECNPT)
607    ISTATUS(Njac) = ISTATUS(Njac) + 1
608
609 !~~~>  Repeat step calculation until current step accepted
610 UntilAccepted: DO
611
612    CALL ros_PrepareMatrix(H,Direction,ros_Gamma(1), &
613           Jac0,Ghimj,Pivot,Singular)
614    IF (Singular) THEN ! More than 5 consecutive failed decompositions
615        CALL ros_ErrorMsg(-8,T,H,IERR)
616        RETURN
617    END IF
618
619 !~~~>   Compute the stages
620 Stage: DO istage = 1, ros_S
621
622       ! Current istage offset. Current istage vector is K(ioffset+1:ioffset+N)
623        ioffset = N*(istage-1)
624
625       ! For the 1st istage the function has been computed previously
626        IF ( istage == 1 ) THEN
627          CALL WCOPY(N,Fcn0,1,Fcn,1)
628       ! istage>1 and a new function evaluation is needed at the current istage
629        ELSEIF ( ros_NewF(istage) ) THEN
630          CALL WCOPY(N,Y,1,Ynew,1)
631          DO j = 1, istage-1
632            CALL WAXPY(N,ros_A((istage-1)*(istage-2)/2+j), &
633             K(N*(j-1)+1),1,Ynew,1)
634          END DO
635          Tau = T + ros_Alpha(istage)*Direction*H
636          CALL FunTemplate(N,Tau,Ynew,Fcn,KMI,KVECNPT)
637          ISTATUS(Nfun) = ISTATUS(Nfun) + 1
638        END IF ! if istage == 1 elseif ros_NewF(istage)
639        CALL WCOPY(N,Fcn,1,K(ioffset+1),1)
640        DO j = 1, istage-1
641          HC = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H)
642          CALL WAXPY(N,HC,K(N*(j-1)+1),1,K(ioffset+1),1)
643        END DO
644        IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN
645          HG = Direction*H*ros_Gamma(istage)
646          CALL WAXPY(N,HG,dFdT,1,K(ioffset+1),1)
647        END IF
648        CALL ros_Solve(Ghimj, Pivot, K(ioffset+1))
649
650    END DO Stage
651
652
653 !~~~>  Compute the new solution
654    CALL WCOPY(N,Y,1,Ynew,1)
655    DO j=1,ros_S
656          CALL WAXPY(N,ros_M(j),K(N*(j-1)+1),1,Ynew,1)
657    END DO
658
659 !~~~>  Compute the error estimation
660    CALL WSCAL(N,ZERO,Yerr,1)
661    DO j=1,ros_S
662         CALL WAXPY(N,ros_E(j),K(N*(j-1)+1),1,Yerr,1)
663    END DO
664    Err = ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol )
665
666 !~~~> New step size is bounded by FacMin <= Hnew/H <= FacMax
667    Fac  = MIN(FacMax,MAX(FacMin,FacSafe/Err**(ONE/ros_ELO)))
668    Hnew = H*Fac
669
670 !~~~>  Check the error magnitude and adjust step size
671    ISTATUS(Nstp) = ISTATUS(Nstp) + 1
672    IF ( (Err <= ONE).OR.(H <= Hmin) ) THEN  !~~~> Accept step
673       ISTATUS(Nacc) = ISTATUS(Nacc) + 1
674       CALL WCOPY(N,Ynew,1,Y,1)
675       T = T + Direction*H
676       Hnew = MAX(Hmin,MIN(Hnew,Hmax))
677       IF (RejectLastH) THEN  ! No step size increase after a rejected step
678          Hnew = MIN(Hnew,H)
679       END IF
680       RSTATUS(Nhexit) = H
681       RSTATUS(Nhnew)  = Hnew
682       RSTATUS(Ntexit) = T
683       RejectLastH = .FALSE.
684       RejectMoreH = .FALSE.
685       H = Hnew
686       EXIT UntilAccepted ! EXIT THE LOOP: WHILE STEP NOT ACCEPTED
687    ELSE           !~~~> Reject step
688       IF (RejectMoreH) THEN
689          Hnew = H*FacRej
690       END IF
691       RejectMoreH = RejectLastH
692       RejectLastH = .TRUE.
693       H = Hnew
694       IF (ISTATUS(Nacc) >= 1)  ISTATUS(Nrej) = ISTATUS(Nrej) + 1
695    END IF ! Err <= 1
696
697    END DO UntilAccepted
698
699    END DO TimeLoop
700
701 !~~~> Succesful exit
702    IERR = 1  !~~~> The integration was successful
703
704   END SUBROUTINE ros_Integrator
705
706
707 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
708   REAL FUNCTION ros_ErrorNorm ( Y, Ynew, Yerr, &
709                                AbsTol, RelTol, VectorTol )
710 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
711 !~~~> Computes the "scaled norm" of the error vector Yerr
712 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
713    IMPLICIT NONE
714
715 ! Input arguments
716    REAL, INTENT(IN) :: Y(N), Ynew(N), &
717           Yerr(N), AbsTol(N), RelTol(N)
718    LOGICAL, INTENT(IN) ::  VectorTol
719 ! Local variables
720    REAL :: Err, Scale, Ymax
721    INTEGER  :: i
722    REAL, PARAMETER :: ZERO = 0.0
723
724    Err = ZERO
725    DO i=1,N
726      Ymax = MAX(ABS(Y(i)),ABS(Ynew(i)))
727      IF (VectorTol) THEN
728        Scale = AbsTol(i)+RelTol(i)*Ymax
729      ELSE
730        Scale = AbsTol(1)+RelTol(1)*Ymax
731      END IF
732      Err = Err+(Yerr(i)/Scale)**2
733    END DO
734    Err  = SQRT(Err/N)
735
736    ros_ErrorNorm = MAX(Err,1.0d-10)
737
738   END FUNCTION ros_ErrorNorm
739
740
741 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
742   SUBROUTINE ros_FunTimeDerivative ( T, Roundoff, Y, &
743                 Fcn0, dFdT,KMI,KVECNPT )
744 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
745 !~~~> The time partial derivative of the function by finite differences
746 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
747    IMPLICIT NONE
748
749 !~~~> Input arguments
750    REAL, INTENT(IN) :: T, Roundoff, Y(N), Fcn0(N)
751 !~~~> Output arguments
752    REAL, INTENT(OUT) :: dFdT(N)
753 !~~~> Local variables
754    REAL :: Delta
755    REAL, PARAMETER :: ONE = 1.0, DeltaMin = 1.0E-6
756 !
757 INTEGER, INTENT(IN) :: KMI      ! model number
758 INTEGER, INTENT(IN) :: KVECNPT
759 !
760
761
762    Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T))
763    CALL FunTemplate(N,T+Delta,Y,dFdT,KMI,KVECNPT)
764    ISTATUS(Nfun) = ISTATUS(Nfun) + 1
765    CALL WAXPY(N,(-ONE),Fcn0,1,dFdT,1)
766    CALL WSCAL(N,(ONE/Delta),dFdT,1)
767
768   END SUBROUTINE ros_FunTimeDerivative
769
770
771 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
772   SUBROUTINE ros_PrepareMatrix ( H, Direction, gam, &
773              Jac0, Ghimj, Pivot, Singular )
774 ! --- --- --- --- --- --- --- --- --- --- --- --- ---
775 !  Prepares the LHS matrix for stage calculations
776 !  1.  Construct Ghimj = 1/(H*ham) - Jac0
777 !      "(Gamma H) Inverse Minus Jacobian"
778 !  2.  Repeat LU decomposition of Ghimj until successful.
779 !       -half the step size if LU decomposition fails and retry
780 !       -exit after 5 consecutive fails
781 ! --- --- --- --- --- --- --- --- --- --- --- --- ---
782    IMPLICIT NONE
783
784 !~~~> Input arguments
785 !#ifdef FULL_ALGEBRA    
786 !##   REAL, INTENT(IN) ::  Jac0(N,N)
787 !#else
788    REAL, INTENT(IN) ::  Jac0(LU_NONZERO)
789 !#endif   
790    REAL, INTENT(IN) ::  gam
791    INTEGER, INTENT(IN) ::  Direction
792 !~~~> Output arguments
793 !#ifdef FULL_ALGEBRA    
794 !##   REAL, INTENT(OUT) :: Ghimj(N,N)
795 !#else
796    REAL, INTENT(OUT) :: Ghimj(LU_NONZERO)
797 !#endif   
798    LOGICAL, INTENT(OUT) ::  Singular
799    INTEGER, INTENT(OUT) ::  Pivot(N)
800 !~~~> Inout arguments
801    REAL, INTENT(INOUT) :: H   ! step size is decreased when LU fails
802 !~~~> Local variables
803    INTEGER  :: i, ising, Nconsecutive
804    REAL :: ghinv
805    REAL, PARAMETER :: ONE  = 1.0, HALF = 0.5
806
807    Nconsecutive = 0
808    Singular = .TRUE.
809
810    DO WHILE (Singular)
811
812 !~~~>    Construct Ghimj = 1/(H*gam) - Jac0
813 !#ifdef FULL_ALGEBRA    
814 !##     CALL WCOPY(N*N,Jac0,1,Ghimj,1)
815 !##     CALL WSCAL(N*N,(-ONE),Ghimj,1)
816 !##     ghinv = ONE/(Direction*H*gam)
817 !##     DO i=1,N
818 !##       Ghimj(i,i) = Ghimj(i,i)+ghinv
819 !##     END DO
820 !#else
821      CALL WCOPY(LU_NONZERO,Jac0,1,Ghimj,1)
822      CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1)
823      ghinv = ONE/(Direction*H*gam)
824      DO i=1,N
825        Ghimj(LU_DIAG(i)) = Ghimj(LU_DIAG(i))+ghinv
826      END DO
827 !#endif   
828 !~~~>    Compute LU decomposition
829      CALL ros_Decomp( Ghimj, Pivot, ising )
830      IF (ising == 0) THEN
831 !~~~>    If successful done
832         Singular = .FALSE.
833      ELSE ! ising .ne. 0
834 !~~~>    If unsuccessful half the step size; if 5 consecutive fails then return
835         ISTATUS(Nsng) = ISTATUS(Nsng) + 1
836         Nconsecutive = Nconsecutive+1
837         Singular = .TRUE.
838         PRINT*,'Warning: LU Decomposition returned ising = ',ising
839         IF (Nconsecutive <= 5) THEN ! Less than 5 consecutive failed decompositions
840            H = H*HALF
841         ELSE  ! More than 5 consecutive failed decompositions
842            RETURN
843         END IF  ! Nconsecutive
844       END IF    ! ising
845
846    END DO ! WHILE Singular
847
848   END SUBROUTINE ros_PrepareMatrix
849
850
851 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
852   SUBROUTINE ros_Decomp( A, Pivot, ising )
853 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
854 !  Template for the LU decomposition
855 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
856    IMPLICIT NONE
857 !~~~> Inout variables
858    REAL, INTENT(INOUT) :: A(LU_NONZERO)
859 !~~~> Output variables
860    INTEGER, INTENT(OUT) :: Pivot(N), ising
861
862 !#ifdef FULL_ALGEBRA    
863 !##   CALL  DGETRF( N, N, A, N, Pivot, ising )
864 !#else   
865    CALL KppDecomp ( A, ising )
866    Pivot(1) = 1
867 !#endif
868    ISTATUS(Ndec) = ISTATUS(Ndec) + 1
869
870   END SUBROUTINE ros_Decomp
871
872
873 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
874   SUBROUTINE ros_Solve( A, Pivot, b )
875 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
876 !  Template for the forward/backward substitution (using pre-computed LU decomposition)
877 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
878    IMPLICIT NONE
879 !~~~> Input variables
880    REAL, INTENT(IN) :: A(LU_NONZERO)
881    INTEGER, INTENT(IN) :: Pivot(N)
882 !~~~> InOut variables
883    REAL, INTENT(INOUT) :: b(N)
884
885 !#ifdef FULL_ALGEBRA    
886 !##   CALL  DGETRS( 'N', N , 1, A, N, Pivot, b, N, 0 )
887 !#else   
888    CALL KppSolveIndirect( A, b )
889 !#endif
890
891    ISTATUS(Nsol) = ISTATUS(Nsol) + 1
892
893   END SUBROUTINE ros_Solve
894
895
896   SUBROUTINE Ros1
897 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
898 ! --- AN L-STABLE METHOD, 1 stage , order 1
899 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
900
901    IMPLICIT NONE
902    REAL g
903
904     g = 1.0
905     rosMethod = RS1
906 !~~~> Name of the method
907     ros_Name = 'ROS-1'
908 !~~~> Number of stages
909     ros_S = 1
910
911 !~~~> The coefficient matrices A and C are strictly lower triangular.
912 !   The lower triangular (subdiagonal) elements are stored in row-wise order:
913 !   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
914 !   The general mapping formula is:
915 !       A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
916 !       C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
917
918     ros_A(1) = 0.0
919     ros_C(1) = 0.0
920 !~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
921 !   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
922     ros_NewF(1) = .TRUE.
923 !~~~> M_i = Coefficients for new step solution
924     ros_M(1) = 1.0
925 ! E_i = Coefficients for error estimator
926     ros_E(1) = 0.0
927 !~~~> ros_ELO = estimator of local order - the minimum between the
928 !    main and the embedded scheme orders plus one
929     ros_ELO = 1.0
930 !~~~> Y_stage_i ~ Y( T + H*Alpha_i )
931     ros_Alpha(1) = 0.0
932 !~~~> Gamma_i = \sum_j  gamma_{i,j}
933     ros_Gamma(1) = g
934
935  END SUBROUTINE Ros1
936
937
938 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
939   SUBROUTINE Ros2
940 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
941 ! --- AN L-STABLE METHOD, 2 stages, order 2
942 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
943
944    IMPLICIT NONE
945    REAL g
946
947     g = 1.0 + 1.0/SQRT(2.0)
948     rosMethod = RS2
949 !~~~> Name of the method
950     ros_Name = 'ROS-2'
951 !~~~> Number of stages
952     ros_S = 2
953
954 !~~~> The coefficient matrices A and C are strictly lower triangular.
955 !   The lower triangular (subdiagonal) elements are stored in row-wise order:
956 !   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
957 !   The general mapping formula is:
958 !       A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
959 !       C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
960
961     ros_A(1) = (1.0)/g
962     ros_C(1) = (-2.0)/g
963 !~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
964 !   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
965     ros_NewF(1) = .TRUE.
966     ros_NewF(2) = .TRUE.
967 !~~~> M_i = Coefficients for new step solution
968     ros_M(1)= (3.0)/(2.0*g)
969     ros_M(2)= (1.0)/(2.0*g)
970 ! E_i = Coefficients for error estimator
971     ros_E(1) = 1.0/(2.0*g)
972     ros_E(2) = 1.0/(2.0*g)
973 !~~~> ros_ELO = estimator of local order - the minimum between the
974 !    main and the embedded scheme orders plus one
975     ros_ELO = 2.0
976 !~~~> Y_stage_i ~ Y( T + H*Alpha_i )
977     ros_Alpha(1) = 0.0
978     ros_Alpha(2) = 1.0
979 !~~~> Gamma_i = \sum_j  gamma_{i,j}
980     ros_Gamma(1) = g
981     ros_Gamma(2) =-g
982
983  END SUBROUTINE Ros2
984
985
986 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
987   SUBROUTINE Ros3
988 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
989 ! --- AN L-STABLE METHOD, 3 stages, order 3, 2 function evaluations
990 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
991
992    IMPLICIT NONE
993    rosMethod = RS3
994 !~~~> Name of the method
995    ros_Name = 'ROS-3'
996 !~~~> Number of stages
997    ros_S = 3
998
999 !~~~> The coefficient matrices A and C are strictly lower triangular.
1000 !   The lower triangular (subdiagonal) elements are stored in row-wise order:
1001 !   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
1002 !   The general mapping formula is:
1003 !       A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
1004 !       C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
1005
1006    ros_A(1)= 1.0
1007    ros_A(2)= 1.0
1008    ros_A(3)= 0.0
1009
1010    ros_C(1) = -0.10156171083877702091975600115545E+01
1011    ros_C(2) =  0.40759956452537699824805835358067E+01
1012    ros_C(3) =  0.92076794298330791242156818474003E+01
1013 !~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
1014 !   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
1015    ros_NewF(1) = .TRUE.
1016    ros_NewF(2) = .TRUE.
1017    ros_NewF(3) = .FALSE.
1018 !~~~> M_i = Coefficients for new step solution
1019    ros_M(1) =  0.1E+01
1020    ros_M(2) =  0.61697947043828245592553615689730E+01
1021    ros_M(3) = -0.42772256543218573326238373806514
1022 ! E_i = Coefficients for error estimator
1023    ros_E(1) =  0.5
1024    ros_E(2) = -0.29079558716805469821718236208017E+01
1025    ros_E(3) =  0.22354069897811569627360909276199
1026 !~~~> ros_ELO = estimator of local order - the minimum between the
1027 !    main and the embedded scheme orders plus 1
1028    ros_ELO = 3.0
1029 !~~~> Y_stage_i ~ Y( T + H*Alpha_i )
1030    ros_Alpha(1)= 0.0
1031    ros_Alpha(2)= 0.43586652150845899941601945119356
1032    ros_Alpha(3)= 0.43586652150845899941601945119356
1033 !~~~> Gamma_i = \sum_j  gamma_{i,j}
1034    ros_Gamma(1)= 0.43586652150845899941601945119356
1035    ros_Gamma(2)= 0.24291996454816804366592249683314
1036    ros_Gamma(3)= 0.21851380027664058511513169485832E+01
1037
1038   END SUBROUTINE Ros3
1039
1040 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1041
1042
1043 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1044   SUBROUTINE Ros4
1045 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1046 !     L-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 4 STAGES
1047 !     L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3
1048 !
1049 !      E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
1050 !      EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
1051 !      SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
1052 !      SPRINGER-VERLAG (1990)
1053 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1054
1055    IMPLICIT NONE
1056
1057    rosMethod = RS4
1058 !~~~> Name of the method
1059    ros_Name = 'ROS-4'
1060 !~~~> Number of stages
1061    ros_S = 4
1062
1063 !~~~> The coefficient matrices A and C are strictly lower triangular.
1064 !   The lower triangular (subdiagonal) elements are stored in row-wise order:
1065 !   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
1066 !   The general mapping formula is:
1067 !       A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
1068 !       C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
1069
1070    ros_A(1) = 0.2000000000000000E+01
1071    ros_A(2) = 0.1867943637803922E+01
1072    ros_A(3) = 0.2344449711399156
1073    ros_A(4) = ros_A(2)
1074    ros_A(5) = ros_A(3)
1075    ros_A(6) = 0.0
1076
1077    ros_C(1) =-0.7137615036412310E+01
1078    ros_C(2) = 0.2580708087951457E+01
1079    ros_C(3) = 0.6515950076447975
1080    ros_C(4) =-0.2137148994382534E+01
1081    ros_C(5) =-0.3214669691237626
1082    ros_C(6) =-0.6949742501781779
1083 !~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
1084 !   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
1085    ros_NewF(1)  = .TRUE.
1086    ros_NewF(2)  = .TRUE.
1087    ros_NewF(3)  = .TRUE.
1088    ros_NewF(4)  = .FALSE.
1089 !~~~> M_i = Coefficients for new step solution
1090    ros_M(1) = 0.2255570073418735E+01
1091    ros_M(2) = 0.2870493262186792
1092    ros_M(3) = 0.4353179431840180
1093    ros_M(4) = 0.1093502252409163E+01
1094 !~~~> E_i  = Coefficients for error estimator
1095    ros_E(1) =-0.2815431932141155
1096    ros_E(2) =-0.7276199124938920E-01
1097    ros_E(3) =-0.1082196201495311
1098    ros_E(4) =-0.1093502252409163E+01
1099 !~~~> ros_ELO  = estimator of local order - the minimum between the
1100 !    main and the embedded scheme orders plus 1
1101    ros_ELO  = 4.0
1102 !~~~> Y_stage_i ~ Y( T + H*Alpha_i )
1103    ros_Alpha(1) = 0.0
1104    ros_Alpha(2) = 0.1145640000000000E+01
1105    ros_Alpha(3) = 0.6552168638155900
1106    ros_Alpha(4) = ros_Alpha(3)
1107 !~~~> Gamma_i = \sum_j  gamma_{i,j}
1108    ros_Gamma(1) = 0.5728200000000000
1109    ros_Gamma(2) =-0.1769193891319233E+01
1110    ros_Gamma(3) = 0.7592633437920482
1111    ros_Gamma(4) =-0.1049021087100450
1112
1113   END SUBROUTINE Ros4
1114
1115 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1116   SUBROUTINE Rodas3
1117 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1118 ! --- A STIFFLY-STABLE METHOD, 4 stages, order 3
1119 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1120
1121    IMPLICIT NONE
1122
1123    rosMethod = RD3
1124 !~~~> Name of the method
1125    ros_Name = 'RODAS-3'
1126 !~~~> Number of stages
1127    ros_S = 4
1128
1129 !~~~> The coefficient matrices A and C are strictly lower triangular.
1130 !   The lower triangular (subdiagonal) elements are stored in row-wise order:
1131 !   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
1132 !   The general mapping formula is:
1133 !       A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
1134 !       C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
1135
1136    ros_A(1) = 0.0
1137    ros_A(2) = 2.0
1138    ros_A(3) = 0.0
1139    ros_A(4) = 2.0
1140    ros_A(5) = 0.0
1141    ros_A(6) = 1.0
1142
1143    ros_C(1) = 4.0
1144    ros_C(2) = 1.0
1145    ros_C(3) =-1.0
1146    ros_C(4) = 1.0
1147    ros_C(5) =-1.0
1148    ros_C(6) =-(8.0/3.0)
1149
1150 !~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
1151 !   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
1152    ros_NewF(1)  = .TRUE.
1153    ros_NewF(2)  = .FALSE.
1154    ros_NewF(3)  = .TRUE.
1155    ros_NewF(4)  = .TRUE.
1156 !~~~> M_i = Coefficients for new step solution
1157    ros_M(1) = 2.0
1158    ros_M(2) = 0.0
1159    ros_M(3) = 1.0
1160    ros_M(4) = 1.0
1161 !~~~> E_i  = Coefficients for error estimator
1162    ros_E(1) = 0.0
1163    ros_E(2) = 0.0
1164    ros_E(3) = 0.0
1165    ros_E(4) = 1.0
1166 !~~~> ros_ELO  = estimator of local order - the minimum between the
1167 !    main and the embedded scheme orders plus 1
1168    ros_ELO  = 3.0
1169 !~~~> Y_stage_i ~ Y( T + H*Alpha_i )
1170    ros_Alpha(1) = 0.0
1171    ros_Alpha(2) = 0.0
1172    ros_Alpha(3) = 1.0
1173    ros_Alpha(4) = 1.0
1174 !~~~> Gamma_i = \sum_j  gamma_{i,j}
1175    ros_Gamma(1) = 0.5
1176    ros_Gamma(2) = 1.5
1177    ros_Gamma(3) = 0.0
1178    ros_Gamma(4) = 0.0
1179
1180   END SUBROUTINE Rodas3
1181
1182 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1183   SUBROUTINE Rodas4
1184 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1185 !     STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 6 STAGES
1186 !
1187 !      E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
1188 !      EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
1189 !      SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
1190 !      SPRINGER-VERLAG (1996)
1191 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1192
1193    IMPLICIT NONE
1194
1195     rosMethod = RD4
1196 !~~~> Name of the method
1197     ros_Name = 'RODAS-4'
1198 !~~~> Number of stages
1199     ros_S = 6
1200
1201 !~~~> Y_stage_i ~ Y( T + H*Alpha_i )
1202     ros_Alpha(1) = 0.000
1203     ros_Alpha(2) = 0.386
1204     ros_Alpha(3) = 0.210
1205     ros_Alpha(4) = 0.630
1206     ros_Alpha(5) = 1.000
1207     ros_Alpha(6) = 1.000
1208
1209 !~~~> Gamma_i = \sum_j  gamma_{i,j}
1210     ros_Gamma(1) = 0.2500000000000000
1211     ros_Gamma(2) =-0.1043000000000000
1212     ros_Gamma(3) = 0.1035000000000000
1213     ros_Gamma(4) =-0.3620000000000023E-01
1214     ros_Gamma(5) = 0.0
1215     ros_Gamma(6) = 0.0
1216
1217 !~~~> The coefficient matrices A and C are strictly lower triangular.
1218 !   The lower triangular (subdiagonal) elements are stored in row-wise order:
1219 !   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
1220 !   The general mapping formula is:  A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
1221 !                  C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
1222
1223     ros_A(1) = 0.1544000000000000E+01
1224     ros_A(2) = 0.9466785280815826
1225     ros_A(3) = 0.2557011698983284
1226     ros_A(4) = 0.3314825187068521E+01
1227     ros_A(5) = 0.2896124015972201E+01
1228     ros_A(6) = 0.9986419139977817
1229     ros_A(7) = 0.1221224509226641E+01
1230     ros_A(8) = 0.6019134481288629E+01
1231     ros_A(9) = 0.1253708332932087E+02
1232     ros_A(10) =-0.6878860361058950
1233     ros_A(11) = ros_A(7)
1234     ros_A(12) = ros_A(8)
1235     ros_A(13) = ros_A(9)
1236     ros_A(14) = ros_A(10)
1237     ros_A(15) = 1.0
1238
1239     ros_C(1) =-0.5668800000000000E+01
1240     ros_C(2) =-0.2430093356833875E+01
1241     ros_C(3) =-0.2063599157091915
1242     ros_C(4) =-0.1073529058151375
1243     ros_C(5) =-0.9594562251023355E+01
1244     ros_C(6) =-0.2047028614809616E+02
1245     ros_C(7) = 0.7496443313967647E+01
1246     ros_C(8) =-0.1024680431464352E+02
1247     ros_C(9) =-0.3399990352819905E+02
1248     ros_C(10) = 0.1170890893206160E+02
1249     ros_C(11) = 0.8083246795921522E+01
1250     ros_C(12) =-0.7981132988064893E+01
1251     ros_C(13) =-0.3152159432874371E+02
1252     ros_C(14) = 0.1631930543123136E+02
1253     ros_C(15) =-0.6058818238834054E+01
1254
1255 !~~~> M_i = Coefficients for new step solution
1256     ros_M(1) = ros_A(7)
1257     ros_M(2) = ros_A(8)
1258     ros_M(3) = ros_A(9)
1259     ros_M(4) = ros_A(10)
1260     ros_M(5) = 1.0
1261     ros_M(6) = 1.0
1262
1263 !~~~> E_i  = Coefficients for error estimator
1264     ros_E(1) = 0.0
1265     ros_E(2) = 0.0
1266     ros_E(3) = 0.0
1267     ros_E(4) = 0.0
1268     ros_E(5) = 0.0
1269     ros_E(6) = 1.0
1270
1271 !~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
1272 !   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
1273     ros_NewF(1) = .TRUE.
1274     ros_NewF(2) = .TRUE.
1275     ros_NewF(3) = .TRUE.
1276     ros_NewF(4) = .TRUE.
1277     ros_NewF(5) = .TRUE.
1278     ros_NewF(6) = .TRUE.
1279
1280 !~~~> ros_ELO  = estimator of local order - the minimum between the
1281 !        main and the embedded scheme orders plus 1
1282     ros_ELO = 4.0
1283
1284   END SUBROUTINE Rodas4
1285
1286 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1287 !   End of the set of internal Rosenbrock subroutines
1288 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1289 END SUBROUTINE Rosenbrock
1290 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1291
1292
1293 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1294 SUBROUTINE FunTemplate( N, T, Y, Ydot, KMI, KVECNPT )
1295 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1296 !  Template for the ODE function call.
1297 !  Updates the rate coefficients (and possibly the fixed species) at each call
1298 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1299  USE MODD_CH_M9_n, ONLY: NEQ
1300  USE MODD_RBK90_JacobianSP_n, ONLY: LU_DIM_SPECIES
1301
1302 !JPP USE MODD_RBK90_Parameters_n, ONLY: NVAR
1303 !JPP USE MODD_RBK90_Global_n, ONLY: FIX, RCONST, TIME
1304 !JPP USE RBK90_Function, ONLY: Fun
1305 !JPP USE RBK90_Rates, ONLY: Update_SUN, Update_RCONST
1306
1307 USE MODI_CH_FCN
1308
1309 !~~~> Input variables
1310     INTEGER           :: N
1311     REAL :: T, Y(N)
1312 !~~~> Output variables
1313     REAL :: Ydot(N)
1314 !~~~> Local variables
1315     REAL :: Told
1316     REAL :: TIME
1317 !
1318 INTEGER, INTENT(IN) :: KMI      ! model number
1319 INTEGER, INTENT(IN) :: KVECNPT
1320 !
1321 !~~~> Local variables
1322     REAL :: ZCONC(KVECNPT,NEQ)
1323     REAL ::  ZFCN(KVECNPT,NEQ)
1324     INTEGER       :: JL, JLSHIFT   ! Loop indexes
1325     INTEGER       :: ISPECIES      ! Ancillary variable
1326
1327
1328     Told = TIME
1329     TIME = T
1330 !JPP   CALL Update_SUN()
1331 !JPP   CALL Update_RCONST()
1332 !JPP   CALL Fun( Y, FIX, RCONST, Ydot )
1333    
1334     JLSHIFT = 0
1335     ZCONC(:,:) = 0.0
1336     DO JL = 1,KVECNPT
1337       ISPECIES = LU_DIM_SPECIES(JL)
1338       ZCONC(JL,1:ISPECIES) = Y(JLSHIFT+1:JLSHIFT+ISPECIES)
1339       JLSHIFT = JLSHIFT + ISPECIES
1340     END DO
1341
1342     CALL CH_FCN( T, ZCONC, ZFCN, KMI, KVECNPT, NEQ)
1343
1344     JLSHIFT = 0
1345     DO JL = 1,KVECNPT
1346       ISPECIES = LU_DIM_SPECIES(JL)
1347       Ydot(JLSHIFT+1:JLSHIFT+ISPECIES) = ZFCN(JL,1:ISPECIES)
1348       JLSHIFT = JLSHIFT + ISPECIES
1349     END DO
1350     TIME = Told
1351
1352 END SUBROUTINE FunTemplate
1353
1354
1355 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1356 SUBROUTINE JacTemplate( N, T, Y, Jcb, KMI, KVECNPT )
1357 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1358 !  Template for the ODE Jacobian call.
1359 !  Updates the rate coefficients (and possibly the fixed species) at each call
1360 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1361  USE MODD_CH_M9_n,     ONLY: NEQ
1362  USE MODD_CH_ROSENBROCK_n,    ONLY: NSPARSEDIM, NSPARSE_ICOL, NSPARSE_IROW, &
1363                           NSPARSEDIM_NAQ, NSPARSE_ICOL_NAQ, NSPARSE_IROW_NAQ
1364  USE MODD_RBK90_Parameters_n, ONLY: LU_NONZERO
1365  USE MODD_RBK90_JacobianSP_n, ONLY: LU_IROW, LU_ICOL, LU_DIM_SPECIES
1366
1367 !JPP USE MODD_RBK90_Global_n, ONLY: FIX, RCONST, TIME
1368 !JPP USE RBK90_LinearAlgebra
1369 !JPP USE RBK90_Rates, ONLY: Update_SUN, Update_RCONST
1370
1371 USE MODI_CH_JAC
1372
1373 !~~~> Input variables
1374     INTEGER           :: N
1375     REAL :: T, Y(N)
1376 !~~~> Output variables
1377 !#ifdef FULL_ALGEBRA    
1378 !##    REAL :: JV(LU_NONZERO), Jcb(NVAR,NVAR)
1379 !#else
1380     REAL :: Jcb(LU_NONZERO)
1381 !#endif   
1382 !
1383 INTEGER, INTENT(IN) :: KMI      ! model number
1384 INTEGER, INTENT(IN) :: KVECNPT
1385 !
1386 !~~~> Local variables
1387     REAL :: Told
1388     REAL :: TIME
1389 !#ifdef FULL_ALGEBRA    
1390 !##    INTEGER :: i, j
1391 !#endif   
1392     REAL :: ZCONC(KVECNPT,NEQ)
1393     REAL ::  ZJAC(KVECNPT,NEQ,NEQ)
1394     INTEGER       :: JL, JLL, JLSHIFT   ! Loop indexes
1395     INTEGER       :: ISPECIES      ! Ancillary variable
1396
1397     Told = TIME
1398     TIME = T
1399 !JPP    CALL Update_SUN()
1400 !JPP    CALL Update_RCONST()
1401 !#ifdef FULL_ALGEBRA    
1402 !##    CALL Jac_SP(Y, FIX, RCONST, JV)
1403 !##    DO j=1,NVAR
1404 !##      DO i=1,NVAR
1405 !##         Jcb(i,j) = 0.0
1406 !##      END DO
1407 !##    END DO
1408 !##    DO i=1,LU_NONZERO
1409 !##       Jcb(LU_IROW(i),LU_ICOL(i)) = JV(i)
1410 !##    END DO
1411 !#else
1412 !JPP    CALL Jac_SP( Y, FIX, RCONST, Jcb )
1413 !#endif   
1414     JLSHIFT = 0
1415     ZCONC(:,:) = 0.0
1416     DO JL = 1,KVECNPT
1417       ISPECIES = LU_DIM_SPECIES(JL)
1418       ZCONC(JL,1:ISPECIES) = Y(JLSHIFT+1:JLSHIFT+ISPECIES)
1419       JLSHIFT = JLSHIFT + ISPECIES
1420     END DO
1421
1422     CALL CH_JAC( T, ZCONC, ZJAC, KMI, KVECNPT, NEQ)
1423
1424     JLSHIFT = 0
1425     DO JL = 1,KVECNPT
1426       ISPECIES = LU_DIM_SPECIES(JL)
1427       IF( ISPECIES==NEQ ) THEN
1428         DO JLL = 1, NSPARSEDIM
1429           Jcb(JLSHIFT+JLL) = ZJAC(JL,NSPARSE_IROW(JLL),NSPARSE_ICOL(JLL))
1430         END DO
1431         JLSHIFT = JLSHIFT + NSPARSEDIM
1432       ELSE
1433         DO JLL = 1, NSPARSEDIM_NAQ
1434           Jcb(JLSHIFT+JLL) =ZJAC(JL,NSPARSE_IROW_NAQ(JLL),NSPARSE_ICOL_NAQ(JLL))
1435         END DO
1436         JLSHIFT = JLSHIFT + NSPARSEDIM_NAQ
1437       END IF
1438     END DO
1439     TIME = Told
1440
1441 END SUBROUTINE JacTemplate
1442
1443 END MODULE MODE_RBK90_Integrator