source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/templates/rosenbrock_vec.f90.save @ 2696

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

Merge of branch palm4u into trunk

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