source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/kpp_seulex.f90 @ 4794

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

Merge of branch palm4u into trunk

File size: 36.9 KB
Line 
1!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
2!  SEULEX - Stiff extrapolation method based on linearly implicit Euler   !
3!  By default the code employs the KPP sparse linear algebra routines     !
4!  Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK)        !
5!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
6
7MODULE KPP_ROOT_Integrator
8
9  USE KPP_ROOT_Precision, ONLY: dp
10  USE KPP_ROOT_Parameters, ONLY: NVAR, LU_NONZERO
11  USE KPP_ROOT_Jacobian, ONLY: LU_DIAG
12  USE KPP_ROOT_LinearAlgebra
13
14  IMPLICIT NONE
15  PUBLIC
16  SAVE
17
18  ! variables from the former COMMON block /CONRA5/ are now here:
19  INTEGER :: NN, NN2, NN3, NN4
20  KPP_REAL :: TSOL, HSOL
21 
22  ! Statistics
23  INTEGER :: Nfun,Njac,Nstp,Nacc,Nrej,Ndec,Nsol,Nsng
24 
25  ! Method parameters
26 
27  ! mz_rs_20050717: TODO: use strings of IERR_NAMES for error messages
28  ! description of the error numbers IERR
29  CHARACTER(LEN=50), PARAMETER, DIMENSION(-4:1) :: IERR_NAMES = (/ &
30    'Matrix is repeatedly singular                     ', & ! -4
31    'Step size too small                               ', & ! -3
32    'More than Max_no_steps steps are needed           ', & ! -2
33    'Insufficient storage for work or iwork            ', & ! -1
34    '                                                  ', & !  0 (not used)
35    'Success                                           ' /) !  1
36
37CONTAINS
38
39!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
40  SUBROUTINE INTEGRATE( TIN, TOUT, &
41    ICNTRL_U, RCNTRL_U, ISTATUS_U, RSTATUS_U, IERR_U )
42!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
43
44    USE KPP_ROOT_Parameters, ONLY: nvar
45    USE KPP_ROOT_Global,     ONLY: atol,rtol,var
46
47    IMPLICIT NONE
48
49    KPP_REAL :: TIN  ! TIN - Start Time
50    KPP_REAL :: TOUT ! TOUT - End Time
51    INTEGER,  INTENT(IN),  OPTIONAL :: ICNTRL_U(20)
52    KPP_REAL, INTENT(IN),  OPTIONAL :: RCNTRL_U(20)
53    INTEGER,  INTENT(OUT), OPTIONAL :: ISTATUS_U(20)
54    KPP_REAL, INTENT(OUT), OPTIONAL :: RSTATUS_U(20)
55    INTEGER,  INTENT(OUT), OPTIONAL :: IERR_U
56
57    INTEGER :: Ncolumns, Ncolumns2, NRDENS
58    PARAMETER (Ncolumns=12,Ncolumns2=2+Ncolumns*(Ncolumns+3)/2,NRDENS=NVAR)
59
60    KPP_REAL :: RCNTRL(20), RSTATUS(20)
61    INTEGER ::       ICNTRL(20), ISTATUS(20)
62    INTEGER :: IERR
63    INTEGER, SAVE :: Ntotal = 0
64    KPP_REAL, SAVE :: H
65
66    H = 0.0_dp
67
68    ICNTRL(1:20) = 0
69    RCNTRL(1:20)  = 0._dp
70    ICNTRL(10)=0    !~~~> OUTPUT ROUTINE IS NOT USED DURING INTEGRATION
71
72    ! if optional parameters are given, and if they are >0,
73    ! they overwrite the default settings
74    IF (PRESENT(ICNTRL_U)) ICNTRL(:) = ICNTRL_U(:)
75    IF (PRESENT(RCNTRL_U)) RCNTRL(:) = RCNTRL_U(:)
76    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~-----
77   
78
79    CALL ATMSEULEX(NVAR,TIN,TOUT,VAR,H,RTOL,ATOL,      &
80                    RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR )
81    Ntotal = Ntotal + Nstp
82!!$    PRINT*,'NSTEPS=',Nstp,' (',Ntotal,')  T=',TIN
83
84
85    Nfun = Nfun + ISTATUS(1)
86    Njac = Njac + ISTATUS(2)
87    Nstp = Nstp + ISTATUS(3)
88    Nacc = Nacc + ISTATUS(4)
89    Nrej = Nrej + ISTATUS(5)
90    Ndec = Ndec + ISTATUS(6)
91    Nsol = Nsol + ISTATUS(7)
92
93    ! if optional parameters are given for output
94    ! use them to store information in them
95    IF (PRESENT(ISTATUS_U)) THEN
96      ISTATUS_U(:) = 0
97      ISTATUS_U(1) = Nfun ! function calls
98      ISTATUS_U(2) = Njac ! jacobian calls
99      ISTATUS_U(3) = Nstp ! steps
100      ISTATUS_U(4) = Nacc ! accepted steps
101      ISTATUS_U(5) = Nrej ! rejected steps (except at the beginning)
102      ISTATUS_U(6) = Ndec ! LU decompositions
103      ISTATUS_U(7) = Nsol ! forward/backward substitutions
104    ENDIF
105    IF (PRESENT(RSTATUS_U)) THEN
106      RSTATUS_U(:) = 0.
107      RSTATUS_U(1) = TOUT ! final time
108    ENDIF
109    IF (PRESENT(IERR_U)) IERR_U = IERR
110
111! mz_rs_20050716: IERR is returned to the user who then decides what to do
112! about it, i.e. either stop the run or ignore it.
113!!$    IF (IERR < 0) THEN
114!!$      PRINT *,'SEULEX: Unsuccessful exit at T=', TIN,' (IERR=',IERR,')'
115!!$      STOP
116!!$    ENDIF
117
118  END SUBROUTINE INTEGRATE
119 
120!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
121  SUBROUTINE ATMSEULEX( N,Tinitial,Tfinal,Y,H,RelTol,AbsTol,     &
122                        RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR )
123
124!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
125!     NUMERICAL SOLUTION OF A STIFF (OR DIFFERENTIAL ALGEBRAIC)
126!     SYSTEM OF FIRST 0RDER ORDINARY DIFFERENTIAL EQUATIONS  MY'=F(T,Y).
127!     THIS IS AN EXTRAPOLATION-ALGORITHM, BASED ON THE
128!     LINEARLY IMPLICIT EULER METHOD (WITH STEP SIZE CONTROL
129!     AND ORDER SELECTION).
130!
131!     AUTHORS: E. HAIRER AND G. WANNER
132!              UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES
133!              CH-1211 GENEVE 24, SWITZERLAND
134!              E-MAIL:  HAIRER@DIVSUN.UNIGE.CH,  WANNER@DIVSUN.UNIGE.CH
135!              INCLUSION OF DENSE OUTPUT BY E. HAIRER AND A. OSTERMANN
136!
137!     THIS CODE IS PART OF THE BOOK:
138!         E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
139!         EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
140!         SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS 14,
141!         SPRINGER-VERLAG (1991)
142!
143!     VERSION OF SEPTEMBER 30, 1995
144!
145!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
146!
147!     INPUT PARAMETERS
148!     ----------------
149!     N           DIMENSION OF THE SYSTEM
150!
151!     T           INITIAL T-VALUE
152!
153!     Y(N)        INITIAL VALUES FOR Y
154!
155!     Tend        FINAL T-VALUE (Tend-T MAY BE POSITIVE OR NEGATIVE)
156!
157!     H           INITIAL STEP SIZE GUESS;
158!                 FOR STIFF EQUATIONS WITH INITIAL TRANSIENT,
159!                 H=1.D0/(NORM OF F'), USUALLY 1.D-2 OR 1.D-3, IS GOOD.
160!                 THIS CHOICE IS NOT VERY IMPORTANT, THE CODE QUICKLY
161!                 ADAPTS ITS STEP SIZE (IF H=0.D0, THE CODE PUTS H=1.D-6
162!
163!     RelTol,AbsTol   RELATIVE AND ABSOLUTE ERROR TOLERANCES. THEY
164!                 CAN BE BOTH SCALARS OR ELSE BOTH VECTORS OF LENGTH N.
165!
166!     JAC         NAME (EXTERNAL) OF THE SUBROUTINE WHICH COMPUTES
167!                 THE PARTIAL DERIVATIVES OF F(T,Y) WITH RESPECT TO Y
168!
169!     SOLOUT      NAME (EXTERNAL) OF SUBROUTINE PROVIDING THE
170!                 NUMERICAL SOLUTION DURING INTEGRATION.
171!                 IF IOUT>=1, IT IS CALLED AFTER EVERY SUCCESSFUL STEP.
172!                 SUPPLY A DUMMY SUBROUTINE IF IOUT=0.
173!                 IT MUST HAVE THE FORM
174!                    SUBROUTINE SOLOUT (NR,TOLD,T,Y,RC,LRC,IC,LIC,N,
175!                                       RPAR,IPAR,IRTRN)
176!                    KPP_REAL T,Y(N),RC(LRC),IC(LIC)
177!                    ....
178!                 SOLOUT FURNISHES THE SOLUTION "Y" AT THE NR-TH
179!                    GRID-POINT "T" (THEREBY THE INITIAL VALUE IS
180!                    THE FIRST GRID-POINT).
181!                 "TOLD" IS THE PRECEEDING GRID-POINT.
182!                 "IRTRN" SERVES TO INTERRUPT THE INTEGRATION. IF IRTRN
183!                    IS SET <0, SEULEX RETURNS TO THE CALLING PROGRAM.
184!                 DO NOT CHANGE THE ENTRIES OF RC(LRC),IC(LIC)!
185!
186!          -----  CONTINUOUS OUTPUT (IF IOUT=2): -----
187!                 DURING CALLS TO "SOLOUT", A CONTINUOUS SOLUTION
188!                 FOR THE INTERVAL [TOLD,T] IS AVAILABLE THROUGH
189!                 THE KPP_REAL FUNCTION
190!                        >>>   CONTEX(I,S,RC,LRC,IC,LIC)   <<<
191!                 WHICH PROVIDES AN APPROXIMATION TO THE I-TH
192!                 COMPONENT OF THE SOLUTION AT THE POINT S. THE VALUE
193!                 S SHOULD LIE IN THE INTERVAL [TOLD,T].
194!
195!     IOUT        GIVES INFORMATION ON THE SUBROUTINE SOLOUT:
196!                    IOUT=0: SUBROUTINE IS NEVER CALLED
197!                    IOUT=1: SUBROUTINE IS USED FOR OUTPUT
198!                    IOUT=2: DENSE OUTPUT IS PERFORMED IN SOLOUT
199!
200!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
201!
202!     SOPHISTICATED SETTING OF PARAMETERS
203!     -----------------------------------
204!              SEVERAL PARAMETERS OF THE CODE ARE TUNED TO MAKE IT CNTRL
205!              WELL. THEY MAY BE DEFINED BY SETTING CNTRL(1),..,CNTRL(13)
206!              AS WELL AS ICNTRL(1),..,ICNTRL(4) DIFFERENT FROM ZERO.
207!              FOR ZERO INPUT, THE CODE CHOOSES DEFAULT VALUES:
208!
209!
210!~~~>     INPUT PARAMETERS:
211!
212!    Note: For input parameters equal to zero the default values of the
213!          corresponding variables are used.
214!~~~> 
215!    ICNTRL(1) = 1: F = F(y)   Independent of T (autonomous)
216!              = 0: F = F(t,y) Depends on T (non-autonomous)
217!
218!    ICNTRL(2) = 0: AbsTol, RelTol are NVAR-dimensional vectors
219!              = 1: AbsTol, RelTol are scalars
220!
221!    ICNTRL(3)  -> not used
222!
223!    ICNTRL(4)  -> maximum number of integration steps
224!        For ICNTRL(4)=0 the default value of 100000 is used
225!
226!    ICNTRL(11)  THE MAXIMUM NUMBER OF COLUMNS IN THE EXTRAPOLATION
227!              TABLE. THE DEFAULT VALUE (FOR ICNTRL(3)=0) IS 12.
228!              IF ICNTRL(3).NE.0 THEN ICNTRL(3) SHOULD BE  >= 3.
229!
230!    ICNTRL(12)  SWITCH FOR THE STEP SIZE SEQUENCE
231!              IF ICNTRL(4) == 1 THEN 1,2,3,4,6,8,12,16,24,32,48,...
232!              IF ICNTRL(4) == 2 THEN 2,3,4,6,8,12,16,24,32,48,64,...
233!              IF ICNTRL(4) == 3 THEN 1,2,3,4,5,6,7,8,9,10,...
234!              IF ICNTRL(4) == 4 THEN 2,3,4,5,6,7,8,9,10,11,...
235!              THE DEFAULT VALUE (FOR ICNTRL(4)=0) IS ICNTRL(4)=2.
236!
237!    ICNTRL(13)  PARAMETER "LAMBDA" OF DENSE OUTPUT; POSSIBLE VALUES
238!              ARE 0 AND 1; DEFAULT ICNTRL(5)=0.
239!
240!    ICNTRL(14)  = NRDENS = NUMBER OF COMPONENTS, FOR WHICH DENSE OUTPUT
241!              IS REQUIRED
242!
243!    ICNTRL(21),...,ICNTRL(NRDENS+20) INDICATE THE COMPONENTS, FOR WHICH
244!              DENSE OUTPUT IS REQUIRED
245!
246!~~~>  Real parameters
247!
248!    RCNTRL(1)  -> Hmin, lower bound for the integration step size
249!                  It is strongly recommended to keep Hmin = ZERO
250!    RCNTRL(2)  -> Hmax, upper bound for the integration step size
251!    RCNTRL(3)  -> Hstart, starting value for the integration step size
252!    RCNTRL(4)  -> FacMin, lower bound on step decrease factor (default=0.2)
253!    RCNTRL(5)  -> FacMax, upper bound on step increase factor (default=6)
254!    RCNTRL(6)  -> FacRej, step decrease factor after multiple rejections
255!                 (default=0.1)
256!    RCNTRL(7)  -> FacSafe, by which the new step is slightly smaller
257!                  than the predicted value  (default=0.9)
258!    RCNTRL(8)  -> ThetaMin. If Newton convergence rate smaller
259!                  than ThetaMin the Jacobian is not recomputed;
260!                  (default=0.001). Increase cntrl(3), to 0.01 say, when
261!                   Jacobian evaluations are costly. for small systems it
262!                   should be smaller.
263!    RCNTRL(9)  -> not used
264!    RCNTRL(10,11) -> FAC1,FAC2 (parameters for step size selection)
265!    RCNTRL(12,13) -> FAC3,FAC4 (parameters for order selection)
266!    RCNTRL(14,15) -> FacSafe1, FacSafe2
267!                     Safety factors for step size prediction
268!                     HNEW=H*FacSafe2*(FacSafe1*TOL/ERR)**(1/(J-1))
269!    RCNTRL(16:19) -> WorkFcn, WorkJac, WorkDec, WorkSol
270!                     estimated computational work
271!
272!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
273!
274!     OUTPUT PARAMETERS
275!     -----------------
276!     T           T-VALUE WHERE THE SOLUTION IS COMPUTED
277!                 (AFTER SUCCESSFUL RETURN T=Tend)
278!
279!     Y(N)        SOLUTION AT T
280!
281!     H           PREDICTED STEP SIZE OF THE LAST ACCEPTED STEP
282!
283!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
284!          DECLARATIONS
285!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
286      IMPLICIT NONE
287      INTEGER :: N, IERR, ITOL, Max_no_steps, Ncolumns, Nsequence, Lambda, &
288                 NRDENS, i, Ncolumns2, NRD, IOUT
289      KPP_REAL :: Y(NVAR),AbsTol(*),RelTol(*)
290      KPP_REAL :: Tinitial, Tfinal, Roundoff, Hmin, Hmax,         &
291                 FacMin, FacMax, FAC1, FAC2, FAC3, FAC4, FacSafe1,     &
292                 FacSafe2, H, Hstart,WorkFcn,WorkJac, WorkDec, WorkSol,&
293                 WorkRow, FacRej, FacSafe, ThetaMin, T
294      LOGICAL :: AUTNMS
295      KPP_REAL :: RCNTRL(20), RSTATUS(20)
296      INTEGER ::       ICNTRL(20), ISTATUS(20)
297      KPP_REAL, PARAMETER :: ZERO = 0.0d0
298
299!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
300!        SETTING THE PARAMETERS
301!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
302   Nfun=0
303   Njac=0
304   Nstp=0
305   Nacc=0
306   Nrej=0
307   Ndec=0
308   Nsol=0
309       
310   IERR = 0
311
312   IF (ICNTRL(1) == 0) THEN
313      AUTNMS = .FALSE. 
314   ELSE
315      AUTNMS = .TRUE.
316   END IF
317
318!~~~>  For Scalar tolerances (ICNTRL(1)/=0) the code uses AbsTol(1) and RelTol(1)
319!~~~>  For Vector tolerances (ICNTRL(1)==0) the code uses AbsTol(1:NVAR) and RelTol(1:NVAR)
320   IF (ICNTRL(2) == 0) THEN
321      ITOL = 1
322   ELSE
323      ITOL = 0
324   END IF
325
326!~~~> Max_no_steps: the maximum number of time steps admitted
327   IF (ICNTRL(4) == 0) THEN
328      Max_no_steps = 100000
329   ELSEIF (ICNTRL(4) > 0) THEN
330      Max_no_steps=ICNTRL(4)
331   ELSE
332      PRINT * ,'User-selected ICNTRL(4)=',ICNTRL(4)
333      CALL SEULEX_ErrorMsg(-1,Tinitial,ZERO,IERR);
334   END IF
335
336!~~~> IOUT = use (or not) the output routine
337   IOUT = ICNTRL(10)
338   IF ( IOUT<0 .OR. IOUT>2 ) THEN
339       PRINT * ,'User-selected ICNTRL(10)=',ICNTRL(10)
340       IOUT = 0
341   END IF
342   
343!~~~> Ncolumns:  maximum number of columns in the extrapolation
344   IF (ICNTRL(11)==0) THEN
345       Ncolumns=12
346   ELSEIF (ICNTRL(11) > 2) THEN
347       Ncolumns=ICNTRL(11)
348   ELSE   
349       PRINT * ,'User-selected ICNTRL(11)=',ICNTRL(11)
350       CALL SEULEX_ErrorMsg(-2,Tinitial,ZERO,IERR);
351   END IF
352
353!~~~> Nsequence: choice of step size sequence
354   IF (ICNTRL(12)==0) THEN
355      Nsequence = 2
356   ELSEIF ( (ICNTRL(12)>0).AND.(ICNTRL(12)<5) ) THEN
357      Nsequence = ICNTRL(4)
358   ELSE
359      PRINT * ,'User-selected ICNTRL(12)=',ICNTRL(12)
360      CALL SEULEX_ErrorMsg(-3,Tinitial,ZERO,IERR)
361   END IF
362
363!~~~> LAMBDA: parameter for dense output
364   LAMBDA = ICNTRL(13)
365   IF ( LAMBDA < 0 .OR. LAMBDA >= 2 ) THEN
366       PRINT * ,'User-selected ICNTRL(13)=',ICNTRL(13)
367       CALL SEULEX_ErrorMsg(-4,Tinitial,ZERO,IERR)
368   END IF
369     
370!~~~>- NRDENS:  number of dense output components
371   NRDENS=ICNTRL(14)
372   IF ( (NRDENS < 0) .OR. (NRDENS > N)  ) THEN
373       PRINT * ,'User-selected ICNTRL(14)=',ICNTRL(14)
374       CALL SEULEX_ErrorMsg(-5,Tinitial,ZERO,IERR)
375   END IF
376
377
378!~~~>  Unit roundoff (1+Roundoff>1)
379      Roundoff = WLAMCH('E')
380
381!~~~>  Lower bound on the step size: (positive value)
382   IF (RCNTRL(1) == ZERO) THEN
383      Hmin = ZERO
384   ELSEIF (RCNTRL(1) > ZERO) THEN
385      Hmin = RCNTRL(1)
386   ELSE
387      PRINT * , 'User-selected RCNTRL(1)=', RCNTRL(1)
388      CALL SEULEX_ErrorMsg(-7,Tinitial,ZERO,IERR)
389      RETURN
390   END IF
391   
392!~~~>  Upper bound on the step size: (positive value)
393   IF (RCNTRL(2) == ZERO) THEN
394      Hmax = ABS(Tfinal-Tinitial)
395   ELSEIF (RCNTRL(2) > ZERO) THEN
396      Hmax = MIN(ABS(RCNTRL(2)),ABS(Tfinal-Tinitial))
397   ELSE
398      PRINT * , 'User-selected RCNTRL(2)=', RCNTRL(2)
399      CALL SEULEX_ErrorMsg(-7,Tinitial,ZERO,IERR)
400      RETURN
401   END IF
402   
403!~~~>  Starting step size: (positive value)
404   IF (RCNTRL(3) == ZERO) THEN
405      Hstart = MAX(Hmin,Roundoff)
406   ELSEIF (RCNTRL(3) > ZERO) THEN
407      Hstart = MIN(ABS(RCNTRL(3)),ABS(Tfinal-Tinitial))
408   ELSE
409      PRINT * , 'User-selected Hstart: RCNTRL(3)=', RCNTRL(3)
410      CALL SEULEX_ErrorMsg(-7,Tinitial,ZERO,IERR)
411      RETURN
412   END IF
413   
414   
415!~~~>  Step size can be changed s.t.  FacMin < Hnew/Hexit < FacMax
416   IF (RCNTRL(4) == ZERO) THEN
417      FacMin = 0.2_dp
418   ELSEIF (RCNTRL(4) > ZERO) THEN
419      FacMin = RCNTRL(4)
420   ELSE
421      PRINT * , 'User-selected FacMin: RCNTRL(4)=', RCNTRL(4)
422      CALL SEULEX_ErrorMsg(-8,Tinitial,ZERO,IERR)
423      RETURN
424   END IF
425   IF (RCNTRL(5) == ZERO) THEN
426      FacMax = 10.0_dp
427   ELSEIF (RCNTRL(5) > ZERO) THEN
428      FacMax = RCNTRL(5)
429   ELSE
430      PRINT * , 'User-selected FacMax: RCNTRL(5)=', RCNTRL(5)
431      CALL SEULEX_ErrorMsg(-8,Tinitial,ZERO,IERR)
432      RETURN
433   END IF
434!~~~>   FacRej: Factor to decrease step after 2 succesive rejections
435   IF (RCNTRL(6) == ZERO) THEN
436      FacRej = 0.1_dp
437   ELSEIF (RCNTRL(6) > ZERO) THEN
438      FacRej = RCNTRL(6)
439   ELSE
440      PRINT * , 'User-selected FacRej: RCNTRL(6)=', RCNTRL(6)
441      CALL SEULEX_ErrorMsg(-8,Tinitial,ZERO,IERR)
442      RETURN
443   END IF
444!~~~>   FacSafe: Safety Factor in the computation of new step size
445   IF (RCNTRL(7) == ZERO) THEN
446      FacSafe = 0.9_dp
447   ELSEIF (RCNTRL(7) > ZERO) THEN
448      FacSafe = RCNTRL(7)
449   ELSE
450      PRINT * , 'User-selected FacSafe: RCNTRL(7)=', RCNTRL(7)
451      CALL SEULEX_ErrorMsg(-8,Tinitial,ZERO,IERR)
452      RETURN
453   END IF
454
455!~~~> ThetaMin: DECIDES WHETHER THE JACOBIAN SHOULD BE RECOMPUTED;
456!              INCREASE WORK(3), TO 0.01 SAY, WHEN JACOBIAN EVALUATIONS
457!              ARE COSTLY. FOR SMALL SYSTEMS WORK(3) SHOULD BE SMALLER.
458   IF(RCNTRL(8) == 0.D0)THEN
459       ThetaMin = 1.0d-3
460   ELSE
461       ThetaMin = RCNTRL(8)
462   END IF
463
464!~~~>  FAC1,FAC2:     PARAMETERS FOR STEP SIZE SELECTION
465!              THE NEW STEP SIZE FOR THE J-TH DIAGONAL ENTRY IS
466!              CHOSEN SUBJECT TO THE RESTRICTION
467!                 FACMIN/WORK(5) <= HNEW(J)/HOLD <= 1/FACMIN
468!              WHERE FACMIN=WORK(4)**(1/(J-1))
469      IF(RCNTRL(10) == 0.D0)THEN
470         FAC1=0.1D0
471      ELSE
472         FAC1=RCNTRL(10)
473      END IF
474      IF(RCNTRL(11) == 0.D0)THEN
475         FAC2=4.0D0
476      ELSE
477         FAC2=RCNTRL(11)
478      END IF
479!~~~>  FAC3, FAC4:   PARAMETERS FOR THE ORDER SELECTION
480!              ORDER IS DECREASED IF    W(K-1) <= W(K)*WORK(6)
481!              ORDER IS INCREASED IF    W(K) <= W(K-1)*WORK(7)
482      IF(RCNTRL(12) == 0.D0)THEN
483         FAC3=0.7D0
484      ELSE
485         FAC3=RCNTRL(12)
486      END IF
487      IF(RCNTRL(13) == 0.D0)THEN
488         FAC4=0.9D0
489      ELSE
490         FAC4=RCNTRL(13)
491      END IF
492!~~~>- FacSafe1, FacSafe2: safety factors for step size prediction
493!             HNEW=H*WORK(9)*(WORK(8)*TOL/ERR)**(1/(J-1))
494      IF(RCNTRL(14) == 0.D0)THEN
495         FacSafe1=0.6D0
496      ELSE
497         FacSafe1=RCNTRL(14)
498      END IF
499      IF(RCNTRL(15) == 0.D0)THEN
500         FacSafe2=0.93D0
501      ELSE
502         FacSafe2=RCNTRL(15)
503      END IF
504     
505!~~~> WorkFcn: estimated computational work for a calls to FCN
506      IF(RCNTRL(16) == 0.D0)THEN
507         WorkFcn=1.D0
508      ELSE
509         WorkFcn=RCNTRL(16)
510      END IF
511!~~~> WorkJac: estimated computational work for calls to JAC
512      IF(RCNTRL(17) == 0.D0)THEN
513         WorkJac=5.D0
514      ELSE
515         WorkJac=RCNTRL(17)
516      END IF
517!~~~> WorkDec: estimated computational work for calls to DEC
518      IF(RCNTRL(18) == 0.D0)THEN
519         WorkDec=1.D0
520      ELSE
521         WorkDec=RCNTRL(18)
522      END IF
523!~~~> WorkSol: estimated computational work for calls to SOL
524      IF(RCNTRL(19) == 0.D0)THEN
525         WorkSol=1.D0
526      ELSE
527         WorkSol=RCNTRL(19)
528      END IF
529      WorkRow=WorkFcn+WorkSol
530
531!~~~>  Check if tolerances are reasonable
532      IF (ITOL == 0) THEN
533         IF (AbsTol(1) <= 0.D0.OR.RelTol(1) <= 10.D0*Roundoff) THEN
534            PRINT * , ' Scalar AbsTol = ',AbsTol(1)
535            PRINT * , ' Scalar RelTol = ',RelTol(1)
536            CALL SEULEX_ErrorMsg(-9,Tinitial,ZERO,IERR)
537         END IF
538      ELSE
539         DO i=1,N
540            IF (AbsTol(i) <= 0.D0.OR.RelTol(i) <= 10.D0*Roundoff) THEN
541              PRINT * , ' AbsTol(',i,') = ',AbsTol(i)
542              PRINT * , ' RelTol(',i,') = ',RelTol(i)
543              CALL SEULEX_ErrorMsg(-9,Tinitial,ZERO,IERR)
544            END IF
545         END DO
546      END IF
547   
548    IF (IERR < 0) RETURN
549       
550!~~~>---- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK -----
551      Ncolumns2=(Ncolumns*(Ncolumns+1))/2
552      NRD=MAX(1,NRDENS)
553
554      T = Tinitial
555!~~~> CALL TO CORE INTEGRATOR
556      CALL SEULEX_Integrator(N,T,Tfinal,Y,Hmax,H,Ncolumns,RelTol,AbsTol,ITOL,  &
557                IOUT,IERR,Max_no_steps,Roundoff,Nsequence,AUTNMS,  &
558                FAC1,FAC2,FAC3,FAC4,ThetaMin,FacSafe1,FacSafe2,WorkJac,  &
559                WorkDec,WorkRow,Ncolumns2,NRD,LAMBDA,Nstp)
560       
561      ISTATUS(1)=Nfun
562      ISTATUS(2)=Njac
563      ISTATUS(3)=Nstp
564      ISTATUS(4)=Nacc
565      ISTATUS(5)=Nrej
566      ISTATUS(6)=Ndec
567      ISTATUS(7)=Nsol
568
569      END SUBROUTINE ATMSEULEX
570     
571!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
572 SUBROUTINE SEULEX_ErrorMsg(Code,T,H,IERR)
573!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
574!    Handles all error messages
575!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
576
577   KPP_REAL, INTENT(IN) :: T, H
578   INTEGER, INTENT(IN)  :: Code
579   INTEGER, INTENT(OUT) :: IERR
580
581   IERR = Code
582   PRINT * , &
583     'Forced exit from SEULEX due to the following error:'
584
585   SELECT CASE (Code)
586    CASE (-1)
587      PRINT * , '--> Improper value for maximal no of steps'
588    CASE (-2)
589      PRINT * , '--> Improper value for maximum no of columns in extrapolation'
590    CASE (-3)
591      PRINT * , '--> Improper value for step size sequence'
592    CASE (-4)
593      PRINT * , '--> Improper value for Lambda (must be 0/1)'
594    CASE (-5)
595      PRINT * , '--> Improper number of dense output components'
596    CASE (-6)
597      PRINT * , '--> Improper parameters for second order equations'
598    CASE (-7)
599      PRINT * , '--> Hmin/Hmax/Hstart must be positive'
600    CASE (-8)
601      PRINT * , '--> FacMin/FacMax/FacRej must be positive'
602    CASE (-9)
603      PRINT * , '--> Improper tolerance values'
604    CASE (-10)
605      PRINT * , '--> No of steps exceeds maximum bound'
606    CASE (-11)
607      PRINT * , '--> Step size too small: T + 10*H = T', &
608            ' or H < Roundoff'
609    CASE (-12)
610      PRINT * , '--> Matrix is repeatedly singular'
611    CASE DEFAULT
612      PRINT *, 'Unknown Error code: ', Code
613   END SELECT
614
615   PRINT *, "T=", T, "and H=", H
616
617 END SUBROUTINE SEULEX_ErrorMsg
618     
619
620!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
621      SUBROUTINE SEULEX_Integrator(N,T,Tend,Y,Hmax,H,Ncolumns,RelTol,AbsTol,ITOL,&
622       IOUT,IERR,Max_no_steps,Roundoff,Nsequence,AUTNMS,             &
623       FAC1,FAC2,FAC3,FAC4,ThetaMin,FacSafe1,FacSafe2,WorkJac,       &
624       WorkDec,WorkRow,Ncolumns2,NRD,LAMBDA,Nstp)
625!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
626!     CORE INTEGRATOR FOR SEULEX
627!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
628!         DECLARATIONS
629!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
630       USE KPP_ROOT_Parameters
631       USE KPP_ROOT_Jacobian
632       IMPLICIT KPP_REAL (A-H,O-Z)
633       IMPLICIT INTEGER (I-N)
634
635       INTEGER ::  N, Ncolumns, Ncolumns2, K, KC, KRIGHT, KLR, KK, KRN,&
636       KOPT, NRD
637       KPP_REAL ::  Y(NVAR),DY(NVAR),FX(NVAR),YHH(NVAR)
638       KPP_REAL ::  DYH(NVAR), DEL(NVAR), WH(NVAR)
639       KPP_REAL ::  SCAL(NVAR), HH(Ncolumns), W(Ncolumns), A(Ncolumns)
640#ifdef FULL_ALGEBRA   
641       KPP_REAL :: FJAC(NVAR,NVAR)
642#else       
643       KPP_REAL :: FJAC(LU_NONZERO)
644#endif       
645       KPP_REAL Table(Ncolumns,N)
646       INTEGER IP(N),NJ(Ncolumns),IPHES(N),ICOMP(NRD)
647       KPP_REAL RelTol(*),AbsTol(*)
648       KPP_REAL FSAFE(Ncolumns2,NRD),FACUL(Ncolumns),E(N,N),DENS((Ncolumns+2)*NRD)
649       LOGICAL REJECT,LAST,ATOV,CALJAC,CALHES,AUTNMS
650
651       KPP_REAL TOLDD,HHH,NNRD
652       COMMON /COSEU/TOLDD,HHH,NNRD,KRIGHT
653
654!~~~> COMPUTE COEFFICIENTS FOR DENSE OUTPUT
655       IF (IOUT == 2) THEN
656          NNRD=NRD
657!~~~> COMPUTE THE FACTORIALS --------
658          FACUL(1)=1.D0
659          DO i=1,Ncolumns-1
660             FACUL(i+1)=i*FACUL(i)
661          END DO
662       END IF
663
664!~~~> DEFINE THE STEP SIZE SEQUENCE
665       IF (Nsequence == 1) THEN
666           NJ(1)=1
667           NJ(2)=2
668           NJ(3)=3
669           DO I=4,Ncolumns
670              NJ(i)=2*NJ(I-2)
671           END DO
672       END IF
673       IF (Nsequence == 2) THEN
674           NJ(1)=2
675           NJ(2)=3
676           DO I=3,Ncolumns
677              NJ(i)=2*NJ(I-2)
678           END DO
679       END IF
680       DO i=1,Ncolumns
681          IF (Nsequence == 3) NJ(i)=I
682          IF (Nsequence == 4) NJ(i)=I+1
683       END DO
684       A(1)=WorkJac+NJ(1)*WorkRow+WorkDec
685       DO I=2,Ncolumns
686          A(i)=A(i-1)+(NJ(i)-1)*WorkRow+WorkDec
687       END DO
688       K=MAX0(3,MIN0(Ncolumns-2,INT(-DLOG10(RelTol(1)+AbsTol(1))*.6D0+1.5D0)))
689       
690       ! T = Tinitial
691       HmaxN = MIN(ABS(Hmax),ABS(Tend-T))
692       IF (ABS(H) <= 10.D0*Roundoff) H=1.0D-6
693       H=MIN(ABS(H),HmaxN)
694       Theta=2*ABS(ThetaMin)
695       ERR=0.D0
696       W(1)=1.D30
697       DO i=1,N
698          IF (ITOL == 0) THEN
699            SCAL(i)=AbsTol(1)+RelTol(1)*DABS(Y(i))
700          ELSE
701            SCAL(i)=AbsTol(i)+RelTol(i)*DABS(Y(i))
702          END IF
703       END DO
704       CALJAC=.FALSE.
705       REJECT=.FALSE.
706       LAST=.FALSE.
707   10  CONTINUE
708       IF (REJECT) Theta=2*ABS(ThetaMin)
709       ATOV=.FALSE.
710!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
711!~~~> IS Tend REACHED IN THE NEXT STEP?
712!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
713       H1=Tend-T
714       IF (H1 <= Roundoff) GO TO 110
715       HOPT=H
716       H=MIN(H,H1,HmaxN)
717       IF (H >= H1-Roundoff) LAST=.TRUE.
718       IF (AUTNMS) THEN
719          CALL FUN_CHEM(T,Y,DY)
720       END IF
721       IF (Theta > ThetaMin.AND..NOT.CALJAC) THEN
722!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
723!  COMPUTATION OF THE JACOBIAN
724!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
725          CALL JAC_CHEM(T,Y,FJAC)
726          CALJAC=.TRUE.
727          CALHES=.FALSE.
728       END IF
729!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
730!~~~> THE FIRST AND LAST STEP
731!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
732      IF (Nstp == 0.OR.LAST) THEN
733        IPT=0
734        Nstp=Nstp+1
735        DO J=1,K
736         KC=J
737         CALL SEUL(J,N,T,Y,DY,FX,FJAC,LFJAC,E,LE,IP,H,Ncolumns, &
738                  HmaxN,Table,SCAL,NJ,HH,W,A,YHH,DYH,DEL,WH,ERR,FacSafe1,FAC, &
739                  FAC1,FAC2,FacSafe2,Theta,Nfun,Ndec,Nsol, &
740                  ERROLD,IPHES,ICOMP,AUTNMS,REJECT, &
741                  ATOV,FSAFE,Ncolumns2,NRD,IOUT,IPT,CALHES)
742         IF (ATOV) GOTO 10
743         IF (J > 1 .AND. ERR <= 1.d0) GOTO 60
744        END DO
745        GO TO 55
746      END IF
747!~~~> BASIC INTEGRATION STEP
748   30 CONTINUE
749      IPT=0
750      Nstp=Nstp+1
751      IF (Nstp >= Max_no_steps) GOTO 120
752      KC=K-1
753      DO J=1,KC
754       CALL SEUL(J,N,T,Y,DY,FX,FJAC,LFJAC,E,LE,IP,H,Ncolumns,&
755                HmaxN,Table,SCAL,NJ,HH,W,A,YHH,DYH,DEL,WH,ERR,FacSafe1,&
756                FAC,FAC1,FAC2,FacSafe2,Theta,Nfun,Ndec,Nsol,&
757                ERROLD,IPHES,ICOMP,AUTNMS,REJECT,&
758                ATOV,FSAFE,Ncolumns2,NRD,IOUT,IPT,CALHES)
759       IF (ATOV) GOTO 10
760      END DO
761!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
762!~~~> CONVERGENCE MONITOR
763!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
764      IF (K == 2.OR.REJECT) GO TO 50
765      IF (ERR <= 1.D0) GO TO 60
766      IF (ERR > DBLE(NJ(K+1)*NJ(K))*4.D0) GO TO 100
767   50 CALL SEUL(K,N,T,Y,DY,FX,FJAC,LFJAC,E,LE,IP,H,Ncolumns,&
768                HmaxN,Table,SCAL,NJ,HH,W,A,YHH,DYH,DEL,WH,ERR,FacSafe1,&
769                FAC,FAC1,FAC2,FacSafe2,Theta,Nfun,Ndec,Nsol,&
770                ERROLD,IPHES,ICOMP,AUTNMS,REJECT,&
771                ATOV,FSAFE,Ncolumns2,NRD,IOUT,IPT,CALHES)
772      IF (ATOV) GOTO 10
773      KC=K
774      IF (ERR <= 1.D0) GO TO 60
775!~~~> HOPE FOR CONVERGENCE IN LINE K+1
776   55 IF (ERR > DBLE(NJ(K+1))*2.D0) GO TO 100
777      KC=K+1
778      CALL SEUL(KC,N,T,Y,DY,FX,FJAC,LFJAC,E,LE,IP,H,Ncolumns,&
779                HmaxN,Table,SCAL,NJ,HH,W,A,YHH,DYH,DEL,WH,ERR,FacSafe1,&
780                FAC,FAC1,FAC2,FacSafe2,Theta,Nfun,Ndec,Nsol,&
781                ERROLD,IPHES,ICOMP,AUTNMS,REJECT,&
782                ATOV,FSAFE,Ncolumns2,NRD,IOUT,IPT,CALHES)
783      IF (ATOV) GOTO 10
784      IF (ERR > 1.D0) GO TO 100
785      !Adi IF ((ERR > 1.D0).and.(H.gt.Hmin)) GO TO 100
786!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
787!~~~> STEP IS ACCEPTED
788!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
789   60 TOLD=T
790      T=T+H
791      IF (IOUT == 2) THEN
792         KRIGHT=KC
793         DO i=1,NRD
794            DENS(i)=Y(ICOMP(i))
795         END DO
796      END IF
797      DO i=1,N
798         T1I=Table(1,I)
799         IF (ITOL == 0) THEN
800            SCAL(i)=AbsTol(1)+RelTol(1)*DABS(T1I)
801         ELSE
802            SCAL(i)=AbsTol(i)+RelTol(i)*DABS(T1I)
803         END IF
804         Y(i)=T1I
805      END DO
806      Nacc=Nacc+1
807      CALJAC=.FALSE.
808      IF (IOUT == 2) THEN
809         TOLDD=TOLD
810         HHH=H
811         DO i=1,NRD
812            DENS(NRD+I)=Y(ICOMP(i))
813         END DO
814         DO KLR=1,KRIGHT-1
815!~~~> COMPUTE DIFFERENCES
816            IF (KLR >= 2) THEN
817               DO KK=KLR,KC
818                  LBEG=((KK+1)*KK)/2
819                  LEND=LBEG-KK+2
820                  DO L=LBEG,LEND,-1
821                     DO i=1,NRD
822                        FSAFE(L,I)=FSAFE(L,I)-FSAFE(L-1,I)
823                     END DO
824                  END DO
825               END DO
826             END IF
827!~~~> COMPUTE DERIVATIVES AT RIGHT END ----
828             DO KK=KLR+LAMBDA,KC
829                FACNJ=NJ(KK)
830                FACNJ=FACNJ**KLR/FACUL(KLR+1)
831                IPT=((KK+1)*KK)/2
832                DO  I=1,NRD
833                   KRN=(KK-LAMBDA+1)*NRD
834                   DENS(KRN+I)=FSAFE(IPT,I)*FACNJ
835                END DO
836             END DO
837             DO J=KLR+LAMBDA+1,KC
838               DBLENJ=NJ(J)
839               DO L=J,KLR+LAMBDA+1,-1
840                 FACTOR=DBLENJ/NJ(L-1)-1.D0
841                 DO i=1,NRD
842                   KRN=(L-LAMBDA+1)*NRD+I
843                   DENS(KRN-NRD)=DENS(KRN)+(DENS(KRN)-DENS(KRN-NRD))/FACTOR
844                 END DO
845               END DO
846             END DO
847         END DO
848!~~~>  COMPUTE THE COEFFICIENTS OF THE INTERPOLATION POLYNOMIAL
849         DO IN=1,NRD
850            DO J=1,KRIGHT
851               II=NRD*J+IN
852               DENS(II)=DENS(II)-DENS(II-NRD)
853            END DO
854         END DO
855      END IF
856!~~~> COMPUTE OPTIMAL ORDER
857      IF (KC == 2) THEN
858         KOPT=3
859         IF (REJECT) KOPT=2
860         GO TO 80
861      END IF
862      IF (KC <= K) THEN
863         KOPT=KC
864         IF (W(KC-1) < W(KC)*FAC3) KOPT=KC-1
865         IF (W(KC) < W(KC-1)*FAC4) KOPT=MIN0(KC+1,Ncolumns-1)
866      ELSE
867         KOPT=KC-1
868         IF (KC > 3.AND.W(KC-2) < W(KC-1)*FAC3) KOPT=KC-2
869         IF (W(KC) < W(KOPT)*FAC4) KOPT=MIN0(KC,Ncolumns-1)
870      END IF
871!~~~> AFTER A REJECTED STEP
872   80 IF (REJECT) THEN
873         K=MIN0(KOPT,KC)
874         H=MIN(H,HH(K))
875         REJECT=.FALSE.
876         GO TO 10
877      END IF
878!~~~> COMPUTE STEP SIZE FOR NEXT STEP
879      IF (KOPT <= KC) THEN
880         H=HH(KOPT)
881      ELSE
882         IF (KC < K.AND.W(KC) < W(KC-1)*FAC4) THEN
883            H=HH(KC)*A(KOPT+1)/A(KC)
884         ELSE
885            H=HH(KC)*A(KOPT)/A(KC)
886         END IF
887      END IF
888      K=KOPT
889      !Adi H = MAX(H, Hmin)
890      GO TO 10
891!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
892!~~~> STEP IS REJECTED
893!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
894  100 K=MIN0(K,KC)
895      IF (K > 2.AND.W(K-1) < W(K)*FAC3) K=K-1
896      Nrej=Nrej+1
897      H=HH(K)
898      LAST=.FALSE.
899      REJECT=.TRUE.
900      IF (CALJAC) GOTO 30
901      GO TO 10
902!~~~> SOLUTION EXIT
903  110 CONTINUE
904      H=HOPT
905      IERR=1
906      RETURN
907!~~~> FAIL EXIT
908  120 WRITE (6,979) T,H
909  979 FORMAT(' EXIT OF SEULEX AT T=',D14.7,'   H=',D14.7)
910      IERR=-1
911      RETURN
912                                                     
913                             
914      END SUBROUTINE SEULEX_Integrator
915
916
917!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
918      SUBROUTINE SEUL(JJ,N,T,Y,DY,FX,FJAC,LFJAC,E,LE,IP,&
919               H,Ncolumns,HmaxN,Table,SCAL,NJ,HH,W,A,YH,DYH,DEL,WH,ERR,FacSafe1, &
920               FAC,FAC1,FAC2,FacSafe2,Theta,Nfun,Ndec,Nsol,&
921               ERROLD,IPHES,ICOMP,                             &
922               AUTNMS,REJECT,ATOV,FSAFE,Ncolumns2,NRD,IOUT,    &
923               IPT,CALHES)
924!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
925!~~~> THIS SUBROUTINE COMPUTES THE J-TH LINE OF THE
926!~~~> EXTRAPOLATION TABLE AND PROVIDES AN ESTIMATE
927!~~~> OF THE OPTIMAL STEP SIZE
928      USE KPP_ROOT_Parameters
929      USE KPP_ROOT_Jacobian
930      IMPLICIT KPP_REAL (A-H,O-Z)
931      IMPLICIT INTEGER (I-N)
932      INTEGER :: Ncolumns, Ncolumns2, N, NRD
933      KPP_REAL :: Y(NVAR),YH(NVAR),DY(NVAR),FX(NVAR),DYH(NVAR)
934      KPP_REAL :: DEL(NVAR),WH(NVAR),SCAL(NVAR),HH(Ncolumns),W(Ncolumns),A(Ncolumns)
935#ifdef FULL_ALGEBRA   
936      KPP_REAL :: FJAC(NVAR,NVAR), E(NVAR,NVAR)
937#else
938      KPP_REAL :: FJAC(LU_NONZERO), E(LU_NONZERO)
939#endif     
940      KPP_REAL :: Table(Ncolumns,NVAR)
941      KPP_REAL :: FSAFE(Ncolumns2,NRD)
942      INTEGER :: IP(N),NJ(Ncolumns),IPHES(N),ICOMP(NRD)
943      LOGICAL ATOV,REJECT,AUTNMS,CALHES
944     
945!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
946!  COMPUTE THE MATRIX E AND ITS DECOMPOSITION
947!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
948      HJ=H/NJ(JJ)
949      HJI=1.D0/HJ
950#ifdef FULL_ALGEBRA   
951      DO j=1,N
952        DO  i=1,N
953           E(i,j)=-FJAC(i,j)
954        END DO
955        E(j,j)=E(j,j)+HJI
956      END DO
957      CALL DGETRF(N,N,E,N,IP,ISING)
958#else
959      DO  i=1,LU_NONZERO
960         E(i)=-FJAC(i)
961      END DO
962      DO j=1,N
963         E(LU_DIAG(j))=E(LU_DIAG(j))+HJI
964      END DO
965      CALL KppDecomp (E,ISING)
966#endif     
967      Ndec=Ndec+1
968      IF (ISING.NE.0) GOTO 79
969!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
970!~~~> STARTING PROCEDURE
971!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
972       IF (.NOT.AUTNMS) THEN
973          CALL FUN_CHEM(T+HJ,Y,DY)
974       END IF
975       DO i=1,N
976          YH(i)=Y(i)
977          DEL(i)=DY(i)
978       END DO
979#ifdef FULL_ALGEBRA     
980       CALL DGETRS ('N',N,1,E,N,IP,DEL,N,ISING)
981#else
982       CALL KppSolve (E,DEL)
983#endif       
984       Nsol=Nsol+1
985       M=NJ(JJ)
986       IF (IOUT == 2.AND.M == JJ) THEN
987          IPT=IPT+1
988          DO i=1,NRD
989             FSAFE(IPT,I)=DEL(ICOMP(i))
990          END DO
991       END IF
992!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
993!~~~> SEMI-IMPLICIT EULER METHOD
994!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
995       IF (M > 1) THEN
996          DO MM=1,M-1
997             DO i=1,N
998                YH(i)=YH(i)+DEL(i)
999             END DO
1000             IF (AUTNMS) THEN
1001                CALL FUN_CHEM(T+HJ*MM,YH,DYH)
1002             ELSE
1003                CALL FUN_CHEM(T+HJ*(MM+1),YH,DYH)
1004             END IF
1005             
1006             IF (MM == 1.AND.JJ <= 2) THEN
1007!~~~> STABILITY CHECK
1008                DEL1=0.D0
1009                DO i=1,N
1010                   DEL1=DEL1+(DEL(i)/SCAL(i))**2
1011                END DO
1012                DEL1=SQRT(DEL1)
1013                IF (.NOT.AUTNMS) THEN
1014                   CALL FUN_CHEM(T+HJ,YH,WH)
1015                   
1016                   DO i=1,N
1017                      DEL(i)=WH(i)-DEL(i)*HJI
1018                   END DO
1019                ELSE
1020                   DO i=1,N
1021                      DEL(i)=DYH(i)-DEL(i)*HJI
1022                   END DO
1023                END IF
1024#ifdef FULL_ALGEBRA     
1025                CALL DGETRS ('N',N,1,E,N,IP,DEL,N,ISING)
1026#else
1027                CALL KppSolve (E,DEL)
1028#endif
1029                Nsol=Nsol+1
1030                DEL2=0.D0
1031                DO i=1,N
1032                   DEL2=DEL2+(DEL(i)/SCAL(i))**2
1033                END DO
1034                DEL2=SQRT(DEL2)
1035                Theta=DEL2/MAX(1.D0,DEL1)
1036                IF (Theta > 1.D0) GOTO 79
1037             END IF
1038#ifdef FULL_ALGEBRA     
1039             CALL DGETRS ('N',N,1,E,N,IP,DYH,N,ISING)
1040#else
1041             CALL KppSolve (E,DYH)
1042#endif             
1043             Nsol=Nsol+1
1044             DO i=1,N
1045                DEL(i)=DYH(i)
1046             END DO
1047             IF (IOUT == 2.AND.MM >= M-JJ) THEN
1048                IPT=IPT+1
1049                DO i=1,NRD
1050                   FSAFE(IPT,i)=DEL(ICOMP(i))
1051                END DO
1052             END IF
1053          END DO
1054       END IF
1055       DO i=1,N
1056          Table(JJ,I)=YH(i)+DEL(i)
1057       END DO
1058!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1059!~~~> POLYNOMIAL EXTRAPOLATION
1060!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1061       IF (JJ == 1) RETURN
1062       DO L=JJ,2,-1
1063          FAC=(DBLE(NJ(JJ))/DBLE(NJ(L-1)))-1.D0
1064          DO i=1,N
1065             Table(L-1,I)=Table(L,I)+(Table(L,I)-Table(L-1,I))/FAC
1066          END DO
1067       END DO
1068       ERR=0.D0
1069       DO i=1,N
1070          ERR=ERR+MIN(ABS((Table(1,I)-Table(2,I)))/SCAL(i),1.D15)**2
1071       END DO
1072       IF (ERR >= 1.D30) GOTO 79
1073       ERR=SQRT(ERR/DBLE(N))
1074       IF (JJ > 2.AND.ERR >= ERROLD) GOTO 79
1075       ERROLD=MAX(4*ERR,1.D0)
1076!~~~> COMPUTE OPTIMAL STEP SIZES
1077       EXPO=1.D0/JJ
1078       FACMIN=FAC1**EXPO
1079       FAC=MIN(FAC2/FACMIN,MAX(FACMIN,(ERR/FacSafe1)**EXPO/FacSafe2))
1080       FAC=1.D0/FAC
1081       HH(JJ)=MIN(H*FAC,HmaxN)
1082       W(JJ)=A(JJ)/HH(JJ)
1083       RETURN
1084   79  ATOV=.TRUE.
1085       H=H*0.5D0
1086       REJECT=.TRUE.
1087       RETURN
1088      END SUBROUTINE SEUL
1089
1090
1091
1092!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1093  SUBROUTINE FUN_CHEM( T, V, FCT )
1094!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1095
1096    USE KPP_ROOT_Parameters
1097    USE KPP_ROOT_Global
1098    USE KPP_ROOT_Function, ONLY: Fun
1099    USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST, Update_PHOTO
1100
1101    IMPLICIT NONE
1102
1103    KPP_REAL :: V(NVAR), FCT(NVAR)
1104    KPP_REAL :: T, TOLD
1105
1106    !TOLD = TIME
1107    !TIME = T
1108    !CALL Update_SUN()
1109    !CALL Update_RCONST()
1110    !CALL Update_PHOTO()
1111    !TIME = TOLD
1112    CALL Fun(V, FIX, RCONST, FCT)
1113    Nfun=Nfun+1
1114
1115  END SUBROUTINE FUN_CHEM
1116
1117
1118!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1119  SUBROUTINE JAC_CHEM ( T, V, Jcb )
1120!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1121
1122    USE KPP_ROOT_Parameters
1123    USE KPP_ROOT_Global
1124    USE KPP_ROOT_Jacobian, ONLY: Jac_SP
1125    USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST, Update_PHOTO
1126
1127    IMPLICIT NONE
1128
1129    KPP_REAL :: V(NVAR), T, TOLD
1130#ifdef FULL_ALGEBRA   
1131    KPP_REAL :: JV(LU_NONZERO), Jcb(NVAR,NVAR)
1132    INTEGER :: i,j
1133#else
1134    KPP_REAL :: Jcb(LU_NONZERO)
1135#endif   
1136
1137    !TOLD = TIME
1138    !TIME = T
1139    !CALL Update_SUN()
1140    !CALL Update_RCONST()
1141    !CALL Update_PHOTO()
1142    !TIME = TOLD
1143   
1144#ifdef FULL_ALGEBRA   
1145    CALL Jac_SP(V, FIX, RCONST, JV)
1146    DO j=1,NVAR
1147      DO i=1,NVAR
1148         Jcb(i,j) = 0.0D0
1149      END DO
1150    END DO
1151    DO i=1,LU_NONZERO
1152       Jcb(LU_IROW(i),LU_ICOL(i)) = JV(i)
1153    END DO
1154#else
1155    CALL Jac_SP(V, FIX, RCONST, Jcb) 
1156#endif   
1157    Njac=Njac+1
1158   
1159  END SUBROUTINE JAC_CHEM
1160
1161!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1162
1163END MODULE KPP_ROOT_Integrator
Note: See TracBrowser for help on using the repository browser.