?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Original WRF subgrid support version from John Michalakes without fire
[wrffire.git] / wrfv2_fire / chem / KPP / kpp / kpp-2.1 / int / rosenbrock.f
blobf7920b8f3b6e16fe91e85688fa1d22336fbd9655
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
600 101 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
A coupled weather-fire forecasting model built on top of Weather Research and Forecasting (WRF).