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 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
7 MODULE KPP_ROOT_Integrator
9 USE KPP_ROOT_Precision
, ONLY
: dp
10 USE KPP_ROOT_Jacobian
, ONLY
: NVAR
, LU_NONZERO
, LU_DIAG
11 !! USE KPP_ROOT_LinearAlgebra
17 ! variables from the former COMMON block /CONRA5/ are now here:
18 INTEGER :: NN
, NN2
, NN3
, NN4
19 KPP_REAL
:: TSOL
, HSOL
22 INTEGER :: Nfun
,Njac
,Nstp
,Nacc
,Nrej
,Ndec
,Nsol
,Nsng
26 ! mz_rs_20050717: TODO: use strings of IERR_NAMES for error messages
27 ! description of the error numbers IERR
28 CHARACTER(LEN
=50), PARAMETER, DIMENSION(-4:1) :: IERR_NAMES
= (/ &
29 'Matrix is repeatedly singular ', & ! -4
30 'Step size too small ', & ! -3
31 'More than Max_no_steps steps are needed ', & ! -2
32 'Insufficient storage for work or iwork ', & ! -1
38 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
39 SUBROUTINE KPP_ROOT_INTEGRATE( TIN
, TOUT
, &
40 FIX
, VAR
, RCONST
, ATOL
, RTOL
, &
41 ICNTRL_U
, RCNTRL_U
, ISTATUS_U
, RSTATUS_U
, IERR_U
)
43 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
45 USE KPP_ROOT_Parameters
, ONLY
: nvar
46 !! USE KPP_ROOT_Global, ONLY: atol,rtol,var
50 KPP_REAL
, INTENT(INOUT
), DIMENSION(NFIX
) :: FIX
51 KPP_REAL
, INTENT(INOUT
), DIMENSION(NVAR
) :: VAR
52 KPP_REAL
, INTENT(IN
), DIMENSION(NSPEC
) :: ATOL
, RTOL
53 KPP_REAL
, INTENT(IN
), DIMENSION(NREACT
) :: RCONST
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
63 INTEGER :: Ncolumns
, Ncolumns2
, NRDENS
64 PARAMETER (Ncolumns
=12,Ncolumns2
=2+Ncolumns
*(Ncolumns
+3)/2,NRDENS
=NVAR
)
66 KPP_REAL
:: RCNTRL(20), RSTATUS(20)
67 INTEGER :: ICNTRL(20), ISTATUS(20)
69 INTEGER, SAVE :: Ntotal
= 0
76 ICNTRL(10)=0 !~~~> OUTPUT ROUTINE IS NOT USED DURING INTEGRATION
78 ! if optional parameters are given, and if they are >0,
79 ! they overwrite the default settings
80 IF (PRESENT(ICNTRL_U
)) ICNTRL(:) = ICNTRL_U(:)
81 IF (PRESENT(RCNTRL_U
)) RCNTRL(:) = RCNTRL_U(:)
82 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~-----
85 CALL ATMSEULEX(NVAR
,TIN
,TOUT
,VAR
,H
,RTOL
,ATOL
, &
86 RCNTRL
,ICNTRL
,RSTATUS
,ISTATUS
,IERR
)
87 Ntotal
= Ntotal
+ Nstp
88 !!$ PRINT*,'NSTEPS=',Nstp,' (',Ntotal,') T=',TIN
91 Nfun
= Nfun
+ ISTATUS(1)
92 Njac
= Njac
+ ISTATUS(2)
93 Nstp
= Nstp
+ ISTATUS(3)
94 Nacc
= Nacc
+ ISTATUS(4)
95 Nrej
= Nrej
+ ISTATUS(5)
96 Ndec
= Ndec
+ ISTATUS(6)
97 Nsol
= Nsol
+ ISTATUS(7)
99 ! if optional parameters are given for output
100 ! use them to store information in them
101 IF (PRESENT(ISTATUS_U
)) THEN
103 ISTATUS_U(1) = Nfun
! function calls
104 ISTATUS_U(2) = Njac
! jacobian calls
105 ISTATUS_U(3) = Nstp
! steps
106 ISTATUS_U(4) = Nacc
! accepted steps
107 ISTATUS_U(5) = Nrej
! rejected steps (except at the beginning)
108 ISTATUS_U(6) = Ndec
! LU decompositions
109 ISTATUS_U(7) = Nsol
! forward/backward substitutions
111 IF (PRESENT(RSTATUS_U
)) THEN
113 RSTATUS_U(1) = TOUT
! final time
115 IF (PRESENT(IERR_U
)) IERR_U
= IERR
117 ! mz_rs_20050716: IERR is returned to the user who then decides what to do
118 ! about it, i.e. either stop the run or ignore it.
119 !!$ IF (IERR < 0) THEN
120 !!$ PRINT *,'SEULEX: Unsuccessful exit at T=', TIN,' (IERR=',IERR,')'
124 END SUBROUTINE KPP_ROOT_ INTEGRATE
126 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
127 SUBROUTINE KPP_ROOT_ATMSEULEX( N
,Tinitial
,Tfinal
,Y
,H
,RelTol
,AbsTol
, &
128 RCNTRL
,ICNTRL
,RSTATUS
,ISTATUS
,IERR
)
130 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
131 ! NUMERICAL SOLUTION OF A STIFF (OR DIFFERENTIAL ALGEBRAIC)
132 ! SYSTEM OF FIRST 0RDER ORDINARY DIFFERENTIAL EQUATIONS MY'=F(T,Y).
133 ! THIS IS AN EXTRAPOLATION-ALGORITHM, BASED ON THE
134 ! LINEARLY IMPLICIT EULER METHOD (WITH STEP SIZE CONTROL
135 ! AND ORDER SELECTION).
137 ! AUTHORS: E. HAIRER AND G. WANNER
138 ! UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES
139 ! CH-1211 GENEVE 24, SWITZERLAND
140 ! E-MAIL: HAIRER@DIVSUN.UNIGE.CH, WANNER@DIVSUN.UNIGE.CH
141 ! INCLUSION OF DENSE OUTPUT BY E. HAIRER AND A. OSTERMANN
143 ! THIS CODE IS PART OF THE BOOK:
144 ! E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
145 ! EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
146 ! SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS 14,
147 ! SPRINGER-VERLAG (1991)
149 ! VERSION OF SEPTEMBER 30, 1995
151 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
155 ! N DIMENSION OF THE SYSTEM
159 ! Y(N) INITIAL VALUES FOR Y
161 ! Tend FINAL T-VALUE (Tend-T MAY BE POSITIVE OR NEGATIVE)
163 ! H INITIAL STEP SIZE GUESS;
164 ! FOR STIFF EQUATIONS WITH INITIAL TRANSIENT,
165 ! H=1.D0/(NORM OF F'), USUALLY 1.D-2 OR 1.D-3, IS GOOD.
166 ! THIS CHOICE IS NOT VERY IMPORTANT, THE CODE QUICKLY
167 ! ADAPTS ITS STEP SIZE (IF H=0.D0, THE CODE PUTS H=1.D-6
169 ! RelTol,AbsTol RELATIVE AND ABSOLUTE ERROR TOLERANCES. THEY
170 ! CAN BE BOTH SCALARS OR ELSE BOTH VECTORS OF LENGTH N.
172 ! JAC NAME (EXTERNAL) OF THE SUBROUTINE WHICH COMPUTES
173 ! THE PARTIAL DERIVATIVES OF F(T,Y) WITH RESPECT TO Y
175 ! SOLOUT NAME (EXTERNAL) OF SUBROUTINE PROVIDING THE
176 ! NUMERICAL SOLUTION DURING INTEGRATION.
177 ! IF IOUT>=1, IT IS CALLED AFTER EVERY SUCCESSFUL STEP.
178 ! SUPPLY A DUMMY SUBROUTINE IF IOUT=0.
179 ! IT MUST HAVE THE FORM
180 ! SUBROUTINE SOLOUT (NR,TOLD,T,Y,RC,LRC,IC,LIC,N,
182 ! KPP_REAL T,Y(N),RC(LRC),IC(LIC)
184 ! SOLOUT FURNISHES THE SOLUTION "Y" AT THE NR-TH
185 ! GRID-POINT "T" (THEREBY THE INITIAL VALUE IS
186 ! THE FIRST GRID-POINT).
187 ! "TOLD" IS THE PRECEEDING GRID-POINT.
188 ! "IRTRN" SERVES TO INTERRUPT THE INTEGRATION. IF IRTRN
189 ! IS SET <0, SEULEX RETURNS TO THE CALLING PROGRAM.
190 ! DO NOT CHANGE THE ENTRIES OF RC(LRC),IC(LIC)!
192 ! ----- CONTINUOUS OUTPUT (IF IOUT=2): -----
193 ! DURING CALLS TO "SOLOUT", A CONTINUOUS SOLUTION
194 ! FOR THE INTERVAL [TOLD,T] IS AVAILABLE THROUGH
195 ! THE KPP_REAL FUNCTION
196 ! >>> CONTEX(I,S,RC,LRC,IC,LIC) <<<
197 ! WHICH PROVIDES AN APPROXIMATION TO THE I-TH
198 ! COMPONENT OF THE SOLUTION AT THE POINT S. THE VALUE
199 ! S SHOULD LIE IN THE INTERVAL [TOLD,T].
201 ! IOUT GIVES INFORMATION ON THE SUBROUTINE SOLOUT:
202 ! IOUT=0: SUBROUTINE IS NEVER CALLED
203 ! IOUT=1: SUBROUTINE IS USED FOR OUTPUT
204 ! IOUT=2: DENSE OUTPUT IS PERFORMED IN SOLOUT
206 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
208 ! SOPHISTICATED SETTING OF PARAMETERS
209 ! -----------------------------------
210 ! SEVERAL PARAMETERS OF THE CODE ARE TUNED TO MAKE IT CNTRL
211 ! WELL. THEY MAY BE DEFINED BY SETTING CNTRL(1),..,CNTRL(13)
212 ! AS WELL AS ICNTRL(1),..,ICNTRL(4) DIFFERENT FROM ZERO.
213 ! FOR ZERO INPUT, THE CODE CHOOSES DEFAULT VALUES:
216 !~~~> INPUT PARAMETERS:
218 ! Note: For input parameters equal to zero the default values of the
219 ! corresponding variables are used.
221 ! ICNTRL(1) = 1: F = F(y) Independent of T (autonomous)
222 ! = 0: F = F(t,y) Depends on T (non-autonomous)
224 ! ICNTRL(2) = 0: AbsTol, RelTol are NVAR-dimensional vectors
225 ! = 1: AbsTol, RelTol are scalars
227 ! ICNTRL(3) -> not used
229 ! ICNTRL(4) -> maximum number of integration steps
230 ! For ICNTRL(4)=0 the default value of 100000 is used
232 ! ICNTRL(11) THE MAXIMUM NUMBER OF COLUMNS IN THE EXTRAPOLATION
233 ! TABLE. THE DEFAULT VALUE (FOR ICNTRL(3)=0) IS 12.
234 ! IF ICNTRL(3).NE.0 THEN ICNTRL(3) SHOULD BE >= 3.
236 ! ICNTRL(12) SWITCH FOR THE STEP SIZE SEQUENCE
237 ! IF ICNTRL(4) == 1 THEN 1,2,3,4,6,8,12,16,24,32,48,...
238 ! IF ICNTRL(4) == 2 THEN 2,3,4,6,8,12,16,24,32,48,64,...
239 ! IF ICNTRL(4) == 3 THEN 1,2,3,4,5,6,7,8,9,10,...
240 ! IF ICNTRL(4) == 4 THEN 2,3,4,5,6,7,8,9,10,11,...
241 ! THE DEFAULT VALUE (FOR ICNTRL(4)=0) IS ICNTRL(4)=2.
243 ! ICNTRL(13) PARAMETER "LAMBDA" OF DENSE OUTPUT; POSSIBLE VALUES
244 ! ARE 0 AND 1; DEFAULT ICNTRL(5)=0.
246 ! ICNTRL(14) = NRDENS = NUMBER OF COMPONENTS, FOR WHICH DENSE OUTPUT
249 ! ICNTRL(21),...,ICNTRL(NRDENS+20) INDICATE THE COMPONENTS, FOR WHICH
250 ! DENSE OUTPUT IS REQUIRED
252 !~~~> Real parameters
254 ! RCNTRL(1) -> Hmin, lower bound for the integration step size
255 ! It is strongly recommended to keep Hmin = ZERO
256 ! RCNTRL(2) -> Hmax, upper bound for the integration step size
257 ! RCNTRL(3) -> Hstart, starting value for the integration step size
258 ! RCNTRL(4) -> FacMin, lower bound on step decrease factor (default=0.2)
259 ! RCNTRL(5) -> FacMax, upper bound on step increase factor (default=6)
260 ! RCNTRL(6) -> FacRej, step decrease factor after multiple rejections
262 ! RCNTRL(7) -> FacSafe, by which the new step is slightly smaller
263 ! than the predicted value (default=0.9)
264 ! RCNTRL(8) -> ThetaMin. If Newton convergence rate smaller
265 ! than ThetaMin the Jacobian is not recomputed;
266 ! (default=0.001). Increase cntrl(3), to 0.01 say, when
267 ! Jacobian evaluations are costly. for small systems it
269 ! RCNTRL(9) -> not used
270 ! RCNTRL(10,11) -> FAC1,FAC2 (parameters for step size selection)
271 ! RCNTRL(12,13) -> FAC3,FAC4 (parameters for order selection)
272 ! RCNTRL(14,15) -> FacSafe1, FacSafe2
273 ! Safety factors for step size prediction
274 ! HNEW=H*FacSafe2*(FacSafe1*TOL/ERR)**(1/(J-1))
275 ! RCNTRL(16:19) -> WorkFcn, WorkJac, WorkDec, WorkSol
276 ! estimated computational work
278 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
282 ! T T-VALUE WHERE THE SOLUTION IS COMPUTED
283 ! (AFTER SUCCESSFUL RETURN T=Tend)
287 ! H PREDICTED STEP SIZE OF THE LAST ACCEPTED STEP
289 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
291 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
293 INTEGER :: N
, IERR
, ITOL
, Max_no_steps
, Ncolumns
, Nsequence
, Lambda
, &
294 NRDENS
, i
, Ncolumns2
, NRD
, IOUT
295 KPP_REAL
:: Y(NVAR
),AbsTol(*),RelTol(*)
296 KPP_REAL
:: Tinitial
, Tfinal
, Roundoff
, Hmin
, Hmax
, &
297 FacMin
, FacMax
, FAC1
, FAC2
, FAC3
, FAC4
, FacSafe1
, &
298 FacSafe2
, H
, Hstart
,WorkFcn
,WorkJac
, WorkDec
, WorkSol
,&
299 WorkRow
, FacRej
, FacSafe
, ThetaMin
, T
301 KPP_REAL
:: RCNTRL(20), RSTATUS(20)
302 INTEGER :: ICNTRL(20), ISTATUS(20)
303 KPP_REAL
, PARAMETER :: ZERO
= 0.0d0
305 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
306 ! SETTING THE PARAMETERS
307 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
318 IF (ICNTRL(1) == 0) THEN
324 !~~~> For Scalar tolerances (ICNTRL(1)/=0) the code uses AbsTol(1) and RelTol(1)
325 !~~~> For Vector tolerances (ICNTRL(1)==0) the code uses AbsTol(1:NVAR) and RelTol(1:NVAR)
326 IF (ICNTRL(2) == 0) THEN
332 !~~~> Max_no_steps: the maximum number of time steps admitted
333 IF (ICNTRL(4) == 0) THEN
334 Max_no_steps
= 100000
335 ELSEIF (ICNTRL(4) > 0) THEN
336 Max_no_steps
=ICNTRL(4)
338 PRINT * ,'User-selected ICNTRL(4)=',ICNTRL(4)
339 CALL SEULEX_ErrorMsg(-1,Tinitial
,ZERO
,IERR
);
342 !~~~> IOUT = use (or not) the output routine
344 IF ( IOUT
<0 .OR
. IOUT
>2 ) THEN
345 PRINT * ,'User-selected ICNTRL(10)=',ICNTRL(10)
349 !~~~> Ncolumns: maximum number of columns in the extrapolation
350 IF (ICNTRL(11)==0) THEN
352 ELSEIF (ICNTRL(11) > 2) THEN
355 PRINT * ,'User-selected ICNTRL(11)=',ICNTRL(11)
356 CALL SEULEX_ErrorMsg(-2,Tinitial
,ZERO
,IERR
);
359 !~~~> Nsequence: choice of step size sequence
360 IF (ICNTRL(12)==0) THEN
362 ELSEIF ( (ICNTRL(12)>0).AND
.(ICNTRL(12)<5) ) THEN
363 Nsequence
= ICNTRL(4)
365 PRINT * ,'User-selected ICNTRL(12)=',ICNTRL(12)
366 CALL SEULEX_ErrorMsg(-3,Tinitial
,ZERO
,IERR
)
369 !~~~> LAMBDA: parameter for dense output
371 IF ( LAMBDA
< 0 .OR
. LAMBDA
>= 2 ) THEN
372 PRINT * ,'User-selected ICNTRL(13)=',ICNTRL(13)
373 CALL SEULEX_ErrorMsg(-4,Tinitial
,ZERO
,IERR
)
376 !~~~>- NRDENS: number of dense output components
378 IF ( (NRDENS
< 0) .OR
. (NRDENS
> N
) ) THEN
379 PRINT * ,'User-selected ICNTRL(14)=',ICNTRL(14)
380 CALL SEULEX_ErrorMsg(-5,Tinitial
,ZERO
,IERR
)
384 !~~~> Unit roundoff (1+Roundoff>1)
385 Roundoff
= WLAMCH('E')
387 !~~~> Lower bound on the step size: (positive value)
388 IF (RCNTRL(1) == ZERO
) THEN
390 ELSEIF (RCNTRL(1) > ZERO
) THEN
393 PRINT * , 'User-selected RCNTRL(1)=', RCNTRL(1)
394 CALL SEULEX_ErrorMsg(-7,Tinitial
,ZERO
,IERR
)
398 !~~~> Upper bound on the step size: (positive value)
399 IF (RCNTRL(2) == ZERO
) THEN
400 Hmax
= ABS(Tfinal
-Tinitial
)
401 ELSEIF (RCNTRL(2) > ZERO
) THEN
402 Hmax
= MIN(ABS(RCNTRL(2)),ABS(Tfinal
-Tinitial
))
404 PRINT * , 'User-selected RCNTRL(2)=', RCNTRL(2)
405 CALL SEULEX_ErrorMsg(-7,Tinitial
,ZERO
,IERR
)
409 !~~~> Starting step size: (positive value)
410 IF (RCNTRL(3) == ZERO
) THEN
411 Hstart
= MAX(Hmin
,Roundoff
)
412 ELSEIF (RCNTRL(3) > ZERO
) THEN
413 Hstart
= MIN(ABS(RCNTRL(3)),ABS(Tfinal
-Tinitial
))
415 PRINT * , 'User-selected Hstart: RCNTRL(3)=', RCNTRL(3)
416 CALL SEULEX_ErrorMsg(-7,Tinitial
,ZERO
,IERR
)
421 !~~~> Step size can be changed s.t. FacMin < Hnew/Hexit < FacMax
422 IF (RCNTRL(4) == ZERO
) THEN
424 ELSEIF (RCNTRL(4) > ZERO
) THEN
427 PRINT * , 'User-selected FacMin: RCNTRL(4)=', RCNTRL(4)
428 CALL SEULEX_ErrorMsg(-8,Tinitial
,ZERO
,IERR
)
431 IF (RCNTRL(5) == ZERO
) THEN
433 ELSEIF (RCNTRL(5) > ZERO
) THEN
436 PRINT * , 'User-selected FacMax: RCNTRL(5)=', RCNTRL(5)
437 CALL SEULEX_ErrorMsg(-8,Tinitial
,ZERO
,IERR
)
440 !~~~> FacRej: Factor to decrease step after 2 succesive rejections
441 IF (RCNTRL(6) == ZERO
) THEN
443 ELSEIF (RCNTRL(6) > ZERO
) THEN
446 PRINT * , 'User-selected FacRej: RCNTRL(6)=', RCNTRL(6)
447 CALL SEULEX_ErrorMsg(-8,Tinitial
,ZERO
,IERR
)
450 !~~~> FacSafe: Safety Factor in the computation of new step size
451 IF (RCNTRL(7) == ZERO
) THEN
453 ELSEIF (RCNTRL(7) > ZERO
) THEN
456 PRINT * , 'User-selected FacSafe: RCNTRL(7)=', RCNTRL(7)
457 CALL SEULEX_ErrorMsg(-8,Tinitial
,ZERO
,IERR
)
461 !~~~> ThetaMin: DECIDES WHETHER THE JACOBIAN SHOULD BE RECOMPUTED;
462 ! INCREASE WORK(3), TO 0.01 SAY, WHEN JACOBIAN EVALUATIONS
463 ! ARE COSTLY. FOR SMALL SYSTEMS WORK(3) SHOULD BE SMALLER.
464 IF(RCNTRL(8) == 0.D0
)THEN
470 !~~~> FAC1,FAC2: PARAMETERS FOR STEP SIZE SELECTION
471 ! THE NEW STEP SIZE FOR THE J-TH DIAGONAL ENTRY IS
472 ! CHOSEN SUBJECT TO THE RESTRICTION
473 ! FACMIN/WORK(5) <= HNEW(J)/HOLD <= 1/FACMIN
474 ! WHERE FACMIN=WORK(4)**(1/(J-1))
475 IF(RCNTRL(10) == 0.D0
)THEN
480 IF(RCNTRL(11) == 0.D0
)THEN
485 !~~~> FAC3, FAC4: PARAMETERS FOR THE ORDER SELECTION
486 ! ORDER IS DECREASED IF W(K-1) <= W(K)*WORK(6)
487 ! ORDER IS INCREASED IF W(K) <= W(K-1)*WORK(7)
488 IF(RCNTRL(12) == 0.D0
)THEN
493 IF(RCNTRL(13) == 0.D0
)THEN
498 !~~~>- FacSafe1, FacSafe2: safety factors for step size prediction
499 ! HNEW=H*WORK(9)*(WORK(8)*TOL/ERR)**(1/(J-1))
500 IF(RCNTRL(14) == 0.D0
)THEN
505 IF(RCNTRL(15) == 0.D0
)THEN
511 !~~~> WorkFcn: estimated computational work for a calls to FCN
512 IF(RCNTRL(16) == 0.D0
)THEN
517 !~~~> WorkJac: estimated computational work for calls to JAC
518 IF(RCNTRL(17) == 0.D0
)THEN
523 !~~~> WorkDec: estimated computational work for calls to DEC
524 IF(RCNTRL(18) == 0.D0
)THEN
529 !~~~> WorkSol: estimated computational work for calls to SOL
530 IF(RCNTRL(19) == 0.D0
)THEN
535 WorkRow
=WorkFcn
+WorkSol
537 !~~~> Check if tolerances are reasonable
539 IF (AbsTol(1) <= 0.D0
.OR
.RelTol(1) <= 10.D0
*Roundoff
) THEN
540 PRINT * , ' Scalar AbsTol = ',AbsTol(1)
541 PRINT * , ' Scalar RelTol = ',RelTol(1)
542 CALL SEULEX_ErrorMsg(-9,Tinitial
,ZERO
,IERR
)
546 IF (AbsTol(i
) <= 0.D0
.OR
.RelTol(i
) <= 10.D0
*Roundoff
) THEN
547 PRINT * , ' AbsTol(',i
,') = ',AbsTol(i
)
548 PRINT * , ' RelTol(',i
,') = ',RelTol(i
)
549 CALL SEULEX_ErrorMsg(-9,Tinitial
,ZERO
,IERR
)
556 !~~~>---- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK -----
557 Ncolumns2
=(Ncolumns
*(Ncolumns
+1))/2
561 !~~~> CALL TO CORE INTEGRATOR
562 CALL SEULEX_Integrator(N
,T
,Tfinal
,Y
,Hmax
,H
,Ncolumns
,RelTol
,AbsTol
,ITOL
, &
563 IOUT
,IERR
,Max_no_steps
,Roundoff
,Nsequence
,AUTNMS
, &
564 FAC1
,FAC2
,FAC3
,FAC4
,ThetaMin
,FacSafe1
,FacSafe2
,WorkJac
, &
565 WorkDec
,WorkRow
,Ncolumns2
,NRD
,LAMBDA
,Nstp
)
575 END SUBROUTINE KPP_ROOT_ATMSEULEX
577 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
578 SUBROUTINE KPP_ROOT_SEULEX_ErrorMsg(Code
,T
,H
,IERR
)
579 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
580 ! Handles all error messages
581 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
583 KPP_REAL
, INTENT(IN
) :: T
, H
584 INTEGER, INTENT(IN
) :: Code
585 INTEGER, INTENT(OUT
) :: IERR
589 'Forced exit from SEULEX due to the following error:'
593 PRINT * , '--> Improper value for maximal no of steps'
595 PRINT * , '--> Improper value for maximum no of columns in extrapolation'
597 PRINT * , '--> Improper value for step size sequence'
599 PRINT * , '--> Improper value for Lambda (must be 0/1)'
601 PRINT * , '--> Improper number of dense output components'
603 PRINT * , '--> Improper parameters for second order equations'
605 PRINT * , '--> Hmin/Hmax/Hstart must be positive'
607 PRINT * , '--> FacMin/FacMax/FacRej must be positive'
609 PRINT * , '--> Improper tolerance values'
611 PRINT * , '--> No of steps exceeds maximum bound'
613 PRINT * , '--> Step size too small: T + 10*H = T', &
616 PRINT * , '--> Matrix is repeatedly singular'
618 PRINT *, 'Unknown Error code: ', Code
621 PRINT *, "T=", T
, "and H=", H
623 END SUBROUTINE KPP_ROOT_SEULEX_ErrorMsg
626 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
627 SUBROUTINE KPP_ROOT_SEULEX_Integrator(N
,T
,Tend
,Y
,Hmax
,H
,Ncolumns
,RelTol
,AbsTol
,ITOL
,&
628 IOUT
,IERR
,Max_no_steps
,Roundoff
,Nsequence
,AUTNMS
, &
629 FAC1
,FAC2
,FAC3
,FAC4
,ThetaMin
,FacSafe1
,FacSafe2
,WorkJac
, &
630 WorkDec
,WorkRow
,Ncolumns2
,NRD
,LAMBDA
,Nstp
)
631 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
632 ! CORE INTEGRATOR FOR SEULEX
633 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
635 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
636 USE KPP_ROOT_Parameters
637 USE KPP_ROOT_Jacobian
638 IMPLICIT KPP_REAL (A
-H
,O
-Z
)
639 IMPLICIT INTEGER (I
-N
)
641 INTEGER :: N
, Ncolumns
, Ncolumns2
, K
, KC
, KRIGHT
, KLR
, KK
, KRN
,&
643 KPP_REAL
:: Y(NVAR
),DY(NVAR
),FX(NVAR
),YHH(NVAR
)
644 KPP_REAL
:: DYH(NVAR
), DEL(NVAR
), WH(NVAR
)
645 KPP_REAL
:: SCAL(NVAR
), HH(Ncolumns
), W(Ncolumns
), A(Ncolumns
)
647 KPP_REAL
:: FJAC(NVAR
,NVAR
)
649 KPP_REAL
:: FJAC(LU_NONZERO
)
651 KPP_REAL
Table(Ncolumns
,N
)
652 INTEGER IP(N
),NJ(Ncolumns
),IPHES(N
),ICOMP(NRD
)
653 KPP_REAL
RelTol(*),AbsTol(*)
654 KPP_REAL
FSAFE(Ncolumns2
,NRD
),FACUL(Ncolumns
),E(N
,N
),DENS((Ncolumns
+2)*NRD
)
655 LOGICAL REJECT
,LAST
,ATOV
,CALJAC
,CALHES
,AUTNMS
657 KPP_REAL TOLDD
,HHH
,NNRD
658 COMMON /COSEU
/TOLDD
,HHH
,NNRD
,KRIGHT
660 !~~~> COMPUTE COEFFICIENTS FOR DENSE OUTPUT
663 !~~~> COMPUTE THE FACTORIALS --------
666 FACUL(i
+1)=i
*FACUL(i
)
670 !~~~> DEFINE THE STEP SIZE SEQUENCE
671 IF (Nsequence
== 1) THEN
679 IF (Nsequence
== 2) THEN
687 IF (Nsequence
== 3) NJ(i
)=I
688 IF (Nsequence
== 4) NJ(i
)=I
+1
690 A(1)=WorkJac
+NJ(1)*WorkRow
+WorkDec
692 A(i
)=A(i
-1)+(NJ(i
)-1)*WorkRow
+WorkDec
694 K
=MAX0(3,MIN0(Ncolumns
-2,INT(-DLOG10(RelTol(1)+AbsTol(1))*.6D0
+1.5D0
)))
697 HmaxN
= MIN(ABS(Hmax
),ABS(Tend
-T
))
698 IF (ABS(H
) <= 10.D0
*Roundoff
) H
=1.0D
-6
700 Theta
=2*ABS(ThetaMin
)
705 SCAL(i
)=AbsTol(1)+RelTol(1)*DABS(Y(i
))
707 SCAL(i
)=AbsTol(i
)+RelTol(i
)*DABS(Y(i
))
714 IF (REJECT
) Theta
=2*ABS(ThetaMin
)
716 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
717 !~~~> IS Tend REACHED IN THE NEXT STEP?
718 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
720 IF (H1
<= Roundoff
) GO
TO 110
723 IF (H
>= H1
-Roundoff
) LAST
=.TRUE
.
725 CALL FUN_CHEM(T
,Y
,DY
)
727 IF (Theta
> ThetaMin
.AND
..NOT
.CALJAC
) THEN
728 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
729 ! COMPUTATION OF THE JACOBIAN
730 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
731 CALL JAC_CHEM(T
,Y
,FJAC
)
735 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
736 !~~~> THE FIRST AND LAST STEP
737 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
738 IF (Nstp
== 0.OR
.LAST
) THEN
743 CALL SEUL(J
,N
,T
,Y
,DY
,FX
,FJAC
,LFJAC
,E
,LE
,IP
,H
,Ncolumns
, &
744 HmaxN
,Table
,SCAL
,NJ
,HH
,W
,A
,YHH
,DYH
,DEL
,WH
,ERR
,FacSafe1
,FAC
, &
745 FAC1
,FAC2
,FacSafe2
,Theta
,Nfun
,Ndec
,Nsol
, &
746 ERROLD
,IPHES
,ICOMP
,AUTNMS
,REJECT
, &
747 ATOV
,FSAFE
,Ncolumns2
,NRD
,IOUT
,IPT
,CALHES
)
749 IF (J
> 1 .AND
. ERR
<= 1.d0
) GOTO 60
753 !~~~> BASIC INTEGRATION STEP
757 IF (Nstp
>= Max_no_steps
) GOTO 120
760 CALL SEUL(J
,N
,T
,Y
,DY
,FX
,FJAC
,LFJAC
,E
,LE
,IP
,H
,Ncolumns
,&
761 HmaxN
,Table
,SCAL
,NJ
,HH
,W
,A
,YHH
,DYH
,DEL
,WH
,ERR
,FacSafe1
,&
762 FAC
,FAC1
,FAC2
,FacSafe2
,Theta
,Nfun
,Ndec
,Nsol
,&
763 ERROLD
,IPHES
,ICOMP
,AUTNMS
,REJECT
,&
764 ATOV
,FSAFE
,Ncolumns2
,NRD
,IOUT
,IPT
,CALHES
)
767 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
768 !~~~> CONVERGENCE MONITOR
769 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
770 IF (K
== 2.OR
.REJECT
) GO
TO 50
771 IF (ERR
<= 1.D0
) GO
TO 60
772 IF (ERR
> DBLE(NJ(K
+1)*NJ(K
))*4.D0
) GO
TO 100
773 50 CALL SEUL(K
,N
,T
,Y
,DY
,FX
,FJAC
,LFJAC
,E
,LE
,IP
,H
,Ncolumns
,&
774 HmaxN
,Table
,SCAL
,NJ
,HH
,W
,A
,YHH
,DYH
,DEL
,WH
,ERR
,FacSafe1
,&
775 FAC
,FAC1
,FAC2
,FacSafe2
,Theta
,Nfun
,Ndec
,Nsol
,&
776 ERROLD
,IPHES
,ICOMP
,AUTNMS
,REJECT
,&
777 ATOV
,FSAFE
,Ncolumns2
,NRD
,IOUT
,IPT
,CALHES
)
780 IF (ERR
<= 1.D0
) GO
TO 60
781 !~~~> HOPE FOR CONVERGENCE IN LINE K+1
782 55 IF (ERR
> DBLE(NJ(K
+1))*2.D0
) GO
TO 100
784 CALL SEUL(KC
,N
,T
,Y
,DY
,FX
,FJAC
,LFJAC
,E
,LE
,IP
,H
,Ncolumns
,&
785 HmaxN
,Table
,SCAL
,NJ
,HH
,W
,A
,YHH
,DYH
,DEL
,WH
,ERR
,FacSafe1
,&
786 FAC
,FAC1
,FAC2
,FacSafe2
,Theta
,Nfun
,Ndec
,Nsol
,&
787 ERROLD
,IPHES
,ICOMP
,AUTNMS
,REJECT
,&
788 ATOV
,FSAFE
,Ncolumns2
,NRD
,IOUT
,IPT
,CALHES
)
790 IF (ERR
> 1.D0
) GO
TO 100
791 !Adi IF ((ERR > 1.D0).and.(H.gt.Hmin)) GO TO 100
792 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
793 !~~~> STEP IS ACCEPTED
794 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
806 SCAL(i
)=AbsTol(1)+RelTol(1)*DABS(T1I
)
808 SCAL(i
)=AbsTol(i
)+RelTol(i
)*DABS(T1I
)
818 DENS(NRD
+I
)=Y(ICOMP(i
))
821 !~~~> COMPUTE DIFFERENCES
828 FSAFE(L
,I
)=FSAFE(L
,I
)-FSAFE(L
-1,I
)
833 !~~~> COMPUTE DERIVATIVES AT RIGHT END ----
836 FACNJ
=FACNJ
**KLR
/FACUL(KLR
+1)
839 KRN
=(KK
-LAMBDA
+1)*NRD
840 DENS(KRN
+I
)=FSAFE(IPT
,I
)*FACNJ
845 DO L
=J
,KLR
+LAMBDA
+1,-1
846 FACTOR
=DBLENJ
/NJ(L
-1)-1.D0
848 KRN
=(L
-LAMBDA
+1)*NRD
+I
849 DENS(KRN
-NRD
)=DENS(KRN
)+(DENS(KRN
)-DENS(KRN
-NRD
))/FACTOR
854 !~~~> COMPUTE THE COEFFICIENTS OF THE INTERPOLATION POLYNOMIAL
858 DENS(II
)=DENS(II
)-DENS(II
-NRD
)
862 !~~~> COMPUTE OPTIMAL ORDER
870 IF (W(KC
-1) < W(KC
)*FAC3
) KOPT
=KC
-1
871 IF (W(KC
) < W(KC
-1)*FAC4
) KOPT
=MIN0(KC
+1,Ncolumns
-1)
874 IF (KC
> 3.AND
.W(KC
-2) < W(KC
-1)*FAC3
) KOPT
=KC
-2
875 IF (W(KC
) < W(KOPT
)*FAC4
) KOPT
=MIN0(KC
,Ncolumns
-1)
877 !~~~> AFTER A REJECTED STEP
884 !~~~> COMPUTE STEP SIZE FOR NEXT STEP
888 IF (KC
< K
.AND
.W(KC
) < W(KC
-1)*FAC4
) THEN
889 H
=HH(KC
)*A(KOPT
+1)/A(KC
)
891 H
=HH(KC
)*A(KOPT
)/A(KC
)
895 !Adi H = MAX(H, Hmin)
897 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
898 !~~~> STEP IS REJECTED
899 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
901 IF (K
> 2.AND
.W(K
-1) < W(K
)*FAC3
) K
=K
-1
914 120 WRITE (6,979) T
,H
915 979 FORMAT(' EXIT OF SEULEX AT T=',D14
.7
,' H=',D14
.7
)
920 END SUBROUTINE KPP_ROOT_SEULEX_Integrator
923 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
924 SUBROUTINE KPP_ROOT_SEUL(JJ
,N
,T
,Y
,DY
,FX
,FJAC
,LFJAC
,E
,LE
,IP
,&
925 H
,Ncolumns
,HmaxN
,Table
,SCAL
,NJ
,HH
,W
,A
,YH
,DYH
,DEL
,WH
,ERR
,FacSafe1
, &
926 FAC
,FAC1
,FAC2
,FacSafe2
,Theta
,Nfun
,Ndec
,Nsol
,&
927 ERROLD
,IPHES
,ICOMP
, &
928 AUTNMS
,REJECT
,ATOV
,FSAFE
,Ncolumns2
,NRD
,IOUT
, &
930 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
931 !~~~> THIS SUBROUTINE COMPUTES THE J-TH LINE OF THE
932 !~~~> EXTRAPOLATION TABLE AND PROVIDES AN ESTIMATE
933 !~~~> OF THE OPTIMAL STEP SIZE
934 USE KPP_ROOT_Parameters
935 USE KPP_ROOT_Jacobian
936 IMPLICIT KPP_REAL (A
-H
,O
-Z
)
937 IMPLICIT INTEGER (I
-N
)
938 INTEGER :: Ncolumns
, Ncolumns2
, N
, NRD
939 KPP_REAL
:: Y(NVAR
),YH(NVAR
),DY(NVAR
),FX(NVAR
),DYH(NVAR
)
940 KPP_REAL
:: DEL(NVAR
),WH(NVAR
),SCAL(NVAR
),HH(Ncolumns
),W(Ncolumns
),A(Ncolumns
)
942 KPP_REAL
:: FJAC(NVAR
,NVAR
), E(NVAR
,NVAR
)
944 KPP_REAL
:: FJAC(LU_NONZERO
), E(LU_NONZERO
)
946 KPP_REAL
:: Table(Ncolumns
,NVAR
)
947 KPP_REAL
:: FSAFE(Ncolumns2
,NRD
)
948 INTEGER :: IP(N
),NJ(Ncolumns
),IPHES(N
),ICOMP(NRD
)
949 LOGICAL ATOV
,REJECT
,AUTNMS
,CALHES
951 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
952 ! COMPUTE THE MATRIX E AND ITS DECOMPOSITION
953 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
963 CALL DGETRF(N
,N
,E
,N
,IP
,ISING
)
969 E(LU_DIAG(j
))=E(LU_DIAG(j
))+HJI
971 CALL KppDecomp (E
,ISING
)
974 IF (ISING
.NE
.0) GOTO 79
975 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
976 !~~~> STARTING PROCEDURE
977 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
978 IF (.NOT
.AUTNMS
) THEN
979 CALL FUN_CHEM(T
+HJ
,Y
,DY
)
986 CALL DGETRS ('N',N
,1,E
,N
,IP
,DEL
,N
,ISING
)
988 CALL KppSolve (E
,DEL
)
992 IF (IOUT
== 2.AND
.M
== JJ
) THEN
995 FSAFE(IPT
,I
)=DEL(ICOMP(i
))
998 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
999 !~~~> SEMI-IMPLICIT EULER METHOD
1000 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1007 CALL FUN_CHEM(T
+HJ
*MM
,YH
,DYH
)
1009 CALL FUN_CHEM(T
+HJ
*(MM
+1),YH
,DYH
)
1012 IF (MM
== 1.AND
.JJ
<= 2) THEN
1013 !~~~> STABILITY CHECK
1016 DEL1
=DEL1
+(DEL(i
)/SCAL(i
))**2
1019 IF (.NOT
.AUTNMS
) THEN
1020 CALL FUN_CHEM(T
+HJ
,YH
,WH
)
1023 DEL(i
)=WH(i
)-DEL(i
)*HJI
1027 DEL(i
)=DYH(i
)-DEL(i
)*HJI
1031 CALL DGETRS ('N',N
,1,E
,N
,IP
,DEL
,N
,ISING
)
1033 CALL KppSolve (E
,DEL
)
1038 DEL2
=DEL2
+(DEL(i
)/SCAL(i
))**2
1041 Theta
=DEL2
/MAX(1.D0
,DEL1
)
1042 IF (Theta
> 1.D0
) GOTO 79
1045 CALL DGETRS ('N',N
,1,E
,N
,IP
,DYH
,N
,ISING
)
1047 CALL KppSolve (E
,DYH
)
1053 IF (IOUT
== 2.AND
.MM
>= M
-JJ
) THEN
1056 FSAFE(IPT
,i
)=DEL(ICOMP(i
))
1062 Table(JJ
,I
)=YH(i
)+DEL(i
)
1064 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1065 !~~~> POLYNOMIAL EXTRAPOLATION
1066 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1069 FAC
=(DBLE(NJ(JJ
))/DBLE(NJ(L
-1)))-1.D0
1071 Table(L
-1,I
)=Table(L
,I
)+(Table(L
,I
)-Table(L
-1,I
))/FAC
1076 ERR
=ERR
+MIN(ABS((Table(1,I
)-Table(2,I
)))/SCAL(i
),1.D15
)**2
1078 IF (ERR
>= 1.D30
) GOTO 79
1079 ERR
=SQRT(ERR
/DBLE(N
))
1080 IF (JJ
> 2.AND
.ERR
>= ERROLD
) GOTO 79
1081 ERROLD
=MAX(4*ERR
,1.D0
)
1082 !~~~> COMPUTE OPTIMAL STEP SIZES
1085 FAC
=MIN(FAC2
/FACMIN
,MAX(FACMIN
,(ERR
/FacSafe1
)**EXPO
/FacSafe2
))
1087 HH(JJ
)=MIN(H
*FAC
,HmaxN
)
1094 END SUBROUTINE KPP_ROOT_SEUL
1098 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1099 SUBROUTINE KPP_ROOT_FUN_CHEM( T
, V
, FCT
)
1100 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1102 USE KPP_ROOT_Parameters
1104 USE KPP_ROOT_Function
, ONLY
: Fun
1105 USE KPP_ROOT_Rates
, ONLY
: Update_SUN
, Update_RCONST
, Update_PHOTO
1109 KPP_REAL
:: V(NVAR
), FCT(NVAR
)
1115 !CALL Update_RCONST()
1116 !CALL Update_PHOTO()
1118 CALL Fun(V
, FIX
, RCONST
, FCT
)
1121 END SUBROUTINE KPP_ROOT_FUN_CHEM
1124 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1125 SUBROUTINE KPP_ROOT_JAC_CHEM ( T
, V
, Jcb
)
1126 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1128 USE KPP_ROOT_Parameters
1130 USE KPP_ROOT_Jacobian
, ONLY
: Jac_SP
1131 USE KPP_ROOT_Rates
, ONLY
: Update_SUN
, Update_RCONST
, Update_PHOTO
1135 KPP_REAL
:: V(NVAR
), T
, TOLD
1137 KPP_REAL
:: JV(LU_NONZERO
), Jcb(NVAR
,NVAR
)
1140 KPP_REAL
:: Jcb(LU_NONZERO
)
1146 !CALL Update_RCONST()
1147 !CALL Update_PHOTO()
1151 CALL Jac_SP(V
, FIX
, RCONST
, JV
)
1158 Jcb(LU_IROW(i
),LU_ICOL(i
)) = JV(i
)
1161 CALL Jac_SP(V
, FIX
, RCONST
, Jcb
)
1165 END SUBROUTINE KPP_ROOT_JAC_CHEM
1167 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1169 END MODULE KPP_ROOT_Integrator