source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int.modified_WCOPY/kpp_radau5.f90 @ 3534

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

Merge of branch palm4u into trunk

File size: 37.8 KB
Line 
1!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
2!  RADAU5 - Runge-Kutta method based on Radau-2A quadrature               !
3!                       (2 stages, order 5)                               !
4!  By default the code employs the KPP sparse linear algebra routines     !
5!  Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK)        !
6!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
7
8MODULE KPP_ROOT_Integrator
9
10  USE KPP_ROOT_Precision, ONLY: dp
11  USE KPP_ROOT_Parameters, ONLY: NVAR, LU_NONZERO
12  USE KPP_ROOT_Jacobian, ONLY: LU_DIAG
13  USE KPP_ROOT_LinearAlgebra
14
15  IMPLICIT NONE
16  PUBLIC
17  SAVE
18 
19  ! Statistics
20  INTEGER :: Nfun, Njac, Nstp, Nacc, Nrej, Ndec, Nsol, Nsng
21 
22  ! Method parameters
23  KPP_REAL :: Transf(3,3), TransfInv(3,3),      &
24                   rkA(3,3), rkB(3), rkC(3), rkE(3), &
25                   rkGamma, rkAlpha, rkBeta
26 
27  ! description of the error numbers IERR
28  CHARACTER(LEN=50), PARAMETER, DIMENSION(-11:1) :: IERR_NAMES = (/ &
29    'Matrix is repeatedly singular                     ', & ! -11
30    'Step size too small: T + 10*H = T or H < Roundoff ', & ! -10
31    'No of steps exceeds maximum bound                 ', & ! -9 
32    'Tolerances are too small                          ', & ! -8 
33    'Improper values for Qmin, Qmax                    ', & ! -7 
34    'Newton stopping tolerance too small               ', & ! -6 
35    'Improper value for ThetaMin                       ', & ! -5 
36    'Improper values for FacMin/FacMax/FacSafe/FacRej  ', & ! -4 
37    'Hmin/Hmax/Hstart must be positive                 ', & ! -3 
38    'Improper value for maximal no of Newton iterations', & ! -2 
39    'Improper value for maximal no of steps            ', & ! -1 
40    '                                                  ', & !  0 (not used)
41    'Success                                           ' /) !  1
42
43CONTAINS
44
45  ! **************************************************************************
46
47  SUBROUTINE INTEGRATE( TIN, TOUT, &
48    ICNTRL_U, RCNTRL_U, ISTATUS_U, RSTATUS_U, IERR_U )
49
50    USE KPP_ROOT_Parameters, ONLY: NVAR
51    USE KPP_ROOT_Global,     ONLY: ATOL,RTOL,VAR
52
53    IMPLICIT NONE
54
55    KPP_REAL :: TIN  ! TIN - Start Time
56    KPP_REAL :: TOUT ! TOUT - End Time
57    INTEGER,  INTENT(IN),  OPTIONAL :: ICNTRL_U(20)
58    KPP_REAL, INTENT(IN),  OPTIONAL :: RCNTRL_U(20)
59    INTEGER,  INTENT(OUT), OPTIONAL :: ISTATUS_U(20)
60    KPP_REAL, INTENT(OUT), OPTIONAL :: RSTATUS_U(20)
61    INTEGER,  INTENT(OUT), OPTIONAL :: IERR_U
62
63    KPP_REAL, SAVE :: H
64    INTEGER :: IERR
65
66    KPP_REAL :: RCNTRL(20), RSTATUS(20)
67    INTEGER :: ICNTRL(20), ISTATUS(20)
68    INTEGER, SAVE :: Ntotal = 0
69
70    H =0.0_dp
71
72    !~~~> fine-tune the integrator:
73    ICNTRL(:)  = 0
74    ICNTRL(2)  = 0 ! 0=vector tolerances, 1=scalar tolerances
75    ICNTRL(5)  = 8 ! Max no. of Newton iterations
76    ICNTRL(6)  = 1 ! Starting values for Newton are interpolated (0) or zero (1)
77    ICNTRL(11) = 1 ! Gustaffson (1) or classic(2) controller
78    RCNTRL(1:20) = 0._dp
79
80    !~~~> if optional parameters are given, and if they are >0,
81    !     then use them to overwrite default settings
82    IF (PRESENT(ICNTRL_U)) ICNTRL(1:20) = ICNTRL_U(1:20)
83    IF (PRESENT(RCNTRL_U)) RCNTRL(1:20) = RCNTRL_U(1:20)
84   
85
86    CALL RADAU5(  NVAR,TIN,TOUT,VAR,H,                  &
87                  RTOL,ATOL,                            &
88                  RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR  )
89
90!!$    Ntotal = Ntotal + Nstp
91!!$    PRINT*,'NSTEPS=',Nstp,' (',Ntotal,')'
92
93    Nfun = Nfun + ISTATUS(1)
94    Njac = Njac + ISTATUS(2)
95    Nstp = Nstp + ISTATUS(3)
96    Nacc = Nacc + ISTATUS(4)
97    Nrej = Nrej + ISTATUS(5)
98    Ndec = Ndec + ISTATUS(6)
99    Nsol = Nsol + ISTATUS(7)
100    Nsng = Nsng + ISTATUS(8)
101
102    ! if optional parameters are given for output
103    ! use them to store information in them
104    IF (PRESENT(ISTATUS_U)) THEN
105      ISTATUS_U(:) = 0
106      ISTATUS_U(1) = Nfun ! function calls
107      ISTATUS_U(2) = Njac ! jacobian calls
108      ISTATUS_U(3) = Nstp ! steps
109      ISTATUS_U(4) = Nacc ! accepted steps
110      ISTATUS_U(5) = Nrej ! rejected steps (except at the beginning)
111      ISTATUS_U(6) = Ndec ! LU decompositions
112      ISTATUS_U(7) = Nsol ! forward/backward substitutions
113    ENDIF
114    IF (PRESENT(RSTATUS_U)) THEN
115      RSTATUS_U(:) = 0.
116      RSTATUS_U(1) = TOUT ! final time
117    ENDIF
118    IF (PRESENT(IERR_U)) IERR_U = IERR
119
120! mz_rs_20050716: IERR is returned to the user who then decides what to do
121! about it, i.e. either stop the run or ignore it.
122!!$    IF (IERR < 0) THEN
123!!$      PRINT *,'RADAU: Unsuccessful exit at T=', TIN,' (IERR=',IERR,')'
124!!$      STOP
125!!$    ENDIF
126
127  END SUBROUTINE INTEGRATE
128
129   SUBROUTINE RADAU5(N,T,Tend,Y,H,RelTol,AbsTol,    &
130                       RCNTRL,ICNTRL,RSTATUS,ISTATUS,IDID)
131
132!~~~>-----------------------------------------------
133!     NUMERICAL SOLUTION OF A STIFF (OR DIFFERENTIAL ALGEBRAIC)
134!     SYSTEM OF FirstStep 0RDER ORDINARY DIFFERENTIAL EQUATIONS
135!                     M*Y'=F(T,Y).
136!     THE SYSTEM CAN BE (LINEARLY) IMPLICIT (MASS-MATRIX M  /=  I)
137!     OR EXPLICIT (M=I).
138!     THE METHOD USED IS AN IMPLICIT RUNGE-KUTTA METHOD (RADAU IIA)
139!     OF ORDER 5 WITH STEP SIZE CONTROL AND CONTINUOUS OUTPUT.
140!     C.F. SECTION IV.8
141!
142!     AUTHORS: E. HAIRER AND G. WANNER
143!              UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES
144!              CH-1211 GENEVE 24, SWITZERLAND
145!              E-MAIL:  HAIRER@DIVSUN.UNIGE.CH,  WANNER@DIVSUN.UNIGE.CH
146!
147!     THIS CODE IS PART OF THE BOOK:
148!         E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
149!         EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
150!         SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS 14,
151!         SPRINGER-VERLAG (1991)
152!
153!     VERSION OF SEPTEMBER 30, 1995
154!
155!     INPUT PARAMETERS
156!     ----------------
157!     N           DIMENSION OF THE SYSTEM
158!
159!     FCN         NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE
160!                 VALUE OF F(T,Y):
161!                 SUBROUTINE FCN(N,T,Y,F)
162!                    KPP_REAL T,Y(N),F(N)
163!                    F(1)=...   ETC.
164!                 RPAR, IPAR (SEE BELOW)
165!
166!     T           INITIAL TIME VALUE
167!
168!     Tend        FINAL T-VALUE (Tend-T MAY BE POSITIVE OR NEGATIVE)
169!
170!     Y(N)        INITIAL VALUES FOR Y
171!
172!     H           INITIALL STEP SIZE GUESS;
173!                 FOR STIFF EQUATIONS WITH INITIALL TRANSIENT,
174!                 H=1.D0/(NORM OF F'), USUALLY 1.D-3 OR 1.D-5, IS GOOD.
175!                 THIS CHOICE IS NOT VERY IMPORTANT, THE STEP SIZE IS
176!                 QUICKLY ADAPTED. (IF H=0.D0, THE CODE PUTS H=1.D-6).
177!
178!     RelTol,AbsTol   RELATIVE AND ABSOLUTE ERROR TOLERANCES.
179!          for ICNTRL(2) = 0: AbsTol, RelTol are N-dimensional vectors
180!                        = 1: AbsTol, RelTol are scalars
181!
182!          -----  CONTINUOUS OUTPUT: -----
183!                 DURING CALLS TO "SOLOUT", A CONTINUOUS SOLUTION
184!                 FOR THE INTERVAL [Told,T] IS AVAILABLE THROUGH
185!                 THE FUNCTION
186!                        >>>   CONTR5(I,S,CONT,LRC)   <<<
187!                 WHICH PROVIDES AN APPROXIMATION TO THE I-TH
188!                 COMPONENT OF THE SOLUTION AT THE POINT S. THE VALUE
189!                 S SHOULD LIE IN THE INTERVAL [Told,T].
190!                 DO NOT CHANGE THE ENTRIES OF CONT(LRC), IF THE
191!                 DENSE OUTPUT FUNCTION IS USED.
192!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
193!
194!~~~>     INPUT PARAMETERS:
195!
196!    Note: For input parameters equal to zero the default values of the
197!          corresponding variables are used.
198!
199!~~~>  Integer input parameters:
200
201!    ICNTRL(1) = not used
202!
203!    ICNTRL(2) = 0: AbsTol, RelTol are NVAR-dimensional vectors
204!              = 1: AbsTol, RelTol are scalars
205!
206!    ICNTRL(3) = not used
207!
208!    ICNTRL(4)  -> maximum number of integration steps
209!        For ICNTRL(4)=0 the default value of 10000 is used
210!
211!    ICNTRL(5)  -> maximum number of Newton iterations
212!        For ICNTRL(5)=0 the default value of 8 is used
213!
214!    ICNTRL(6)  -> starting values of Newton iterations:
215!        ICNTRL(6)=0 : starting values are obtained from
216!                      the extrapolated collocation solution
217!                      (the default)
218!        ICNTRL(6)=1 : starting values are zero
219!
220!    ICNTRL(11) -> switch for step size strategy
221!              ICNTRL(8) == 1:  mod. predictive controller (Gustafsson)
222!              ICNTRL(8) == 2:  classical step size control
223!              the default value (for iwork(8)=0) is iwork(8)=1.
224!              the choice iwork(8) == 1 seems to produce safer results;
225!              for simple problems, the choice iwork(8) == 2 produces
226!              often slightly faster runs
227!              ( currently unused )
228!
229!~~~>  Real input parameters:
230!
231!    RCNTRL(1)  -> not used
232!
233!    RCNTRL(2)  -> Hmax, upper bound for the integration step size
234!
235!    RCNTRL(3)  -> not used
236!
237!    RCNTRL(4)  -> FacMin, lower bound on step decrease factor (default=0.2)
238!
239!    RCNTRL(5)  -> FacMax, upper bound on step increase factor (default=6)
240!
241!    RCNTRL(6)  -> FacRej, step decrease factor after multiple rejections
242!                 (default=0.1)
243!
244!    RCNTRL(7)  -> FacSafe, by which the new step is slightly smaller
245!                  than the predicted value  (default=0.9)
246!
247!    RCNTRL(8)  -> ThetaMin. If Newton convergence rate smaller
248!                  than ThetaMin the Jacobian is not recomputed;
249!                  (default=0.001)
250!
251!    RCNTRL(9)  -> NewtonTol, stopping criterion for Newton's method
252!                  (default=0.03)
253!
254!    RCNTRL(10) -> Qmin
255!
256!    RCNTRL(11) -> Qmax. If Qmin < Hnew/Hold < Qmax, then the
257!                  step size is kept constant and the LU factorization
258!                  reused (default Qmin=1, Qmax=1.2)
259!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
260!
261!
262!    OUTPUT PARAMETERS
263!    -----------------
264!    T           T-VALUE FOR WHICH THE SOLUTION HAS BEEN COMPUTED
265!                     (AFTER SUCCESSFUL RETURN T=Tend).
266!
267!    Y(N)        NUMERICAL SOLUTION AT T
268!
269!    H           PREDICTED STEP SIZE OF THE LastStep ACCEPTED STEP
270!
271!    IDID        REPORTS ON SUCCESSFULNESS UPON RETURN:
272!
273!    ISTATUS(1) = No. of function calls
274!    ISTATUS(2) = No. of jacobian calls
275!    ISTATUS(3) = No. of steps
276!    ISTATUS(4) = No. of accepted steps
277!    ISTATUS(5) = No. of rejected steps (except at the beginning)
278!    ISTATUS(6) = No. of LU decompositions
279!    ISTATUS(7) = No. of forward/backward substitutions
280!    ISTATUS(8) = No. of singular matrix decompositions
281!
282!    RSTATUS(1)  -> Texit, the time corresponding to the
283!                   computed Y upon return
284!    RSTATUS(2)  -> Hexit, last accepted step before exit
285!                   For multiple restarts, use Hexit as Hstart
286!                   in the subsequent run
287!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
288
289      IMPLICIT NONE
290     
291      INTEGER :: N
292      KPP_REAL :: Y(N),AbsTol(N),RelTol(N),RCNTRL(20),RSTATUS(20)
293      INTEGER :: ICNTRL(20), ISTATUS(20)
294      LOGICAL :: StartNewton, Gustafsson
295      INTEGER :: IDID, ITOL
296      KPP_REAL :: H,Tend,T
297
298      !~~~> Control arguments
299      INTEGER :: Max_no_steps, NewtonMaxit
300      KPP_REAL :: Hstart,Hmin,Hmax,Qmin,Qmax
301      KPP_REAL :: Roundoff, ThetaMin,TolNewton
302      KPP_REAL :: FacSafe,FacMin,FacMax,FacRej
303      !~~~> Local variables
304      INTEGER :: i
305      KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0
306 
307      !~~~> variables from the former COMMON block /CONRA5/
308      ! KPP_REAL :: Tsol, Hsol
309   
310!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
311!        SETTING THE PARAMETERS
312!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
313       Nfun=0
314       Njac=0
315       Nstp=0
316       Nacc=0
317       Nrej=0
318       Ndec=0
319       Nsol=0
320       IDID = 0
321       
322!~~~> ICNTRL(1) - autonomous system - not used       
323!~~~> ITOL: 1 for vector and 0 for scalar AbsTol/RelTol
324      IF (ICNTRL(2) == 0) THEN
325         ITOL = 1
326      ELSE
327         ITOL = 0
328      END IF
329!~~~> ICNTRL(3) - method selection - not used       
330!~~~> Max_no_steps: the maximal number of time steps
331      IF (ICNTRL(4) == 0) THEN
332         Max_no_steps = 10000
333      ELSE
334         Max_no_steps=ICNTRL(4)
335         IF (Max_no_steps <= 0) THEN
336            WRITE(6,*) 'ICNTRL(4)=',ICNTRL(4)
337            CALL RAD_ErrorMsg(-1,T,ZERO,IDID)
338         END IF
339      END IF
340!~~~> NewtonMaxit    MAXIMAL NUMBER OF NEWTON ITERATIONS
341      IF (ICNTRL(5) == 0) THEN
342         NewtonMaxit = 8
343      ELSE
344         NewtonMaxit=ICNTRL(5)
345         IF (NewtonMaxit <= 0) THEN
346            WRITE(6,*) 'ICNTRL(5)=',ICNTRL(5)
347            CALL RAD_ErrorMsg(-2,T,ZERO,IDID)
348          END IF
349      END IF
350!~~~> StartNewton:  Use extrapolation for starting values of Newton iterations
351      IF (ICNTRL(6) == 0) THEN
352         StartNewton = .TRUE.
353      ELSE
354         StartNewton = .FALSE.
355      END IF
356!~~~> Gustafsson: step size controller
357      IF(ICNTRL(11) == 0)THEN
358         Gustafsson=.TRUE.
359      ELSE
360         Gustafsson=.FALSE.
361      END IF
362
363
364!~~~> Roundoff   SMALLEST NUMBER SATISFYING 1.0d0+Roundoff>1.0d0
365      Roundoff=WLAMCH('E');
366
367!~~~> RCNTRL(1) = Hmin - not used
368      Hmin = ZERO
369!~~~> Hmax = maximal step size
370      IF (RCNTRL(2) == ZERO) THEN
371         Hmax=Tend-T
372      ELSE
373         Hmax=MAX(ABS(RCNTRL(7)),ABS(Tend-T))
374      END IF
375!~~~> RCNTRL(3) = Hstart - not used
376      Hstart = ZERO
377!~~~> FacMin: lower bound on step decrease factor
378      IF(RCNTRL(4) == ZERO)THEN
379         FacMin = 0.2d0
380      ELSE
381         FacMin = RCNTRL(4)
382      END IF
383!~~~> FacMax: upper bound on step increase factor
384      IF(RCNTRL(5) == ZERO)THEN
385         FacMax = 8.D0
386      ELSE
387         FacMax = RCNTRL(5)
388      END IF
389!~~~> FacRej: step decrease factor after 2 consecutive rejections
390      IF(RCNTRL(6) == ZERO)THEN
391         FacRej = 0.1d0
392      ELSE
393         FacRej = RCNTRL(6)
394      END IF
395!~~~> FacSafe:  by which the new step is slightly smaller
396!               than the predicted value
397      IF (RCNTRL(7) == ZERO) THEN
398         FacSafe=0.9d0
399      ELSE
400         FacSafe=RCNTRL(7)
401      END IF
402      IF ( (FacMax < ONE) .OR. (FacMin > ONE) .OR. &
403           (FacSafe <= 0.001D0) .OR. (FacSafe >= 1.0d0) ) THEN
404            WRITE(6,*)'RCNTRL(4:7)=',RCNTRL(4:7)
405            CALL RAD_ErrorMsg(-4,T,ZERO,IDID)
406      END IF
407
408!~~~> ThetaMin:  decides whether the Jacobian should be recomputed
409      IF (RCNTRL(8) == ZERO) THEN
410         ThetaMin = 1.0d-3
411      ELSE
412         ThetaMin=RCNTRL(8)
413         IF (ThetaMin <= 0.0d0 .OR. ThetaMin >= 1.0d0) THEN
414            WRITE(6,*) 'RCNTRL(8)=', RCNTRL(8)
415            CALL RAD_ErrorMsg(-5,T,ZERO,IDID)
416         END IF
417      END IF
418!~~~> TolNewton:  stopping crierion for Newton's method
419      IF (RCNTRL(9) == ZERO) THEN
420         TolNewton = 3.0d-2
421      ELSE
422         TolNewton = RCNTRL(9)
423         IF (TolNewton <= Roundoff) THEN
424            WRITE(6,*) 'RCNTRL(9)=',RCNTRL(9)
425            CALL RAD_ErrorMsg(-6,T,ZERO,IDID)
426         END IF
427      END IF
428!~~~> Qmin AND Qmax: IF Qmin < Hnew/Hold < Qmax, STEP SIZE = CONST.
429      IF (RCNTRL(10) == ZERO) THEN
430         Qmin=1.D0
431      ELSE
432         Qmin=RCNTRL(10)
433      END IF
434      IF (RCNTRL(11) == ZERO) THEN
435         Qmax=1.2D0
436      ELSE
437         Qmax=RCNTRL(11)
438      END IF
439      IF (Qmin > ONE .OR. Qmax < ONE) THEN
440         WRITE(6,*) 'RCNTRL(10:11)=',Qmin,Qmax
441         CALL RAD_ErrorMsg(-7,T,ZERO,IDID)
442      END IF
443!~~~> Check if tolerances are reasonable
444      IF (ITOL == 0) THEN
445          IF (AbsTol(1) <= ZERO.OR.RelTol(1) <= 10.d0*Roundoff) THEN
446              WRITE (6,*) 'AbsTol/RelTol=',AbsTol,RelTol 
447              CALL RAD_ErrorMsg(-8,T,ZERO,IDID)
448          END IF
449      ELSE
450          DO i=1,N
451          IF (AbsTol(i) <= ZERO.OR.RelTol(i) <= 10.d0*Roundoff) THEN
452              WRITE (6,*) 'AbsTol/RelTol(',i,')=',AbsTol(i),RelTol(i)
453              CALL RAD_ErrorMsg(-8,T,ZERO,IDID)
454          END IF
455          END DO
456      END IF
457
458!~~~> WHEN A FAIL HAS OCCURED, RETURN
459      IF (IDID < 0) RETURN
460
461     
462!~~~>  CALL TO CORE INTEGRATOR ------------
463      CALL RAD_Integrator( N,T,Y,Tend,Hmax,H,RelTol,AbsTol,ITOL,IDID, &
464        Max_no_steps,Roundoff,FacSafe,ThetaMin,TolNewton,Qmin,Qmax,   &
465        NewtonMaxit,StartNewton,Gustafsson,FacMin,FacMax,FacRej )
466       
467      ISTATUS(1)=Nfun
468      ISTATUS(2)=Njac
469      ISTATUS(3)=Nstp
470      ISTATUS(4)=Nacc
471      ISTATUS(5)=Nrej
472      ISTATUS(6)=Ndec
473      ISTATUS(7)=Nsol
474      ISTATUS(8)=Nsng
475
476   ! END SUBROUTINE RADAU5
477   CONTAINS ! INTERNAL PROCEDURES TO RADAU5
478
479
480!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
481   SUBROUTINE RAD_Integrator( N,T,Y,Tend,Hmax,H,RelTol,AbsTol,ITOL,IDID, &
482        Max_no_steps,Roundoff,FacSafe,ThetaMin,TolNewton,Qmin,Qmax,      &
483        NewtonMaxit,StartNewton,Gustafsson,FacMin,FacMax,FacRej )
484!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
485!     CORE INTEGRATOR FOR RADAU5
486!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
487
488      IMPLICIT NONE
489      INTEGER :: N
490      KPP_REAL  Y(NVAR),Z1(NVAR),Z2(NVAR),Z3(NVAR),Y0(NVAR),&
491                     SCAL(NVAR),F1(NVAR),F2(NVAR),F3(NVAR),      &
492                     CONT(N,4),AbsTol(NVAR),RelTol(NVAR)
493
494#ifdef FULL_ALGEBRA
495      KPP_REAL  :: FJAC(NVAR,NVAR), E1(NVAR,NVAR)
496      DOUBLE COMPLEX :: E2(NVAR,NVAR)   
497#else
498      KPP_REAL  :: FJAC(LU_NONZERO), E1(LU_NONZERO)
499      DOUBLE COMPLEX :: E2(LU_NONZERO)   
500#endif
501                 
502      !~~~> Local variables
503      KPP_REAL  :: TMP(NVAR), T, Tend, Tdirection, &
504                 H, Hmax, HmaxN, Hacc, Hnew, Hopt, Hold, &
505                 Fac, FacMin, Facmax, FacSafe, FacRej, FacGus, FacConv, &
506                 Theta, ThetaMin, TolNewton, ERR, ERRACC,   &
507                 Qmin, Qmax, DYNO, Roundoff, &
508                 AK, AK1, AK2, AK3, C3Q, &
509                 Qnewton, DYTH, THQ, THQOLD, DYNOLD, &
510                 DENOM, C1Q, C2Q, ALPHA, BETA, GAMMA, CFAC, ACONT3, QT
511      INTEGER :: IP1(NVAR),IP2(NVAR), ITOL, IDID, Max_no_steps, &
512                 NewtonIter, NewtonMaxit, ISING
513      LOGICAL :: REJECT, FirstStep, FreshJac, LastStep, &
514                 Gustafsson, StartNewton, NewtonDone
515      ! KPP_REAL, PARAMETER  :: ONE = 1.0d0, ZERO = 0.0d0
516     
517!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
518!  INITIALISATIONS
519!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
520
521      CALL RAD_Coefficients
522           
523      Tdirection=SIGN(1.D0,Tend-T)
524      HmaxN=MIN(ABS(Hmax),ABS(Tend-T))
525      H=MIN(ABS(Hmin),ABS(Hstart))
526      H=MIN(ABS(H),HmaxN)
527      IF (ABS(H) <= 10.D0*Roundoff) H=1.0D-6
528      H=SIGN(H,Tdirection)
529      Hold=H
530      REJECT=.FALSE.
531      FirstStep=.TRUE.
532      LastStep=.FALSE.
533      FreshJac=.FALSE.; Theta=1.0d0
534      IF ((T+H*1.0001D0-Tend)*Tdirection >= 0.D0) THEN
535         H=Tend-T
536         LastStep=.TRUE.
537      END IF
538      FacConv=1.D0
539      CFAC=FacSafe*(1+2*NewtonMaxit)
540      Nsng=0
541!      Told=T
542      CALL RAD_ErrorScale(N,ITOL,AbsTol,RelTol,Y,SCAL)
543      CALL FUN_CHEM(T,Y,Y0)
544     
545!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
546!~~~>  Time loop begins
547!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
548Tloop: DO WHILE ( (Tend-T)*Tdirection - Roundoff > ZERO )
549
550      !~~~> COMPUTE JACOBIAN MATRIX ANALYTICALLY
551      IF ( (.NOT.FreshJac)  .AND. (Theta > ThetaMin) ) THEN
552        CALL JAC_CHEM(T,Y,FJAC)
553        FreshJac=.TRUE.
554      END IF
555     
556      !~~~> Compute the matrices E1 and E2 and their decompositions
557      GAMMA  = rkGamma/H
558      ALPHA = rkAlpha/H
559      BETA  = rkBeta/H
560      CALL RAD_DecompReal(N,FJAC,GAMMA,E1,IP1,ISING)
561      IF (ISING /= 0) THEN
562          Nsng=Nsng+1
563          IF (Nsng >= 5) THEN
564            CALL RAD_ErrorMsg(-12,T,H,IDID); RETURN
565          END IF
566          H=H*0.5D0; REJECT=.TRUE.; LastStep=.FALSE.
567          CYCLE Tloop
568      END IF   
569      CALL RAD_DecompCmplx(N,FJAC,ALPHA,BETA,E2,IP2,ISING)
570      IF (ISING /= 0) THEN
571          Nsng=Nsng+1
572          IF (Nsng >= 5) THEN
573            CALL RAD_ErrorMsg(-12,T,H,IDID); RETURN
574          END IF
575          H=H*0.5D0; REJECT=.TRUE.; LastStep=.FALSE.
576          CYCLE Tloop
577      END IF
578
579   30 CONTINUE
580      Nstp=Nstp+1
581      IF (Nstp > Max_no_steps) THEN
582        PRINT*,'Max number of time steps is ',Max_no_steps
583        CALL RAD_ErrorMsg(-9,T,H,IDID); RETURN
584      END IF
585      IF (0.1D0*ABS(H) <= ABS(T)*Roundoff)  THEN
586        CALL RAD_ErrorMsg(-10,T,H,IDID); RETURN
587      END IF
588 
589!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
590!  STARTING VALUES FOR NEWTON ITERATION
591!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
592      IF ( FirstStep .OR. (.NOT.StartNewton) ) THEN
593         CALL Set2zero(N,Z1)
594         CALL Set2zero(N,Z2)
595         CALL Set2zero(N,Z3)
596         CALL Set2zero(N,F1)
597         CALL Set2zero(N,F2)
598         CALL Set2zero(N,F3)
599      ELSE
600         C3Q=H/Hold
601         C1Q=rkC(1)*C3Q
602         C2Q=rkC(2)*C3Q
603         DO i=1,N
604            AK1=CONT(i,2)
605            AK2=CONT(i,3)
606            AK3=CONT(i,4)
607            Z1(i)=C1Q*(AK1+(C1Q-rkC(2)+ONE)*(AK2+(C1Q-rkC(1)+ONE)*AK3))
608            Z2(i)=C2Q*(AK1+(C2Q-rkC(2)+ONE)*(AK2+(C2Q-rkC(1)+ONE)*AK3))
609            Z3(i)=C3Q*(AK1+(C3Q-rkC(2)+ONE)*(AK2+(C3Q-rkC(1)+ONE)*AK3))
610         END DO
611         !  F(1,2,3) = TransfInv x Z(1,2,3)
612         CALL RAD_Transform(N,TransfInv,Z1,Z2,Z3,F1,F2,F3)
613      END IF
614!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
615!  LOOP FOR THE SIMPLIFIED NEWTON ITERATION
616!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
617           
618            FacConv = MAX(FacConv,Roundoff)**0.8D0
619            Theta=ABS(ThetaMin)
620
621NewtonLoop:DO  NewtonIter = 1, NewtonMaxit 
622 
623            !~~~>  The right-hand side
624            DO i=1,N
625               TMP(i)=Y(i)+Z1(i)
626            END DO
627            CALL FUN_CHEM(T+rkC(1)*H,TMP,Z1)
628            DO i=1,N
629               TMP(i)=Y(i)+Z2(i)
630            END DO
631            CALL FUN_CHEM(T+rkC(2)*H,TMP,Z2)
632            DO i=1,N
633               TMP(i)=Y(i)+Z3(i)
634            END DO
635            CALL FUN_CHEM(T+rkC(3)*H,TMP,Z3)
636
637            !~~~> Solve the linear systems
638            !  Z(1,2,3) = TransfInv x Z(1,2,3)
639            CALL RAD_Transform(N,TransfInv,Z1,Z2,Z3,Z1,Z2,Z3)
640            CALL RAD_Solve( N,FJAC,GAMMA,ALPHA,BETA,E1,E2, &
641                        Z1,Z2,Z3,F1,F2,F3,CONT,IP1,IP2,ISING )
642            Nsol=Nsol+1
643           
644            DYNO=0.0d0
645            DO i=1,N
646               DENOM=SCAL(i)
647               DYNO=DYNO+(Z1(i)/DENOM)**2+(Z2(i)/DENOM)**2+(Z3(i)/DENOM)**2
648            END DO
649            DYNO=SQRT(DYNO/(3*N))
650           
651            !~~~> Bad convergence or number of iterations too large
652            IF ( (NewtonIter > 1) .AND. (NewtonIter < NewtonMaxit) ) THEN
653                THQ=DYNO/DYNOLD
654                IF (NewtonIter == 2) THEN
655                   Theta=THQ
656                ELSE
657                   Theta=SQRT(THQ*THQOLD)
658                END IF
659                THQOLD=THQ
660                IF (Theta < 0.99d0) THEN
661                    FacConv=Theta/(1.0d0-Theta)
662                    DYTH=FacConv*DYNO*Theta**(NewtonMaxit-1-NewtonIter)/TolNewton
663                    IF (DYTH >= 1.0d0) THEN
664                         Qnewton=MAX(1.0D-4,MIN(20.0d0,DYTH))
665                         FAC=.8D0*Qnewton**(-1.0d0/(4.0d0+NewtonMaxit-1-NewtonIter))
666                         H=FAC*H
667                         REJECT=.TRUE.
668                         LastStep=.FALSE.
669                         CYCLE Tloop
670                    END IF
671                ELSE ! Non-convergence of Newton
672                   H=H*0.5D0; REJECT=.TRUE.; LastStep=.FALSE.
673                   CYCLE Tloop
674                END IF
675            END IF
676            DYNOLD=MAX(DYNO,Roundoff) 
677            CALL WAXPY(N,ONE,Z1,1,F1,1) ! F1 <- F1 + Z1
678            CALL WAXPY(N,ONE,Z2,1,F2,1) ! F2 <- F2 + Z2
679            CALL WAXPY(N,ONE,Z3,1,F3,1) ! F3 <- F3 + Z3
680            !  Z(1,2,3) = Transf x F(1,2,3)
681            CALL RAD_Transform(N,Transf,F1,F2,F3,Z1,Z2,Z3)
682            NewtonDone = (FacConv*DYNO <= TolNewton)
683            IF (NewtonDone) EXIT NewtonLoop
684           
685       END DO NewtonLoop
686           
687       IF (.NOT.NewtonDone) THEN
688           CALL RAD_ErrorMsg(-8,T,H,IDID);
689           H=H*0.5D0; REJECT=.TRUE.; LastStep=.FALSE.
690           CYCLE Tloop
691       END IF
692
693           
694!~~~> ERROR ESTIMATION
695      CALL  RAD_ErrorEstimate(N,FJAC,H,Y0,Y,T, &
696               E1,Z1,Z2,Z3,IP1,SCAL,ERR,       &
697               FirstStep,REJECT,GAMMA)
698!~~~> COMPUTATION OF Hnew
699      Fac  = ERR**(-0.25d0)*   &
700             MIN(FacSafe,(NewtonIter+2*NewtonMaxit)/CFAC)
701      Fac  = MIN(FacMax,MAX(FacMin,Fac))
702      Hnew = Fac*H
703
704!~~~> IS THE ERROR SMALL ENOUGH ?
705accept:IF (ERR < ONE) THEN !~~~> STEP IS ACCEPTED
706         FirstStep=.FALSE.
707         Nacc=Nacc+1
708         IF (Gustafsson) THEN
709            !~~~> Predictive controller of Gustafsson
710            !~~~> Currently not implemented
711            IF (Nacc > 1) THEN
712               FacGus=FacSafe*(H/Hacc)*(ERR**2/ERRACC)**(-0.25d0)
713               FacGus=MIN(FacMax,MAX(FacMin,FacGus))
714               Fac=MIN(Fac,FacGus)
715               Hnew=H*Fac
716            END IF
717            Hacc=H
718            ERRACC=MAX(1.0D-2,ERR)
719         END IF
720         ! Told = T
721         Hold = H
722         T=T+H
723         DO i=1,N
724            Y(i)=Y(i)+Z3(i)
725            CONT(i,2)=(Z2(i)-Z3(i))/(rkC(2)-ONE)
726            AK=(Z1(i)-Z2(i))/(rkC(1)-rkC(2))
727            ACONT3=Z1(i)/rkC(1)
728            ACONT3=(AK-ACONT3)/rkC(2)
729            CONT(i,3)=(AK-CONT(i,2))/(rkC(1)-ONE)
730            CONT(i,4)=CONT(i,3)-ACONT3
731         END DO
732         CALL RAD_ErrorScale(N,ITOL,AbsTol,RelTol,Y,SCAL)
733         FreshJac=.FALSE.
734         IF (LastStep) THEN
735            H=Hopt
736            IDID=1
737            RETURN
738         END IF
739         CALL FUN_CHEM(T,Y,Y0)
740         Hnew=Tdirection*MIN(ABS(Hnew),HmaxN)
741         Hopt=Hnew
742         Hopt=MIN(H,Hnew)
743         IF (REJECT) Hnew=Tdirection*MIN(ABS(Hnew),ABS(H))
744         REJECT=.FALSE.
745         IF ((T+Hnew/Qmin-Tend)*Tdirection >= 0.D0) THEN
746            H=Tend-T
747            LastStep=.TRUE.
748         ELSE
749            QT=Hnew/H
750            IF ( (Theta<=ThetaMin) .AND. (QT>=Qmin) &
751                            .AND. (QT<=Qmax) ) GOTO 30
752            H=Hnew
753         END IF
754         CYCLE Tloop
755      ELSE accept !~~~> STEP IS REJECTED
756         REJECT=.TRUE.
757         LastStep=.FALSE.
758         IF (FirstStep) THEN
759             H=H*FacRej
760         ELSE
761             H=Hnew
762         END IF
763         IF (Nacc >= 1) Nrej=Nrej+1
764         CYCLE Tloop
765      END IF accept
766     
767     
768      END DO Tloop
769     
770
771   END SUBROUTINE RAD_Integrator
772
773!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
774 SUBROUTINE RAD_ErrorMsg(Code,T,H,IERR)
775!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
776!    Handles all error messages
777!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
778
779   IMPLICIT NONE
780   KPP_REAL, INTENT(IN) :: T, H
781   INTEGER, INTENT(IN)  :: Code
782   INTEGER, INTENT(OUT) :: IERR
783
784   IERR = Code
785   PRINT * , &
786     'Forced exit from RADAU5 due to the following error:'
787   IF ((Code>=-11).AND.(Code<=-1)) THEN
788     PRINT *, IERR_NAMES(Code)
789   ELSE
790     PRINT *, 'Unknown Error code: ', Code
791   ENDIF
792
793   PRINT *, "T=", T, "and H=", H
794
795 END SUBROUTINE RAD_ErrorMsg
796
797
798!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
799 SUBROUTINE RAD_ErrorScale(N,ITOL,AbsTol,RelTol,Y,SCAL)
800!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
801!    Handles all error messages
802!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
803   IMPLICIT NONE
804   INTEGER, INTENT(IN)  :: N, ITOL
805   KPP_REAL, INTENT(IN) :: AbsTol(*), RelTol(*), Y(N)
806   KPP_REAL, INTENT(OUT) :: SCAL(N)
807   INTEGER :: i
808   
809   IF (ITOL==0) THEN
810       DO i=1,N
811          SCAL(i)=AbsTol(1)+RelTol(1)*ABS(Y(i))
812       END DO
813   ELSE
814       DO i=1,N
815          SCAL(i)=AbsTol(i)+RelTol(i)*ABS(Y(i))
816       END DO
817   END IF
818     
819 END SUBROUTINE RAD_ErrorScale
820 
821!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
822  SUBROUTINE RAD_Transform(N,Tr,Z1,Z2,Z3,F1,F2,F3)
823!~~~>                 F = Tr x Z
824!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
825      IMPLICIT NONE
826      INTEGER :: N, i
827      KPP_REAL :: Tr(3,3),Z1(N),Z2(N),Z3(N),F1(N),F2(N),F3(N)
828      KPP_REAL :: x1, x2, x3
829      DO i=1,N
830          x1 = Z1(i); x2 = Z2(i); x3 = Z3(i)
831          F1(i) = Tr(1,1)*x1 + Tr(1,2)*x2 + Tr(1,3)*x3
832          F2(i) = Tr(2,1)*x1 + Tr(2,2)*x2 + Tr(2,3)*x3
833          F3(i) = Tr(3,1)*x1 + Tr(3,2)*x2 + Tr(3,3)*x3
834      END DO
835  END SUBROUTINE RAD_Transform
836 
837!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
838  SUBROUTINE RAD_Coefficients
839!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
840      IMPLICIT NONE
841      KPP_REAL :: s2,s3,s6,x1,x2,x3,x4,y1,y2,y3,y4,y5
842     
843      s2 = SQRT(2.0d0);
844      s3 = SQRT(3.0d0);
845      s6 = SQRT(6.0d0);
846      x1 = 3.d0**(1.d0/3.d0);
847      x2 = 3.d0**(2.d0/3.d0);
848      x3 = 3.d0**(1.d0/6.d0);
849      x4 = 3.d0**(5.d0/6.d0);
850     
851      rkA(1,1) =  11.d0/45.d0-7.d0/360.d0*s6
852      rkA(1,2) =  37.d0/225.d0-169.d0/1800.d0*s6
853      rkA(1,3) =  -2.d0/225.d0+s6/75
854      rkA(2,1) =  37.d0/225.d0+169.d0/1800.d0*s6
855      rkA(2,2) =  11.d0/45.d0+7.d0/360.d0*s6
856      rkA(2,3) =  -2.d0/225.d0-s6/75
857      rkA(3,1) =  4.d0/9.d0-s6/36
858      rkA(3,2) =  4.d0/9.d0+s6/36
859      rkA(3,3) =  1.d0/9.d0
860     
861      rkB(1)   =  4.d0/9.d0-s6/36
862      rkB(2)   =  4.d0/9.d0+s6/36
863      rkB(3)   =  1.d0/9.d0
864     
865      rkC(1)   =  2.d0/5.d0-s6/10
866      rkC(2)   =  2.d0/5.d0+s6/10
867      rkC(3)   =  1.d0
868
869      ! Error estimation
870      rkE(1) = -(13.d0+7.d0*s6)/3.d0
871      rkE(2) =  (-13.d0+7.d0*s6)/3.d0
872      rkE(3) = -1.d0/3.d0
873
874      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
875      !~~~> Diagonalize the RK matrix:
876      ! TransfInv * inv(rkA) * Transf =
877      !       |  rkGamma      0           0       |
878      !       |      0      rkAlpha   -rkBeta   |
879      !       |      0      rkBeta     rkAlpha  |
880      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
881      rkGamma =  3-x1+x2
882      rkAlpha =  x1/2-x2/2+3
883      rkBeta  =  -x4/2-3.d0/2.d0*x3
884
885      y1 = 36.d0/625.d0*s6
886      y2 = 129.d0/2500.d0*x1
887      y3 = 111.d0/2500.d0*x3*s2
888      Transf(1,1) = -31.d0/1250.d0*s6*x1+37.d0/1250.d0*s6*x2-y1 &
889                    +129.d0/1250.d0*x1-33.d0/1250.d0*x2+49.d0/625.d0
890      Transf(1,2) =  -y1-y2-y3 &
891                    +31.d0/2500.d0*x4*s2+33.d0/2500.d0*x2+49.d0/625.d0
892      Transf(1,3) =  3.d0/2500.d0*x3*(-33-43*x2+31*x3*s2+37*s3*s2)
893      Transf(2,1) =  31.d0/1250.d0*s6*x1-37.d0/1250.d0*s6*x2+y1 &
894                    +129.d0/1250.d0*x1-33.d0/1250.d0*x2+49.d0/625.d0
895      Transf(2,2) =  y1-y2+y3&
896                    -31.d0/2500.d0*x4*s2+33.d0/2500.d0*x2+49.d0/625.d0
897      Transf(2,3) =  -3.d0/2500.d0*x3*(33+43*x2+31*x3*s2+37*s3*s2)
898      Transf(3,1) =  1.d0
899      Transf(3,2) =  1.d0
900      Transf(3,3) =  0.d0
901     
902      y1 = 11.d0/36.d0*x3*s2 + 43.d0/108.d0*x4*s2
903      y2 = 11.d0/36.d0*s2*x2 - 43.d0/36.d0*s2*x1
904      y3 = 31.d0/54.d0*x1 + 37.d0/54.d0*x2
905      y4 = 31.d0/54.d0*x4-37.d0/18.d0*x3
906      y5 = -x2/27+5.d0/27.d0*x1
907      TransfInv(1,1) =  y1 + y3
908      TransfInv(1,2) = -y1 + y3
909      TransfInv(1,3) =  y5 + 1.d0/3.d0
910      TransfInv(2,1) = -y1 - y3
911      TransfInv(2,2) =  y1 - y3
912      TransfInv(2,3) = -y5 + 2.d0/3.d0
913      TransfInv(3,1) =  y4 - y2
914      TransfInv(3,2) =  y4 + y2
915      TransfInv(3,3) =  x3/9+5.d0/27.d0*x4
916     
917   END SUBROUTINE RAD_Coefficients
918
919
920!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
921   SUBROUTINE RAD_DecompReal(N,FJAC,GAMMA,E1,IP1,ISING)
922!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
923      IMPLICIT NONE
924     
925      INTEGER :: N, ISING
926      KPP_REAL :: GAMMA
927#ifdef FULL_ALGEBRA     
928      KPP_REAL :: FJAC(NVAR,NVAR),E1(NVAR,NVAR)
929#else     
930      KPP_REAL :: FJAC(LU_NONZERO),E1(LU_NONZERO)
931#endif     
932      INTEGER :: IP1(N), i, j
933
934#ifdef FULL_ALGEBRA     
935      DO j=1,N
936         DO  i=1,N
937            E1(i,j)=-FJAC(i,j)
938         END DO
939         E1(j,j)=E1(j,j)+GAMMA
940      END DO
941      CALL DGETRF(N,N,E1,N,IP1,ISING) 
942#else     
943      DO i=1,LU_NONZERO
944         E1(i)=-FJAC(i)
945      END DO
946      DO i=1,NVAR
947         j = LU_DIAG(i); E1(j)=E1(j)+GAMMA
948      END DO
949      CALL KppDecomp(E1,ISING)
950#endif     
951     
952   END SUBROUTINE RAD_DecompReal
953
954
955!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
956   SUBROUTINE RAD_DecompCmplx(N,FJAC,ALPHA,BETA,E2,IP2,ISING)
957!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
958      IMPLICIT NONE
959      INTEGER :: N, ISING
960#ifdef FULL_ALGEBRA     
961      KPP_REAL ::  FJAC(N,N)
962      DOUBLE COMPLEX :: E2(N,N)
963#else     
964      KPP_REAL ::  FJAC(LU_NONZERO)
965      DOUBLE COMPLEX :: E2(LU_NONZERO)
966#endif     
967      KPP_REAL :: ALPHA, BETA
968      INTEGER :: IP2(N), i, j
969     
970#ifdef FULL_ALGEBRA     
971      DO j=1,N
972        DO i=1,N
973          E2(i,j) = DCMPLX( -FJAC(i,j), 0.0d0 )
974        END DO
975        E2(j,j) = E2(j,j) + DCMPLX( ALPHA, BETA )
976      END DO
977      CALL ZGETRF(N,N,E2,N,IP2,ISING)     
978#else 
979      DO i=1,LU_NONZERO
980         E2(i) = DCMPLX( -FJAC(i), 0.0d0 )
981      END DO
982      DO i=1,NVAR
983         j = LU_DIAG(i); E2(j)=E2(j)+DCMPLX( ALPHA, BETA )
984      END DO
985      CALL KppDecompCmplx(E2,ISING)   
986#endif     
987      Ndec=Ndec+1
988     
989   END SUBROUTINE RAD_DecompCmplx
990
991
992!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
993   SUBROUTINE RAD_Solve(N,FJAC,GAMMA,ALPHA,BETA,E1,E2,&
994               Z1,Z2,Z3,F1,F2,F3,CONT,IP1,IP2,ISING)
995!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
996      IMPLICIT NONE
997      INTEGER :: N,IP1(NVAR),IP2(NVAR),ISING
998#ifdef FULL_ALGEBRA     
999      KPP_REAL  :: FJAC(NVAR,NVAR), E1(NVAR,NVAR)
1000      DOUBLE COMPLEX :: E2(NVAR,NVAR)
1001#else     
1002      KPP_REAL  :: FJAC(LU_NONZERO), E1(LU_NONZERO)
1003      DOUBLE COMPLEX :: E2(LU_NONZERO)
1004#endif     
1005      KPP_REAL :: Z1(N),Z2(N),Z3(N),   &
1006               F1(N),F2(N),F3(N),CONT(N),   &
1007               GAMMA,ALPHA,BETA
1008      DOUBLE COMPLEX :: BC(N)
1009      INTEGER :: i,j
1010      KPP_REAL :: S2, S3
1011!
1012      DO i=1,N
1013         S2=-F2(i)
1014         S3=-F3(i)
1015         Z1(i)=Z1(i)-F1(i)*GAMMA
1016         Z2(i)=Z2(i)+S2*ALPHA-S3*BETA
1017         Z3(i)=Z3(i)+S3*ALPHA+S2*BETA
1018      END DO
1019#ifdef FULL_ALGEBRA     
1020      CALL DGETRS ('N',N,1,E1,N,IP1,Z1,N,0) 
1021#else     
1022      CALL KppSolve (E1,Z1)
1023#endif     
1024     
1025      DO j=1,N
1026        BC(j) = DCMPLX(Z2(j),Z3(j))
1027      END DO
1028#ifdef FULL_ALGEBRA     
1029      CALL ZGETRS ('N',N,1,E2,N,IP2,BC,N,0) 
1030#else     
1031      CALL KppSolveCmplx (E2,BC)
1032#endif     
1033      DO j=1,N
1034        Z2(j) = DBLE( BC(j) )
1035        Z3(j) = AIMAG( BC(j) )
1036      END DO
1037
1038   END SUBROUTINE RAD_Solve
1039
1040
1041!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1042   SUBROUTINE RAD_ErrorEstimate(N,FJAC,H,Y0,Y,T,&
1043               E1,Z1,Z2,Z3,IP1,SCAL,ERR,        &
1044               FirstStep,REJECT,GAMMA)
1045!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1046      IMPLICIT NONE
1047     
1048      INTEGER :: N
1049#ifdef FULL_ALGEBRA     
1050      KPP_REAL :: FJAC(NVAR,NVAR),  E1(NVAR,NVAR)
1051#else     
1052      KPP_REAL :: FJAC(LU_NONZERO), E1(LU_NONZERO)
1053#endif     
1054      KPP_REAL :: SCAL(N),Z1(N),Z2(N),Z3(N),F1(N),F2(N), &
1055               Y0(N),Y(N),TMP(N),T,H,GAMMA
1056      INTEGER :: IP1(N), i
1057      LOGICAL FirstStep,REJECT
1058      KPP_REAL :: HEE1,HEE2,HEE3,ERR
1059
1060      HEE1 = rkE(1)/H
1061      HEE2 = rkE(2)/H
1062      HEE3 = rkE(3)/H
1063
1064      DO  i=1,N
1065         F2(i)=HEE1*Z1(i)+HEE2*Z2(i)+HEE3*Z3(i)
1066         TMP(i)=F2(i)+Y0(i)
1067      END DO
1068
1069#ifdef FULL_ALGEBRA     
1070      CALL DGETRS ('N',N,1,E1,N,IP1,TMP,N,0) 
1071#else     
1072      CALL KppSolve (E1, TMP)
1073#endif     
1074
1075      ERR=0.D0
1076      DO  i=1,N
1077         ERR=ERR+(TMP(i)/SCAL(i))**2
1078      END DO
1079      ERR=MAX(SQRT(ERR/N),1.D-10)
1080!
1081      IF (ERR < 1.D0) RETURN
1082firej:IF (FirstStep.OR.REJECT) THEN
1083          DO i=1,N
1084             TMP(i)=Y(i)+TMP(i)
1085          END DO
1086          CALL FUN_CHEM(T,TMP,F1)
1087          DO i=1,N
1088             TMP(i)=F1(i)+F2(i)
1089          END DO
1090
1091#ifdef FULL_ALGEBRA     
1092          CALL DGETRS ('N',N,1,E1,N,IP1,TMP,N,0) 
1093#else     
1094          CALL KppSolve (E1, TMP)
1095#endif     
1096          ERR=0.D0
1097          DO i=1,N
1098             ERR=ERR+(TMP(i)/SCAL(i))**2
1099          END DO
1100          ERR=MAX(SQRT(ERR/N),1.0d-10)
1101       END IF firej
1102 
1103   END SUBROUTINE RAD_ErrorEstimate
1104
1105!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1106!   KPP_REAL FUNCTION CONTR5(I,N,T,CONT)
1107!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1108!!     THIS FUNCTION CAN BE USED FOR CONTINUOUS OUTPUT. IT PROVIDES AN
1109!!     APPROXIMATION TO THE I-TH COMPONENT OF THE SOLUTION AT T.
1110!!     IT GIVES THE VALUE OF THE COLLOCATION POLYNOMIAL, DEFINED FOR
1111!!     THE STEP SUCCESSFULLY COMPUTED STEP (BY RADAU5).
1112!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1113!      IMPLICIT NONE
1114!      INTEGER :: I, N
1115!      KPP_REAL :: T, CONT(N,4)
1116!      KPP_REAL :: S
1117!      KPP_REAL, PARAMETER :: ONE = 1.0d0
1118!      S=(T-Tsol)/Hsol
1119!      CONTR5=CONT(i,1)+S* &
1120!        (CONT(i,2)+(S-rkC(2)+ONE)*(CONT(i,3)+(S-rkC(1)+ONE)*CONT(i,4)))
1121!   END FUNCTION CONTR5
1122 
1123!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1124  END SUBROUTINE RADAU5 ! AND ALL ITS INTERNAL PROCEDURES
1125!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1126
1127!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1128  SUBROUTINE FUN_CHEM(T, V, FCT)
1129!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1130
1131    USE KPP_ROOT_Parameters
1132    USE KPP_ROOT_Global
1133    USE KPP_ROOT_Function, ONLY: Fun
1134    USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST, Update_PHOTO
1135
1136    IMPLICIT NONE
1137
1138    KPP_REAL :: V(NVAR), FCT(NVAR)
1139    KPP_REAL :: T, Told
1140
1141    !Told = TIME
1142    !TIME = T
1143    !CALL Update_SUN()
1144    !CALL Update_RCONST()
1145    !CALL Update_PHOTO()
1146    !TIME = Told
1147   
1148    CALL Fun(V, FIX, RCONST, FCT)
1149   
1150    Nfun=Nfun+1
1151
1152  END SUBROUTINE FUN_CHEM
1153
1154!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1155  SUBROUTINE JAC_CHEM (T, V, JF)
1156!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1157
1158    USE KPP_ROOT_Parameters
1159    USE KPP_ROOT_Global
1160    USE KPP_ROOT_JacobianSP
1161    USE KPP_ROOT_Jacobian, ONLY: Jac_SP
1162    USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST, Update_PHOTO
1163
1164    IMPLICIT NONE
1165
1166    KPP_REAL :: V(NVAR), T, Told
1167#ifdef FULL_ALGEBRA   
1168    KPP_REAL :: JV(LU_NONZERO), JF(NVAR,NVAR)
1169    INTEGER :: i, j 
1170#else
1171    KPP_REAL :: JF(LU_NONZERO)
1172#endif   
1173
1174    !Told = TIME
1175    !TIME = T
1176    !CALL Update_SUN()
1177    !CALL Update_RCONST()
1178    !CALL Update_PHOTO()
1179    !TIME = Told
1180   
1181#ifdef FULL_ALGEBRA   
1182    CALL Jac_SP(V, FIX, RCONST, JV)
1183    DO j=1,NVAR
1184      DO i=1,NVAR
1185         JF(i,j) = 0.0d0
1186      END DO
1187    END DO
1188    DO i=1,LU_NONZERO
1189       JF(LU_IROW(i),LU_ICOL(i)) = JV(i)
1190    END DO
1191#else
1192    CALL Jac_SP(V, FIX, RCONST, JF) 
1193#endif   
1194   
1195    Njac=Njac+1
1196
1197  END SUBROUTINE JAC_CHEM
1198
1199!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1200
1201END MODULE KPP_ROOT_Integrator
Note: See TracBrowser for help on using the repository browser.