source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int.modified_WCOPY/rosenbrock.f @ 4650

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

Merge of branch palm4u into trunk

File size: 45.3 KB
Line 
1      SUBROUTINE INTEGRATE( TIN, TOUT )
2
3      IMPLICIT NONE     
4      INCLUDE 'KPP_ROOT_Parameters.h'
5      INCLUDE 'KPP_ROOT_Global.h'
6      INTEGER Nstp, Nacc, Nrej, Nsng, IERR
7      SAVE Nstp, Nacc, Nrej, Nsng
8
9! TIN - Start Time
10      KPP_REAL TIN
11! TOUT - End Time
12      KPP_REAL TOUT
13      INTEGER i
14
15      KPP_REAL  RPAR(20)
16      INTEGER  IPAR(20)
17      EXTERNAL  FunTemplate, JacTemplate
18
19
20      DO i=1,20
21        IPAR(i) = 0
22        RPAR(i) = 0.0d0
23      ENDDO
24     
25     
26      IPAR(1) = 0       ! non-autonomous
27      IPAR(2) = 1       ! vector tolerances
28      RPAR(3) = STEPMIN ! starting step
29      IPAR(4) = 5       ! choice of the method
30
31      CALL Rosenbrock(VAR,TIN,TOUT,
32     &            ATOL,RTOL,
33     &            FunTemplate,JacTemplate,
34     &            RPAR,IPAR,IERR)
35
36               
37      Nstp = Nstp + IPAR(13)
38      Nacc = Nacc + IPAR(14)
39      Nrej = Nrej + IPAR(15)
40      Nsng = Nsng + IPAR(18)
41      PRINT*,'Step=',Nstp,' Acc=',Nacc,' Rej=',Nrej,
42     &      ' Singular=',Nsng
43
44
45      IF (IERR.LT.0) THEN
46        print *,'Rosenbrock: Unsucessful step at T=',
47     &          TIN,' (IERR=',IERR,')'
48      ENDIF
49     
50      TIN = RPAR(11) ! Exit time
51      STEPMIN = RPAR(12)
52     
53      RETURN
54      END
55
56
57!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
58      SUBROUTINE Rosenbrock(Y,Tstart,Tend,
59     &            AbsTol,RelTol,
60     &            ode_Fun,ode_Jac ,
61     &            RPAR,IPAR,IERR)
62!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
63!     
64!    Solves the system y'=F(t,y) using a Rosenbrock method defined by:
65!
66!     G = 1/(H*gamma(1)) - ode_Jac(t0,Y0)
67!     T_i = t0 + Alpha(i)*H
68!     Y_i = Y0 + \sum_{j=1}^{i-1} A(i,j)*K_j
69!     G * K_i = ode_Fun( T_i, Y_i ) + \sum_{j=1}^S C(i,j)/H * K_j +
70!                  gamma(i)*dF/dT(t0, Y0)
71!     Y1 = Y0 + \sum_{j=1}^S M(j)*K_j
72!
73!    For details on Rosenbrock methods and their implementation consult:
74!         E. Hairer and G. Wanner
75!         "Solving ODEs II. Stiff and differential-algebraic problems".
76!         Springer series in computational mathematics, Springer-Verlag, 1996. 
77!    The codes contained in the book inspired this implementation.             
78!
79!    (C)  Adrian Sandu, August 2004
80!    Virginia Polytechnic Institute and State University   
81!    Contact: sandu@cs.vt.edu
82!    This implementation is part of KPP - the Kinetic PreProcessor
83!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
84!       
85!~~~>      INPUT ARGUMENTS:
86!       
87!-     Y(NVAR)       = vector of initial conditions (at T=Tstart)
88!-    [Tstart,Tend]  = time range of integration
89!        (if Tstart>Tend the integration is performed backwards in time) 
90!-    RelTol, AbsTol = user precribed accuracy
91!-    SUBROUTINE ode_Fun( T, Y, Ydot ) = ODE function,
92!                                         returns Ydot = Y' = F(T,Y)
93!-    SUBROUTINE ode_Fun( T, Y, Ydot ) = Jacobian of the ODE function,
94!                                         returns Jcb = dF/dY
95!-    IPAR(1:10)    = integer inputs parameters
96!-    RPAR(1:10)    = real inputs parameters
97!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
98!
99!~~~>     OUTPUT ARGUMENTS: 
100!       
101!-    Y(NVAR)       -> vector of final states (at T->Tend)
102!-    IPAR(11:20)   -> integer output parameters
103!-    RPAR(11:20)   -> real output parameters
104!-    IERR          -> job status upon return
105!          - succes (positive value) or failure (negative value) -
106!                    =  1 : Success
107!                    = -1 : Improper value for maximal no of steps
108!                    = -2 : Selected Rosenbrock method not implemented
109!                    = -3 : Hmin/Hmax/Hstart must be positive
110!                    = -4 : FacMin/FacMax/FacRej must be positive
111!                    = -5 : Improper tolerance values
112!                    = -6 : No of steps exceeds maximum bound
113!                    = -7 : Step size too small
114!                    = -8 : Matrix is repeatedly singular
115!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
116!
117!~~~>     INPUT PARAMETERS:
118!
119!    Note: For input parameters equal to zero the default values of the
120!          corresponding variables are used.
121!
122!    IPAR(1)   = 1: F = F(y)   Independent of T (AUTONOMOUS)
123!              = 0: F = F(t,y) Depends on T (NON-AUTONOMOUS)
124!    IPAR(2)   = 0: AbsTol, RelTol are NVAR-dimensional vectors
125!              = 1:  AbsTol, RelTol are scalars
126!    IPAR(3)  -> maximum number of integration steps
127!              For IPAR(3)=0) the default value of 100000 is used
128!
129!    IPAR(4)  -> selection of a particular Rosenbrock method
130!              = 0 :  default method is Rodas3
131!              = 1 :  method is  Ros2
132!              = 2 :  method is  Ros3
133!              = 3 :  method is  Ros4
134!              = 4 :  method is  Rodas3
135!              = 5:   method is  Rodas4
136!
137!    RPAR(1)  -> Hmin, lower bound for the integration step size
138!                It is strongly recommended to keep Hmin = ZERO
139!    RPAR(2)  -> Hmax, upper bound for the integration step size
140!    RPAR(3)  -> Hstart, starting value for the integration step size
141!               
142!    RPAR(4)  -> FacMin, lower bound on step decrease factor (default=0.2)
143!    RPAR(5)  -> FacMin,upper bound on step increase factor (default=6)
144!    RPAR(6)  -> FacRej, step decrease factor after multiple rejections
145!                        (default=0.1)
146!    RPAR(7)  -> FacSafe, by which the new step is slightly smaller
147!                  than the predicted value  (default=0.9)
148!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
149!
150!~~~>     OUTPUT PARAMETERS:
151!
152!    Note: each call to Rosenbrock adds the corrent no. of fcn calls
153!         to previous value of IPAR(11), and similar for the other params.
154!         Set IPAR(11:20) = 0 before call to avoid this accumulation.
155!
156!    IPAR(11) = No. of function calls
157!    IPAR(12) = No. of jacobian calls
158!    IPAR(13) = No. of steps
159!    IPAR(14) = No. of accepted steps
160!    IPAR(15) = No. of rejected steps (except at the beginning)
161!    IPAR(16) = No. of LU decompositions
162!    IPAR(17) = No. of forward/backward substitutions
163!    IPAR(18) = No. of singular matrix decompositions
164!
165!    RPAR(11)  -> Texit, the time corresponding to the
166!                        computed Y upon return
167!    RPAR(12)  -> Hexit, last accepted step before exit
168!    For multiple restarts, use Hexit as Hstart in the following run
169!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
170
171      IMPLICIT NONE
172      INCLUDE 'KPP_ROOT_Parameters.h'
173      INCLUDE 'KPP_ROOT_Sparse.h'
174     
175      KPP_REAL Tstart,Tend
176      KPP_REAL Y(KPP_NVAR),AbsTol(KPP_NVAR),RelTol(KPP_NVAR)
177      INTEGER IPAR(20)
178      KPP_REAL RPAR(20)
179      INTEGER IERR
180!~~~>  The method parameters     
181      INTEGER Smax
182      PARAMETER (Smax = 6)
183      INTEGER  Method, ros_S
184      KPP_REAL ros_M(Smax), ros_E(Smax)
185      KPP_REAL ros_A(Smax*(Smax-1)/2), ros_C(Smax*(Smax-1)/2)
186      KPP_REAL ros_Alpha(Smax), ros_Gamma(Smax), ros_ELO
187      LOGICAL  ros_NewF(Smax)
188      CHARACTER*12 ros_Name
189!~~~>  Local variables     
190      KPP_REAL Roundoff,FacMin,FacMax,FacRej,FacSafe
191      KPP_REAL Hmin, Hmax, Hstart, Hexit
192      KPP_REAL Texit
193      INTEGER i, UplimTol, Max_no_steps
194      LOGICAL Autonomous, VectorTol
195!~~~>   Statistics on the work performed
196      INTEGER Nfun,Njac,Nstp,Nacc,Nrej,Ndec,Nsol,Nsng
197      COMMON /Statistics/ Nfun,Njac,Nstp,Nacc,Nrej,
198     &       Ndec,Nsol,Nsng
199!~~~>   Parameters
200      KPP_REAL ZERO, ONE, DeltaMin 
201      PARAMETER (ZERO = 0.0d0)
202      PARAMETER (ONE  = 1.0d0)
203      PARAMETER (DeltaMin = 1.0d-5)
204!~~~>   Functions
205      EXTERNAL ode_Fun, ode_Jac
206      KPP_REAL WLAMCH, ros_ErrorNorm
207      EXTERNAL WLAMCH, ros_ErrorNorm
208
209!~~~>  Initialize statistics
210      Nfun = IPAR(11)
211      Njac = IPAR(12)
212      Nstp = IPAR(13)
213      Nacc = IPAR(14)
214      Nrej = IPAR(15)
215      Ndec = IPAR(16)
216      Nsol = IPAR(17)
217      Nsng = IPAR(18)
218     
219!~~~>  Autonomous or time dependent ODE. Default is time dependent.
220      Autonomous = .NOT.(IPAR(1).EQ.0)
221
222!~~~>  For Scalar tolerances (IPAR(2).NE.0)  the code uses AbsTol(1) and RelTol(1)
223!      For Vector tolerances (IPAR(2).EQ.0) the code uses AbsTol(1:NVAR) and RelTol(1:NVAR)
224      IF (IPAR(2).EQ.0) THEN
225         VectorTol = .TRUE.
226         UplimTol  = KPP_NVAR
227      ELSE
228         VectorTol = .FALSE.
229         UplimTol  = 1
230      END IF
231     
232!~~~>   The maximum number of steps admitted
233      IF (IPAR(3).EQ.0) THEN
234         Max_no_steps = 100000
235      ELSEIF (Max_no_steps.GT.0) THEN
236         Max_no_steps=IPAR(3)
237      ELSE
238         WRITE(6,*)'User-selected max no. of steps: IPAR(3)=',IPAR(3)
239         CALL ros_ErrorMsg(-1,Tstart,ZERO,IERR)
240         RETURN         
241      END IF
242
243!~~~>  The particular Rosenbrock method chosen
244      IF (IPAR(4).EQ.0) THEN
245         Method = 3
246      ELSEIF ( (IPAR(4).GE.1).AND.(IPAR(4).LE.5) ) THEN
247         Method = IPAR(4)
248      ELSE 
249         WRITE (6,*) 'User-selected Rosenbrock method: IPAR(4)=', Method
250         CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR)
251         RETURN         
252      END IF
253     
254!~~~>  Unit roundoff (1+Roundoff>1) 
255      Roundoff = WLAMCH('E')
256
257!~~~>  Lower bound on the step size: (positive value)
258      IF (RPAR(1).EQ.ZERO) THEN
259         Hmin = ZERO
260      ELSEIF (RPAR(1).GT.ZERO) THEN
261         Hmin = RPAR(1)
262      ELSE       
263         WRITE (6,*) 'User-selected Hmin: RPAR(1)=', RPAR(1)
264         CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
265         RETURN         
266      END IF
267!~~~>  Upper bound on the step size: (positive value)
268      IF (RPAR(2).EQ.ZERO) THEN
269         Hmax = ABS(Tend-Tstart)
270      ELSEIF (RPAR(2).GT.ZERO) THEN
271         Hmax = MIN(ABS(RPAR(2)),ABS(Tend-Tstart))
272      ELSE       
273         WRITE (6,*) 'User-selected Hmax: RPAR(2)=', RPAR(2)
274         CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
275         RETURN         
276      END IF
277!~~~>  Starting step size: (positive value)
278      IF (RPAR(3).EQ.ZERO) THEN
279         Hstart = MAX(Hmin,DeltaMin)
280      ELSEIF (RPAR(3).GT.ZERO) THEN
281         Hstart = MIN(ABS(RPAR(3)),ABS(Tend-Tstart))
282      ELSE       
283         WRITE (6,*) 'User-selected Hstart: RPAR(3)=', RPAR(3)
284         CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
285         RETURN         
286      END IF
287!~~~>  Step size can be changed s.t.  FacMin < Hnew/Hexit < FacMax
288      IF (RPAR(4).EQ.ZERO) THEN
289         FacMin = 0.2d0
290      ELSEIF (RPAR(4).GT.ZERO) THEN
291         FacMin = RPAR(4)
292      ELSE       
293         WRITE (6,*) 'User-selected FacMin: RPAR(4)=', RPAR(4)
294         CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
295         RETURN         
296      END IF
297      IF (RPAR(5).EQ.ZERO) THEN
298         FacMax = 6.0d0
299      ELSEIF (RPAR(5).GT.ZERO) THEN
300         FacMax = RPAR(5)
301      ELSE       
302         WRITE (6,*) 'User-selected FacMax: RPAR(5)=', RPAR(5)
303         CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
304         RETURN         
305      END IF
306!~~~>   FacRej: Factor to decrease step after 2 succesive rejections
307      IF (RPAR(6).EQ.ZERO) THEN
308         FacRej = 0.1d0
309      ELSEIF (RPAR(6).GT.ZERO) THEN
310         FacRej = RPAR(6)
311      ELSE       
312         WRITE (6,*) 'User-selected FacRej: RPAR(6)=', RPAR(6)
313         CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
314         RETURN         
315      END IF
316!~~~>   FacSafe: Safety Factor in the computation of new step size
317      IF (RPAR(7).EQ.ZERO) THEN
318         FacSafe = 0.9d0
319      ELSEIF (RPAR(7).GT.ZERO) THEN
320         FacSafe = RPAR(7)
321      ELSE       
322         WRITE (6,*) 'User-selected FacSafe: RPAR(7)=', RPAR(7)
323         CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
324         RETURN         
325      END IF
326!~~~>  Check if tolerances are reasonable
327       DO i=1,UplimTol
328         IF ( (AbsTol(i).LE.ZERO) .OR. (RelTol(i).LE.10.d0*Roundoff)
329     &          .OR. (RelTol(i).GE.1.0d0) ) THEN
330            WRITE (6,*) ' AbsTol(',i,') = ',AbsTol(i)
331            WRITE (6,*) ' RelTol(',i,') = ',RelTol(i)
332            CALL ros_ErrorMsg(-5,Tstart,ZERO,IERR)
333            RETURN
334         END IF
335       END DO
336     
337 
338!~~~>   Initialize the particular Rosenbrock method
339
340      IF (Method .EQ. 1) THEN
341         CALL Ros2(ros_S, ros_A, ros_C, ros_M, ros_E, 
342     &             ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
343      ELSEIF (Method .EQ. 2) THEN
344         CALL Ros3(ros_S, ros_A, ros_C, ros_M, ros_E, 
345     &             ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
346      ELSEIF (Method .EQ. 3) THEN
347         CALL Ros4(ros_S, ros_A, ros_C, ros_M, ros_E, 
348     &             ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
349      ELSEIF (Method .EQ. 4) THEN
350         CALL Rodas3(ros_S, ros_A, ros_C, ros_M, ros_E, 
351     &             ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
352      ELSEIF (Method .EQ. 5) THEN
353         CALL Rodas4(ros_S, ros_A, ros_C, ros_M, ros_E, 
354     &             ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
355      ELSE
356         WRITE (6,*) 'Unknown Rosenbrock method: IPAR(4)=', Method
357         CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR) 
358         RETURN       
359      END IF
360
361!~~~>  CALL Rosenbrock method   
362      CALL RosenbrockIntegrator(Y,Tstart,Tend,Texit,
363     &      AbsTol,RelTol,
364     &      ode_Fun,ode_Jac ,
365!  Rosenbrock method coefficients     
366     &      ros_S, ros_M, ros_E, ros_A, ros_C, 
367     &      ros_Alpha, ros_Gamma, ros_ELO, ros_NewF,
368!  Integration parameters
369     &      Autonomous, VectorTol, Max_no_steps,
370     &      Roundoff, Hmin, Hmax, Hstart, Hexit,
371     &      FacMin, FacMax, FacRej, FacSafe, 
372!  Error indicator
373     &      IERR)
374
375
376!~~~>  Collect run statistics
377      IPAR(11) = Nfun
378      IPAR(12) = Njac
379      IPAR(13) = Nstp
380      IPAR(14) = Nacc
381      IPAR(15) = Nrej
382      IPAR(16) = Ndec
383      IPAR(17) = Nsol
384      IPAR(18) = Nsng
385!~~~> Last T and H
386      RPAR(11) = Texit
387      RPAR(12) = Hexit       
388     
389      RETURN
390      END  !  SUBROUTINE Rosenbrock
391
392     
393!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
394      SUBROUTINE RosenbrockIntegrator(Y,Tstart,Tend,T,
395     &      AbsTol,RelTol,
396     &      ode_Fun,ode_Jac ,
397!~~~> Rosenbrock method coefficients     
398     &      ros_S, ros_M, ros_E, ros_A, ros_C, 
399     &      ros_Alpha, ros_Gamma, ros_ELO, ros_NewF,
400!~~~> Integration parameters
401     &      Autonomous, VectorTol, Max_no_steps,
402     &      Roundoff, Hmin, Hmax, Hstart, Hexit,
403     &      FacMin, FacMax, FacRej, FacSafe, 
404!~~~> Error indicator
405     &      IERR)
406!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
407!   Template for the implementation of a generic Rosenbrock method
408!            defined by ros_S (no of stages) 
409!            and its coefficients ros_{A,C,M,E,Alpha,Gamma}
410!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
411
412      IMPLICIT NONE
413      INCLUDE 'KPP_ROOT_Parameters.h'
414      INCLUDE 'KPP_ROOT_Sparse.h'
415     
416!~~~> Input: the initial condition at Tstart; Output: the solution at T     
417      KPP_REAL Y(KPP_NVAR)
418!~~~> Input: integration interval   
419      KPP_REAL Tstart,Tend         
420!~~~> Output: time at which the solution is returned (T=Tend if success)   
421      KPP_REAL T         
422!~~~> Input: tolerances           
423      KPP_REAL  AbsTol(KPP_NVAR), RelTol(KPP_NVAR)
424!~~~> Input: ode function and its Jacobian           
425      EXTERNAL ode_Fun, ode_Jac
426!~~~> Input: The Rosenbrock method parameters     
427      INTEGER  ros_S
428      KPP_REAL ros_M(ros_S), ros_E(ros_S) 
429      KPP_REAL ros_A(ros_S*(ros_S-1)/2), ros_C(ros_S*(ros_S-1)/2)
430      KPP_REAL ros_Alpha(ros_S), ros_Gamma(ros_S), ros_ELO
431      LOGICAL ros_NewF(ros_S)
432!~~~> Input: integration parameters     
433      LOGICAL Autonomous, VectorTol
434      KPP_REAL Hstart, Hmin, Hmax
435      INTEGER Max_no_steps
436      KPP_REAL Roundoff, FacMin, FacMax, FacRej, FacSafe 
437!~~~> Output: last accepted step     
438      KPP_REAL Hexit 
439!~~~> Output: Error indicator
440      INTEGER IERR
441! ~~~~ Local variables           
442      KPP_REAL Ynew(KPP_NVAR), Fcn0(KPP_NVAR), Fcn(KPP_NVAR),
443     &       K(KPP_NVAR*ros_S), dFdT(KPP_NVAR),
444     &       Jac0(KPP_LU_NONZERO), Ghimj(KPP_LU_NONZERO)
445      KPP_REAL H, Hnew, HC, HG, Fac, Tau 
446      KPP_REAL Err, Yerr(KPP_NVAR)
447      INTEGER Pivot(KPP_NVAR), Direction, ioffset, j, istage
448      LOGICAL RejectLastH, RejectMoreH, Singular
449!~~~>  Local parameters
450      KPP_REAL ZERO, ONE, DeltaMin 
451      PARAMETER (ZERO = 0.0d0)
452      PARAMETER (ONE  = 1.0d0)
453      PARAMETER (DeltaMin = 1.0d-5)
454!~~~>  Locally called functions
455      KPP_REAL WLAMCH, ros_ErrorNorm
456      EXTERNAL WLAMCH, ros_ErrorNorm
457!~~~>  Statistics on the work performed
458      INTEGER Nfun,Njac,Nstp,Nacc,Nrej,Ndec,Nsol,Nsng
459      COMMON /Statistics/ Nfun,Njac,Nstp,Nacc,Nrej,
460     &       Ndec,Nsol,Nsng
461!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
462
463     
464!~~~>  INITIAL PREPARATIONS
465      T = Tstart
466      Hexit = 0.0d0
467      H = MIN(Hstart,Hmax) 
468      IF (ABS(H).LE.10.d0*Roundoff) THEN
469           H = DeltaMin
470      END IF       
471     
472      IF (Tend .GE. Tstart) THEN
473        Direction = +1
474      ELSE
475        Direction = -1
476      END IF           
477
478      RejectLastH=.FALSE.
479      RejectMoreH=.FALSE.
480   
481!~~~> Time loop begins below
482
483      DO WHILE ( (Direction.GT.0).AND.((T-Tend)+Roundoff.LE.ZERO)
484     &     .OR. (Direction.LT.0).AND.((Tend-T)+Roundoff.LE.ZERO) ) 
485         
486      IF ( Nstp.GT.Max_no_steps ) THEN  ! Too many steps
487        CALL ros_ErrorMsg(-6,T,H,IERR)
488        RETURN
489      END IF
490      IF ( ((T+0.1d0*H).EQ.T).OR.(H.LE.Roundoff) ) THEN  ! Step size too small
491        CALL ros_ErrorMsg(-7,T,H,IERR)
492        RETURN
493      END IF
494     
495!~~~>  Limit H if necessary to avoid going beyond Tend   
496      Hexit = H
497      H = MIN(H,ABS(Tend-T))
498
499!~~~>   Compute the function at current time
500      CALL ode_Fun(T,Y,Fcn0)
501
502!~~~>  Compute the function derivative with respect to T
503      IF (.NOT.Autonomous) THEN
504         CALL ros_FunTimeDerivative ( T, Roundoff, Y, 
505     &                       Fcn0, ode_Fun, dFdT )
506      END IF
507 
508!~~~>   Compute the Jacobian at current time
509      CALL ode_Jac(T,Y,Jac0)
510 
511!~~~>  Repeat step calculation until current step accepted
512      DO WHILE (.TRUE.) ! WHILE STEP NOT ACCEPTED
513
514     
515      CALL ros_PrepareMatrix(H,Direction,ros_Gamma(1),
516     &              Jac0,Ghimj,Pivot,Singular)
517      IF (Singular) THEN ! More than 5 consecutive failed decompositions
518         CALL ros_ErrorMsg(-8,T,H,IERR)
519         RETURN
520      END IF
521
522!~~~>   Compute the stages
523      DO istage = 1, ros_S
524         
525         ! Current istage offset. Current istage vector is K(ioffset+1:ioffset+KPP_NVAR)
526         ioffset = KPP_NVAR*(istage-1)
527         
528         ! For the 1st istage the function has been computed previously
529         IF ( istage.EQ.1 ) THEN
530           CALL WCOPY(KPP_NVAR,Fcn0,1,Fcn,1)
531         ! istage>1 and a new function evaluation is needed at the current istage
532         ELSEIF ( ros_NewF(istage) ) THEN
533           CALL WCOPY(KPP_NVAR,Y,1,Ynew,1)
534           DO j = 1, istage-1
535             CALL WAXPY(KPP_NVAR,ros_A((istage-1)*(istage-2)/2+j),
536     &                   K(KPP_NVAR*(j-1)+1),1,Ynew,1) 
537           END DO
538           Tau = T + ros_Alpha(istage)*Direction*H
539           CALL ode_Fun(Tau,Ynew,Fcn)
540         END IF ! if istage.EQ.1 elseif ros_NewF(istage)
541         CALL WCOPY(KPP_NVAR,Fcn,1,K(ioffset+1),1)
542         DO j = 1, istage-1
543           HC = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H)
544           CALL WAXPY(KPP_NVAR,HC,K(KPP_NVAR*(j-1)+1),1,K(ioffset+1),1)
545         END DO
546         IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN
547           HG = Direction*H*ros_Gamma(istage)
548           CALL WAXPY(KPP_NVAR,HG,dFdT,1,K(ioffset+1),1)
549         END IF
550         CALL SolveTemplate(Ghimj, Pivot, K(ioffset+1))
551         
552      END DO  ! istage     
553           
554
555!~~~>  Compute the new solution
556      CALL WCOPY(KPP_NVAR,Y,1,Ynew,1)
557      DO j=1,ros_S
558         CALL WAXPY(KPP_NVAR,ros_M(j),K(KPP_NVAR*(j-1)+1),1,Ynew,1)
559      END DO
560
561!~~~>  Compute the error estimation
562      CALL WSCAL(KPP_NVAR,ZERO,Yerr,1)
563      DO j=1,ros_S       
564        CALL WAXPY(KPP_NVAR,ros_E(j),K(KPP_NVAR*(j-1)+1),1,Yerr,1)
565      END DO
566      Err = ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol )
567
568!~~~> New step size is bounded by FacMin <= Hnew/H <= FacMax
569      Fac  = MIN(FacMax,MAX(FacMin,FacSafe/Err**(ONE/ros_ELO)))
570      Hnew = H*Fac 
571
572!~~~>  Check the error magnitude and adjust step size
573      Nstp = Nstp+1
574      IF ( (Err.LE.ONE).OR.(H.LE.Hmin) ) THEN  !~~~> Accept step
575         Nacc = Nacc+1
576         CALL WCOPY(KPP_NVAR,Ynew,1,Y,1)
577         T = T + Direction*H
578         Hnew = MAX(Hmin,MIN(Hnew,Hmax))
579         IF (RejectLastH) THEN  ! No step size increase after a rejected step
580            Hnew = MIN(Hnew,H) 
581         END IF   
582         RejectLastH = .FALSE. 
583         RejectMoreH = .FALSE.
584         H = Hnew
585         GOTO 101  ! EXIT THE LOOP: WHILE STEP NOT ACCEPTED
586      ELSE                 !~~~> Reject step
587         IF (RejectMoreH) THEN
588            Hnew=H*FacRej
589         END IF   
590         RejectMoreH = RejectLastH
591         RejectLastH = .TRUE.
592         H = Hnew
593         IF (Nacc.GE.1) THEN
594            Nrej = Nrej+1
595         END IF   
596      END IF ! Err <= 1
597
598      END DO ! LOOP: WHILE STEP NOT ACCEPTED
599
600101   CONTINUE
601
602      END DO ! Time loop   
603     
604!~~~> Succesful exit
605      IERR = 1  !~~~> The integration was successful
606
607      RETURN
608      END  !  SUBROUTINE RosenbrockIntegrator
609 
610     
611!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
612      KPP_REAL FUNCTION ros_ErrorNorm ( Y, Ynew, Yerr, 
613     &                      AbsTol, RelTol, VectorTol )
614!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
615!~~~> Computes the "scaled norm" of the error vector Yerr
616!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
617      IMPLICIT NONE     
618      INCLUDE 'KPP_ROOT_Parameters.h'
619
620! Input arguments   
621      KPP_REAL Y(KPP_NVAR), Ynew(KPP_NVAR), Yerr(KPP_NVAR)
622      KPP_REAL AbsTol(KPP_NVAR), RelTol(KPP_NVAR)   
623      LOGICAL  VectorTol
624! Local variables     
625      KPP_REAL Err, Scale, Ymax, ZERO   
626      INTEGER i
627      PARAMETER (ZERO = 0.0d0)
628     
629      Err = ZERO
630      DO i=1,KPP_NVAR
631        Ymax = MAX(ABS(Y(i)),ABS(Ynew(i)))
632        IF (VectorTol) THEN
633          Scale = AbsTol(i)+RelTol(i)*Ymax
634        ELSE
635          Scale = AbsTol(1)+RelTol(1)*Ymax
636        END IF
637        Err = Err+(Yerr(i)/Scale)**2
638      END DO
639      Err  = SQRT(Err/KPP_NVAR)
640
641      ros_ErrorNorm = Err
642     
643      RETURN
644      END ! FUNCTION ros_ErrorNorm
645
646!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
647      SUBROUTINE ros_FunTimeDerivative ( T, Roundoff, Y, 
648     &                       Fcn0, ode_Fun, dFdT )
649!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
650!~~~> The time partial derivative of the function by finite differences
651!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
652      IMPLICIT NONE     
653      INCLUDE 'KPP_ROOT_Parameters.h'
654
655!~~~> Input arguments   
656      KPP_REAL T, Roundoff, Y(KPP_NVAR), Fcn0(KPP_NVAR) 
657      EXTERNAL ode_Fun 
658!~~~> Output arguments     
659      KPP_REAL dFdT(KPP_NVAR)   
660!~~~> Global variables     
661      INTEGER Nfun,Njac,Nstp,Nacc,Nrej,Ndec,Nsol,Nsng
662      COMMON /Statistics/ Nfun,Njac,Nstp,Nacc,Nrej,
663     &       Ndec,Nsol,Nsng     
664!~~~> Local variables     
665      KPP_REAL Delta, DeltaMin, ONE     
666      PARAMETER ( DeltaMin = 1.0d-6 )   
667      PARAMETER ( ONE = 1.0d0 )
668     
669      Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T))
670      CALL ode_Fun(T+Delta,Y,dFdT)
671      CALL WAXPY(KPP_NVAR,(-ONE),Fcn0,1,dFdT,1)
672      CALL WSCAL(KPP_NVAR,(ONE/Delta),dFdT,1)
673
674      RETURN
675      END ! SUBROUTINE ros_FunTimeDerivative
676
677
678!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
679      SUBROUTINE ros_PrepareMatrix ( H, Direction, gam, 
680     &                    Jac0, Ghimj, Pivot, Singular )
681! --- --- --- --- --- --- --- --- --- --- --- --- ---
682!  Prepares the LHS matrix for stage calculations
683!  1.  Construct Ghimj = 1/(H*ham) - Jac0
684!         "(Gamma H) Inverse Minus Jacobian"
685!  2.  Repeat LU decomposition of Ghimj until successful.
686!          -half the step size if LU decomposition fails and retry
687!          -exit after 5 consecutive fails
688! --- --- --- --- --- --- --- --- --- --- --- --- ---
689      IMPLICIT NONE     
690      INCLUDE 'KPP_ROOT_Parameters.h'
691      INCLUDE 'KPP_ROOT_Sparse.h'
692     
693!~~~> Input arguments     
694      KPP_REAL gam, Jac0(KPP_LU_NONZERO)
695      INTEGER  Direction
696!~~~> Output arguments     
697      KPP_REAL Ghimj(KPP_LU_NONZERO)
698      LOGICAL  Singular
699      INTEGER  Pivot(KPP_NVAR)
700!~~~> Inout arguments     
701      KPP_REAL H      ! step size is decreased when LU fails
702!~~~> Global variables     
703      INTEGER Nfun,Njac,Nstp,Nacc,Nrej,Ndec,Nsol,Nsng
704      COMMON /Statistics/ Nfun,Njac,Nstp,Nacc,Nrej,
705     &       Ndec,Nsol,Nsng     
706!~~~> Local variables     
707      INTEGER  i, ising, Nconsecutive
708      KPP_REAL  ghinv, ONE, HALF
709      PARAMETER ( ONE  = 1.0d0 )
710      PARAMETER ( HALF = 0.5d0 )
711     
712      Nconsecutive = 0
713      Singular = .TRUE.
714     
715      DO WHILE (Singular)
716     
717!~~~>    Construct Ghimj = 1/(H*ham) - Jac0
718        CALL WCOPY(KPP_LU_NONZERO,Jac0,1,Ghimj,1)
719        CALL WSCAL(KPP_LU_NONZERO,(-ONE),Ghimj,1)
720        ghinv = ONE/(Direction*H*gam)
721        DO i=1,KPP_NVAR
722          Ghimj(LU_DIAG(i)) = Ghimj(LU_DIAG(i))+ghinv
723        END DO
724!~~~>    Compute LU decomposition
725        CALL DecompTemplate( Ghimj, Pivot, ising )
726        IF (ising .EQ. 0) THEN
727!~~~>    If successful done
728          Singular = .FALSE. 
729        ELSE ! ising .ne. 0
730!~~~>    If unsuccessful half the step size; if 5 consecutive fails then return
731          Nsng = Nsng+1
732          Nconsecutive = Nconsecutive+1
733          Singular = .TRUE. 
734          PRINT*,'Warning: LU Decomposition returned ising = ',ising
735          IF (Nconsecutive.LE.5) THEN ! Less than 5 consecutive failed decompositions
736            H = H*HALF
737          ELSE  ! More than 5 consecutive failed decompositions
738            RETURN
739          END IF  ! Nconsecutive
740        END IF  ! ising
741         
742      END DO ! WHILE Singular
743
744      RETURN
745      END ! SUBROUTINE ros_PrepareMatrix
746
747
748!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
749      SUBROUTINE ros_ErrorMsg(Code,T,H,IERR)
750      KPP_REAL T, H
751      INTEGER IERR, Code
752     
753      IERR = Code
754      WRITE(6,*) 
755     &   'Forced exit from Rosenbrock due to the following error:' 
756     
757      IF    (Code .EQ. -1) THEN   
758         WRITE(6,*) '--> Improper value for maximal no of steps'
759      ELSEIF (Code .EQ. -2) THEN   
760         WRITE(6,*) '--> Selected Rosenbrock method not implemented'
761      ELSEIF (Code .EQ. -3) THEN   
762         WRITE(6,*) '--> Hmin/Hmax/Hstart must be positive'
763      ELSEIF (Code .EQ. -4) THEN   
764         WRITE(6,*) '--> FacMin/FacMax/FacRej must be positive'
765      ELSEIF (Code .EQ. -5) THEN
766         WRITE(6,*) '--> Improper tolerance values'
767      ELSEIF (Code .EQ. -6) THEN
768         WRITE(6,*) '--> No of steps exceeds maximum bound'
769      ELSEIF (Code .EQ. -7) THEN
770         WRITE(6,*) '--> Step size too small: T + 10*H = T',
771     &                ' or H < Roundoff'
772      ELSEIF (Code .EQ. -8) THEN   
773         WRITE(6,*) '--> Matrix is repeatedly singular'
774      ELSE
775         WRITE(6,102) 'Unknown Error code: ',Code
776      END IF
777     
778  102 FORMAT('       ',A,I4)   
779      WRITE(6,103) T, H
780     
781  103 FORMAT('        T=',E15.7,' and H=',E15.7)   
782     
783      RETURN
784      END
785     
786     
787
788!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
789      SUBROUTINE Ros2 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,
790     &                ros_Gamma,ros_NewF,ros_ELO,ros_Name)
791!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
792! --- AN L-STABLE METHOD, 2 stages, order 2
793!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
794      IMPLICIT NONE
795      INTEGER S
796      PARAMETER (S=2)
797      INTEGER  ros_S
798      KPP_REAL ros_M(S), ros_E(S), ros_A(S*(S-1)/2), ros_C(S*(S-1)/2)
799      KPP_REAL ros_Alpha(S), ros_Gamma(S), ros_ELO
800      LOGICAL  ros_NewF(S)
801      CHARACTER*12 ros_Name
802      DOUBLE PRECISION g
803     
804       g = 1.0d0 + 1.0d0/SQRT(2.0d0)
805     
806!~~~> Name of the method
807       ros_Name = 'ROS-2'     
808!~~~> Number of stages
809       ros_S = 2
810     
811!~~~> The coefficient matrices A and C are strictly lower triangular.
812!   The lower triangular (subdiagonal) elements are stored in row-wise order:
813!   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
814!   The general mapping formula is:
815!             A(i,j) = ros_A( (i-1)*(i-2)/2 + j )       
816!             C(i,j) = ros_C( (i-1)*(i-2)/2 + j ) 
817     
818       ros_A(1) = (1.d0)/g
819       ros_C(1) = (-2.d0)/g
820!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
821!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
822       ros_NewF(1) = .TRUE.
823       ros_NewF(2) = .TRUE.
824!~~~> M_i = Coefficients for new step solution
825       ros_M(1)= (3.d0)/(2.d0*g)
826       ros_M(2)= (1.d0)/(2.d0*g)
827! E_i = Coefficients for error estimator       
828       ros_E(1) = 1.d0/(2.d0*g)
829       ros_E(2) = 1.d0/(2.d0*g)
830!~~~> ros_ELO = estimator of local order - the minimum between the
831!       main and the embedded scheme orders plus one
832       ros_ELO = 2.0d0       
833!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
834       ros_Alpha(1) = 0.0d0
835       ros_Alpha(2) = 1.0d0 
836!~~~> Gamma_i = \sum_j  gamma_{i,j}       
837       ros_Gamma(1) = g
838       ros_Gamma(2) =-g
839     
840      RETURN
841      END ! SUBROUTINE Ros2
842
843
844!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
845      SUBROUTINE Ros3 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,
846     &               ros_Gamma,ros_NewF,ros_ELO,ros_Name)
847!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
848! --- AN L-STABLE METHOD, 3 stages, order 3, 2 function evaluations
849!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
850      IMPLICIT NONE
851      INTEGER S
852      PARAMETER (S=3)
853      INTEGER  ros_S
854      KPP_REAL ros_M(S), ros_E(S), ros_A(S*(S-1)/2), ros_C(S*(S-1)/2)
855      KPP_REAL ros_Alpha(S), ros_Gamma(S), ros_ELO
856      LOGICAL  ros_NewF(S)
857      CHARACTER*12 ros_Name
858     
859!~~~> Name of the method
860      ros_Name = 'ROS-3'     
861!~~~> Number of stages
862      ros_S = 3
863     
864!~~~> The coefficient matrices A and C are strictly lower triangular.
865!   The lower triangular (subdiagonal) elements are stored in row-wise order:
866!   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
867!   The general mapping formula is:
868!             A(i,j) = ros_A( (i-1)*(i-2)/2 + j )       
869!             C(i,j) = ros_C( (i-1)*(i-2)/2 + j ) 
870     
871      ros_A(1)= 1.d0
872      ros_A(2)= 1.d0
873      ros_A(3)= 0.d0
874
875      ros_C(1) = -0.10156171083877702091975600115545d+01
876      ros_C(2) =  0.40759956452537699824805835358067d+01
877      ros_C(3) =  0.92076794298330791242156818474003d+01
878!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
879!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
880      ros_NewF(1) = .TRUE.
881      ros_NewF(2) = .TRUE.
882      ros_NewF(3) = .FALSE.
883!~~~> M_i = Coefficients for new step solution
884      ros_M(1) =  0.1d+01
885      ros_M(2) =  0.61697947043828245592553615689730d+01
886      ros_M(3) = -0.42772256543218573326238373806514d+00
887! E_i = Coefficients for error estimator       
888      ros_E(1) =  0.5d+00
889      ros_E(2) = -0.29079558716805469821718236208017d+01
890      ros_E(3) =  0.22354069897811569627360909276199d+00
891!~~~> ros_ELO = estimator of local order - the minimum between the
892!       main and the embedded scheme orders plus 1
893      ros_ELO = 3.0d0       
894!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
895      ros_Alpha(1)= 0.0d+00
896      ros_Alpha(2)= 0.43586652150845899941601945119356d+00
897      ros_Alpha(3)= 0.43586652150845899941601945119356d+00
898!~~~> Gamma_i = \sum_j  gamma_{i,j}       
899      ros_Gamma(1)= 0.43586652150845899941601945119356d+00
900      ros_Gamma(2)= 0.24291996454816804366592249683314d+00
901      ros_Gamma(3)= 0.21851380027664058511513169485832d+01
902      RETURN
903      END ! SUBROUTINE Ros3
904
905!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
906
907
908!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
909      SUBROUTINE Ros4 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,
910     &               ros_Gamma,ros_NewF,ros_ELO,ros_Name)
911!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
912!     L-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 4 STAGES
913!     L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3
914!
915!         E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
916!         EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
917!         SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
918!         SPRINGER-VERLAG (1990)               
919!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
920
921      IMPLICIT NONE
922      INTEGER S
923      PARAMETER (S=4)
924      INTEGER  ros_S
925      KPP_REAL ros_M(S), ros_E(S), ros_A(S*(S-1)/2), ros_C(S*(S-1)/2)
926      KPP_REAL ros_Alpha(S), ros_Gamma(S), ros_ELO
927      LOGICAL  ros_NewF(S)
928      CHARACTER*12 ros_Name
929     
930!~~~> Name of the method
931      ros_Name = 'ROS-4'     
932!~~~> Number of stages
933      ros_S = 4
934     
935!~~~> The coefficient matrices A and C are strictly lower triangular.
936!   The lower triangular (subdiagonal) elements are stored in row-wise order:
937!   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
938!   The general mapping formula is:
939!             A(i,j) = ros_A( (i-1)*(i-2)/2 + j )       
940!             C(i,j) = ros_C( (i-1)*(i-2)/2 + j ) 
941     
942      ros_A(1) = 0.2000000000000000d+01
943      ros_A(2) = 0.1867943637803922d+01
944      ros_A(3) = 0.2344449711399156d+00
945      ros_A(4) = ros_A(2)
946      ros_A(5) = ros_A(3)
947      ros_A(6) = 0.0D0
948
949      ros_C(1) =-0.7137615036412310d+01
950      ros_C(2) = 0.2580708087951457d+01
951      ros_C(3) = 0.6515950076447975d+00
952      ros_C(4) =-0.2137148994382534d+01
953      ros_C(5) =-0.3214669691237626d+00
954      ros_C(6) =-0.6949742501781779d+00
955!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
956!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
957      ros_NewF(1)  = .TRUE.
958      ros_NewF(2)  = .TRUE.
959      ros_NewF(3)  = .TRUE.
960      ros_NewF(4)  = .FALSE.
961!~~~> M_i = Coefficients for new step solution
962      ros_M(1) = 0.2255570073418735d+01
963      ros_M(2) = 0.2870493262186792d+00
964      ros_M(3) = 0.4353179431840180d+00
965      ros_M(4) = 0.1093502252409163d+01
966!~~~> E_i  = Coefficients for error estimator       
967      ros_E(1) =-0.2815431932141155d+00
968      ros_E(2) =-0.7276199124938920d-01
969      ros_E(3) =-0.1082196201495311d+00
970      ros_E(4) =-0.1093502252409163d+01
971!~~~> ros_ELO  = estimator of local order - the minimum between the
972!       main and the embedded scheme orders plus 1
973      ros_ELO  = 4.0d0       
974!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
975      ros_Alpha(1) = 0.D0
976      ros_Alpha(2) = 0.1145640000000000d+01
977      ros_Alpha(3) = 0.6552168638155900d+00
978      ros_Alpha(4) = ros_Alpha(3)
979!~~~> Gamma_i = \sum_j  gamma_{i,j}       
980      ros_Gamma(1) = 0.5728200000000000d+00
981      ros_Gamma(2) =-0.1769193891319233d+01
982      ros_Gamma(3) = 0.7592633437920482d+00
983      ros_Gamma(4) =-0.1049021087100450d+00
984      RETURN
985      END ! SUBROUTINE Ros4
986     
987!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
988      SUBROUTINE Rodas3 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,
989     &                ros_Gamma,ros_NewF,ros_ELO,ros_Name)
990!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
991! --- A STIFFLY-STABLE METHOD, 4 stages, order 3
992!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
993      IMPLICIT NONE
994      INTEGER S
995      PARAMETER (S=4)
996      INTEGER  ros_S
997      KPP_REAL ros_M(S), ros_E(S), ros_A(S*(S-1)/2), ros_C(S*(S-1)/2)
998      KPP_REAL ros_Alpha(S), ros_Gamma(S), ros_ELO
999      LOGICAL  ros_NewF(S)
1000      CHARACTER*12 ros_Name
1001     
1002!~~~> Name of the method
1003      ros_Name = 'RODAS-3'     
1004!~~~> Number of stages
1005      ros_S = 4
1006     
1007!~~~> The coefficient matrices A and C are strictly lower triangular.
1008!   The lower triangular (subdiagonal) elements are stored in row-wise order:
1009!   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
1010!   The general mapping formula is:
1011!             A(i,j) = ros_A( (i-1)*(i-2)/2 + j )       
1012!             C(i,j) = ros_C( (i-1)*(i-2)/2 + j ) 
1013 
1014      ros_A(1) = 0.0d+00
1015      ros_A(2) = 2.0d+00
1016      ros_A(3) = 0.0d+00
1017      ros_A(4) = 2.0d+00
1018      ros_A(5) = 0.0d+00
1019      ros_A(6) = 1.0d+00
1020
1021      ros_C(1) = 4.0d+00
1022      ros_C(2) = 1.0d+00
1023      ros_C(3) =-1.0d+00
1024      ros_C(4) = 1.0d+00
1025      ros_C(5) =-1.0d+00 
1026      ros_C(6) =-(8.0d+00/3.0d+00) 
1027               
1028!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
1029!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
1030      ros_NewF(1)  = .TRUE.
1031      ros_NewF(2)  = .FALSE.
1032      ros_NewF(3)  = .TRUE.
1033      ros_NewF(4)  = .TRUE.
1034!~~~> M_i = Coefficients for new step solution
1035      ros_M(1) = 2.0d+00
1036      ros_M(2) = 0.0d+00
1037      ros_M(3) = 1.0d+00
1038      ros_M(4) = 1.0d+00
1039!~~~> E_i  = Coefficients for error estimator       
1040      ros_E(1) = 0.0d+00
1041      ros_E(2) = 0.0d+00
1042      ros_E(3) = 0.0d+00
1043      ros_E(4) = 1.0d+00
1044!~~~> ros_ELO  = estimator of local order - the minimum between the
1045!       main and the embedded scheme orders plus 1
1046      ros_ELO  = 3.0d+00       
1047!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
1048      ros_Alpha(1) = 0.0d+00
1049      ros_Alpha(2) = 0.0d+00
1050      ros_Alpha(3) = 1.0d+00
1051      ros_Alpha(4) = 1.0d+00
1052!~~~> Gamma_i = \sum_j  gamma_{i,j}       
1053      ros_Gamma(1) = 0.5d+00
1054      ros_Gamma(2) = 1.5d+00
1055      ros_Gamma(3) = 0.0d+00
1056      ros_Gamma(4) = 0.0d+00
1057      RETURN
1058      END ! SUBROUTINE Rodas3
1059   
1060!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
1061      SUBROUTINE Rodas4 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,
1062     &                 ros_Gamma,ros_NewF,ros_ELO,ros_Name)
1063!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
1064!     STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 6 STAGES
1065!
1066!         E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
1067!         EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
1068!         SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
1069!         SPRINGER-VERLAG (1996)               
1070!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
1071
1072      IMPLICIT NONE
1073      INTEGER S
1074      PARAMETER (S=6)
1075      INTEGER  ros_S
1076      KPP_REAL ros_M(S), ros_E(S), ros_A(S*(S-1)/2), ros_C(S*(S-1)/2)
1077      KPP_REAL ros_Alpha(S), ros_Gamma(S), ros_ELO
1078      LOGICAL  ros_NewF(S)
1079      CHARACTER*12 ros_Name
1080
1081!~~~> Name of the method
1082       ros_Name = 'RODAS-4'     
1083!~~~> Number of stages
1084       ros_S = 6
1085
1086!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
1087       ros_Alpha(1) = 0.000d0
1088       ros_Alpha(2) = 0.386d0
1089       ros_Alpha(3) = 0.210d0 
1090       ros_Alpha(4) = 0.630d0
1091       ros_Alpha(5) = 1.000d0
1092       ros_Alpha(6) = 1.000d0
1093       
1094!~~~> Gamma_i = \sum_j  gamma_{i,j}       
1095       ros_Gamma(1) = 0.2500000000000000d+00
1096       ros_Gamma(2) =-0.1043000000000000d+00
1097       ros_Gamma(3) = 0.1035000000000000d+00
1098       ros_Gamma(4) =-0.3620000000000023d-01
1099       ros_Gamma(5) = 0.0d0
1100       ros_Gamma(6) = 0.0d0
1101
1102!~~~> The coefficient matrices A and C are strictly lower triangular.
1103!   The lower triangular (subdiagonal) elements are stored in row-wise order:
1104!   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
1105!   The general mapping formula is:  A(i,j) = ros_A( (i-1)*(i-2)/2 + j )       
1106!                                    C(i,j) = ros_C( (i-1)*(i-2)/2 + j ) 
1107     
1108       ros_A(1) = 0.1544000000000000d+01
1109       ros_A(2) = 0.9466785280815826d+00
1110       ros_A(3) = 0.2557011698983284d+00
1111       ros_A(4) = 0.3314825187068521d+01
1112       ros_A(5) = 0.2896124015972201d+01
1113       ros_A(6) = 0.9986419139977817d+00
1114       ros_A(7) = 0.1221224509226641d+01
1115       ros_A(8) = 0.6019134481288629d+01
1116       ros_A(9) = 0.1253708332932087d+02
1117       ros_A(10) =-0.6878860361058950d+00
1118       ros_A(11) = ros_A(7)
1119       ros_A(12) = ros_A(8)
1120       ros_A(13) = ros_A(9)
1121       ros_A(14) = ros_A(10)
1122       ros_A(15) = 1.0d+00
1123
1124       ros_C(1) =-0.5668800000000000d+01
1125       ros_C(2) =-0.2430093356833875d+01
1126       ros_C(3) =-0.2063599157091915d+00
1127       ros_C(4) =-0.1073529058151375d+00
1128       ros_C(5) =-0.9594562251023355d+01
1129       ros_C(6) =-0.2047028614809616d+02
1130       ros_C(7) = 0.7496443313967647d+01
1131       ros_C(8) =-0.1024680431464352d+02
1132       ros_C(9) =-0.3399990352819905d+02
1133       ros_C(10) = 0.1170890893206160d+02
1134       ros_C(11) = 0.8083246795921522d+01
1135       ros_C(12) =-0.7981132988064893d+01
1136       ros_C(13) =-0.3152159432874371d+02
1137       ros_C(14) = 0.1631930543123136d+02
1138       ros_C(15) =-0.6058818238834054d+01
1139
1140!~~~> M_i = Coefficients for new step solution
1141       ros_M(1) = ros_A(7)
1142       ros_M(2) = ros_A(8)
1143       ros_M(3) = ros_A(9)
1144       ros_M(4) = ros_A(10)
1145       ros_M(5) = 1.0d+00
1146       ros_M(6) = 1.0d+00
1147
1148!~~~> E_i  = Coefficients for error estimator       
1149       ros_E(1) = 0.0d+00
1150       ros_E(2) = 0.0d+00
1151       ros_E(3) = 0.0d+00
1152       ros_E(4) = 0.0d+00
1153       ros_E(5) = 0.0d+00
1154       ros_E(6) = 1.0d+00
1155
1156!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
1157!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
1158       ros_NewF(1) = .TRUE.
1159       ros_NewF(2) = .TRUE.
1160       ros_NewF(3) = .TRUE.
1161       ros_NewF(4) = .TRUE.
1162       ros_NewF(5) = .TRUE.
1163       ros_NewF(6) = .TRUE.
1164     
1165!~~~> ros_ELO  = estimator of local order - the minimum between the
1166!       main and the embedded scheme orders plus 1
1167       ros_ELO = 4.0d0
1168     
1169      RETURN
1170      END ! SUBROUTINE Rodas4
1171
1172     
1173
1174!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
1175      SUBROUTINE DecompTemplate( A, Pivot, ising )
1176!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
1177!  Template for the LU decomposition   
1178!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
1179      INCLUDE 'KPP_ROOT_Parameters.h'
1180      INCLUDE 'KPP_ROOT_Global.h'
1181!~~~> Inout variables     
1182      KPP_REAL A(KPP_LU_NONZERO)
1183!~~~> Output variables     
1184      INTEGER Pivot(KPP_NVAR), ising
1185!~~~> Collect statistics     
1186      INTEGER Nfun,Njac,Nstp,Nacc,Nrej,Ndec,Nsol,Nsng
1187      COMMON /Statistics/ Nfun,Njac,Nstp,Nacc,Nrej,
1188     &       Ndec,Nsol,Nsng
1189     
1190      CALL KppDecomp ( A, ising )
1191!~~~> Note: for a full matrix use Lapack:
1192!        CALL  DGETRF( KPP_NVAR, KPP_NVAR, A, KPP_NVAR, Pivot, ising )
1193   
1194      Ndec = Ndec + 1
1195
1196      END ! SUBROUTINE DecompTemplate
1197 
1198!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
1199      SUBROUTINE SolveTemplate( A, Pivot, b )
1200!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
1201!  Template for the forward/backward substitution (using pre-computed LU decomposition)   
1202!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
1203      INCLUDE 'KPP_ROOT_Parameters.h'
1204      INCLUDE 'KPP_ROOT_Global.h'
1205!~~~> Input variables     
1206      KPP_REAL A(KPP_LU_NONZERO)
1207      INTEGER Pivot(KPP_NVAR)
1208!~~~> InOut variables     
1209      KPP_REAL b(KPP_NVAR)
1210!~~~> Collect statistics     
1211      INTEGER Nfun,Njac,Nstp,Nacc,Nrej,Ndec,Nsol,Nsng
1212      COMMON /Statistics/ Nfun,Njac,Nstp,Nacc,Nrej,
1213     &       Ndec,Nsol,Nsng
1214     
1215      CALL KppSolve( A, b )
1216!~~~> Note: for a full matrix use Lapack:
1217!        NRHS = 1
1218!        CALL  DGETRS( 'N', KPP_NVAR , NRHS, A, KPP_NVAR, Pivot, b, KPP_NVAR, INFO )
1219     
1220      Nsol = Nsol+1
1221
1222      END ! SUBROUTINE SolveTemplate
1223
1224!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
1225      SUBROUTINE FunTemplate( T, Y, Ydot )
1226!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
1227!  Template for the ODE function call.
1228!  Updates the rate coefficients (and possibly the fixed species) at each call   
1229!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
1230      INCLUDE 'KPP_ROOT_Parameters.h'
1231      INCLUDE 'KPP_ROOT_Global.h'
1232!~~~> Input variables     
1233      KPP_REAL T, Y(KPP_NVAR)
1234!~~~> Output variables     
1235      KPP_REAL Ydot(KPP_NVAR)
1236!~~~> Local variables
1237      KPP_REAL Told     
1238!~~~> Collect statistics     
1239      INTEGER Nfun,Njac,Nstp,Nacc,Nrej,Ndec,Nsol,Nsng
1240      COMMON /Statistics/ Nfun,Njac,Nstp,Nacc,Nrej,
1241     &       Ndec,Nsol,Nsng
1242
1243      Told = TIME
1244      TIME = T
1245      CALL Update_SUN()
1246      CALL Update_RCONST()
1247      CALL Fun( Y, FIX, RCONST, Ydot )
1248      TIME = Told
1249     
1250      Nfun = Nfun+1
1251     
1252      RETURN
1253      END ! SUBROUTINE FunTemplate
1254
1255 
1256!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
1257      SUBROUTINE JacTemplate( T, Y, Jcb )
1258!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
1259!  Template for the ODE Jacobian call.
1260!  Updates the rate coefficients (and possibly the fixed species) at each call   
1261!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
1262      INCLUDE 'KPP_ROOT_Parameters.h'
1263      INCLUDE 'KPP_ROOT_Global.h'
1264!~~~> Input variables     
1265      KPP_REAL T, Y(KPP_NVAR)
1266!~~~> Output variables     
1267      KPP_REAL Jcb(KPP_LU_NONZERO)
1268!~~~> Local variables
1269      KPP_REAL Told     
1270!~~~> Collect statistics     
1271      INTEGER Nfun,Njac,Nstp,Nacc,Nrej,Ndec,Nsol,Nsng
1272      COMMON /Statistics/ Nfun,Njac,Nstp,Nacc,Nrej,
1273     &       Ndec,Nsol,Nsng
1274
1275      Told = TIME
1276      TIME = T   
1277      CALL Update_SUN()
1278      CALL Update_RCONST()
1279      CALL Jac_SP( Y, FIX, RCONST, Jcb )
1280      TIME = Told
1281     
1282      Njac = Njac+1
1283
1284      RETURN
1285      END !  SUBROUTINE JacTemplate                                                                         
1286
Note: See TracBrowser for help on using the repository browser.