source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int.modified_WCOPY/runge_kutta_adj.f90 @ 2696

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

Merge of branch palm4u into trunk

File size: 95.0 KB
Line 
1!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
2!  RungeKuttaADJ - Adjoint Model of Fully Implicit 3-stage Runge-Kutta:   !
3!          * Radau-2A   quadrature (order 5)                              !
4!          * Radau-1A   quadrature (order 5)                              !
5!          * Lobatto-3C quadrature (order 4)                              !
6!          * Gauss      quadrature (order 6)                              !
7!  By default the code employs the KPP sparse linear algebra routines     !
8!  Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK)        !
9!                                                                         !
10!    (C)  Adrian Sandu, August 2005                                       !
11!    Virginia Polytechnic Institute and State University                  !
12!    Contact: sandu@cs.vt.edu                                             !
13!    Revised by Philipp Miehe and Adrian Sandu, May 2006                  !
14!    This implementation is part of KPP - the Kinetic PreProcessor        !
15!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
16
17MODULE KPP_ROOT_Integrator
18
19  USE KPP_ROOT_Precision
20  USE KPP_ROOT_Parameters, ONLY: NVAR, NSPEC, NFIX, LU_NONZERO
21  USE KPP_ROOT_Global, ONLY: FIX, RCONST, TIME
22  USE KPP_ROOT_Jacobian
23  USE KPP_ROOT_LinearAlgebra
24
25  IMPLICIT NONE
26  PUBLIC
27  SAVE
28
29!~~~>  Statistics on the work performed by the Runge-Kutta method
30  INTEGER :: Nfun,Njac,Nstp,Nacc,Nrej,Ndec,Nsol,Nsng
31  INTEGER, PARAMETER :: ifun=1, ijac=2, istp=3, iacc=4, &
32    irej=5, idec=6, isol=7, isng=8, itexit=1, ihacc=2, ihnew=3
33 
34CONTAINS
35
36  ! **************************************************************************
37
38  SUBROUTINE INTEGRATE_ADJ(NADJ, Y, Lambda, TIN, TOUT,  &
39             ATOL_adj, RTOL_adj,                        &
40             ICNTRL_U, RCNTRL_U, ISTATUS_U, RSTATUS_U, IERR_U )
41
42    USE KPP_ROOT_Parameters, ONLY: NVAR
43    USE KPP_ROOT_Global,     ONLY: ATOL,RTOL
44
45    IMPLICIT NONE
46
47!~~~> Y - Concentrations
48   KPP_REAL  :: Y(NVAR)
49!~~~> NADJ - No. of cost functionals for which adjoints
50!                are evaluated simultaneously
51!            If single cost functional is considered (like in
52!                most applications) simply set NADJ = 1     
53   INTEGER NADJ
54!~~~> Lambda - Sensitivities of concentrations
55!     Note: Lambda (1:NVAR,j) contains sensitivities of
56!           the j-th cost functional w.r.t. Y(1:NVAR), j=1...NADJ
57    KPP_REAL :: Lambda(NVAR,NADJ)
58!~~~> Tolerances for adjoint calculations
59!     (used for full continuous adjoint, and for controlling
60!     the convergence of iterations for solving the discrete adjoint)   
61    KPP_REAL, INTENT(IN)  :: ATOL_adj(NVAR,NADJ), RTOL_adj(NVAR,NADJ)
62    KPP_REAL :: TIN  ! TIN - Start Time
63    KPP_REAL :: TOUT ! TOUT - End Time
64    INTEGER,       INTENT(IN),  OPTIONAL :: ICNTRL_U(20)
65    KPP_REAL, INTENT(IN),  OPTIONAL :: RCNTRL_U(20)
66    INTEGER,       INTENT(OUT), OPTIONAL :: ISTATUS_U(20)
67    KPP_REAL, INTENT(OUT), OPTIONAL :: RSTATUS_U(20)
68    INTEGER,       INTENT(OUT), OPTIONAL :: IERR_U
69
70    INTEGER :: IERR
71
72    KPP_REAL :: RCNTRL(20), RSTATUS(20), T1, T2
73    INTEGER       :: ICNTRL(20), ISTATUS(20)
74    INTEGER, SAVE :: Ntotal = 0
75
76
77    ICNTRL(1:20) = 0
78    RCNTRL(1:20) = 0.0_dp
79
80    !~~~> fine-tune the integrator:
81    ICNTRL(2)  = 0   ! 0=vector tolerances, 1=scalar tolerances
82    ICNTRL(5)  = 8   ! Max no. of Newton iterations
83    ICNTRL(6)  = 0   ! Starting values for Newton are: 0=interpolated, 1=zero
84    ICNTRL(7)  = 2   ! Adj. system solved by: 1=iteration, 2=direct, 3=adaptive
85    ICNTRL(8)  = 0   ! Adj. LU decomp: 0=compute, 1=save from fwd
86    ICNTRL(9)  = 2   ! Adjoint: 1=none, 2=discrete, 3=full continuous, 4=simplified continuous
87    ICNTRL(10) = 0   ! Error estimator: 0=classic, 1=SDIRK
88    ICNTRL(11) = 1   ! Step controller: 1=Gustaffson, 2=classic
89
90
91    !~~~> if optional parameters are given, and if they are >0,
92    !     then use them to overwrite default settings
93    IF (PRESENT(ICNTRL_U)) THEN
94      WHERE(ICNTRL_U(:) > 0) ICNTRL(:) = ICNTRL_U(:)
95    END IF
96    IF (PRESENT(RCNTRL_U)) THEN
97      WHERE(RCNTRL_U(:) > 0) RCNTRL(:) = RCNTRL_U(:)
98    END IF
99
100   
101    T1 = TIN; T2 = TOUT
102    CALL RungeKuttaADJ(NVAR, Y, NADJ, Lambda, T1, T2,   &
103                       RTOL, ATOL, ATOL_adj, RTOL_adj,  &
104                       RCNTRL, ICNTRL, RSTATUS, ISTATUS, IERR  )
105
106    Ntotal = Ntotal + ISTATUS(istp)
107!    PRINT*,'NSTEPS=',ISTATUS(istp),' (',Ntotal,')','  O3=', VAR(ind_O3)
108
109    ! if optional parameters are given for output
110    ! use them to store information in them
111    IF (PRESENT(ISTATUS_U)) ISTATUS_U(:) = ISTATUS(:)
112    IF (PRESENT(RSTATUS_U)) RSTATUS_U(:) = RSTATUS(:)
113    IF (PRESENT(IERR_U)) IERR_U = IERR
114
115    IF (IERR < 0) THEN
116      PRINT *,'Runge-Kutta-ADJ: Unsuccessful exit at T=', TIN,' (IERR=',IERR,')'
117    ENDIF
118
119  END SUBROUTINE INTEGRATE_ADJ
120
121
122!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
123  SUBROUTINE RungeKuttaADJ( N, Y, NADJ, Lambda,Tstart,Tend,      &
124                         RelTol, AbsTol, RelTol_adj, AbsTol_adj, &
125                         RCNTRL, ICNTRL, RSTATUS, ISTATUS, IERR )
126!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
127!
128!  This implementation is based on the book and the code Radau5:
129!
130!         E. HAIRER AND G. WANNER
131!         "SOLVING ORDINARY DIFFERENTIAL EQUATIONS II.
132!              STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS."
133!         SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS 14,
134!         SPRINGER-VERLAG (1991)
135!
136!         UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES
137!         CH-1211 GENEVE 24, SWITZERLAND
138!         E-MAIL:  HAIRER@DIVSUN.UNIGE.CH,  WANNER@DIVSUN.UNIGE.CH
139!
140!   Methods:
141!          * Radau-2A   quadrature (order 5)                             
142!          * Radau-1A   quadrature (order 5)                             
143!          * Lobatto-3C quadrature (order 4)                             
144!          * Gauss      quadrature (order 6)                             
145!                                                                         
146!   (C)  Adrian Sandu, August 2005                                       
147!   Virginia Polytechnic Institute and State University                 
148!   Contact: sandu@cs.vt.edu                                             
149!   Revised by Philipp Miehe and Adrian Sandu, May 2006                 
150!   This implementation is part of KPP - the Kinetic PreProcessor       
151!
152!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
153!
154!~~~>   INPUT ARGUMENTS:
155!       ----------------
156!
157!    Note: For input parameters equal to zero the default values of the
158!          corresponding variables are used.
159!
160!     N           Dimension of the system
161!     T           Initial time value
162!
163!     Tend        Final T value (Tend-T may be positive or negative)
164!
165!     Y(N)        Initial values for Y
166!
167!     RelTol,AbsTol   Relative and absolute error tolerances.
168!          for ICNTRL(2) = 0: AbsTol, RelTol are N-dimensional vectors
169!                        = 1: AbsTol, RelTol are scalars
170!
171!~~~>  Integer input parameters:
172
173!    ICNTRL(1) = not used
174!
175!    ICNTRL(2) = 0: AbsTol, RelTol are NVAR-dimensional vectors
176!              = 1: AbsTol, RelTol are scalars
177!
178!    ICNTRL(3) = RK method selection       
179!              = 1:  Radau-2A    (the default)
180!              = 2:  Lobatto-3C
181!              = 3:  Gauss
182!              = 4:  Radau-1A
183!
184!    ICNTRL(4)  -> maximum number of integration steps
185!        For ICNTRL(4)=0 the default value of 10000 is used
186!
187!    ICNTRL(5)  -> maximum number of Newton iterations
188!        For ICNTRL(5)=0 the default value of 8 is used
189!
190!    ICNTRL(6)  -> starting values of Newton iterations:
191!        ICNTRL(6)=0 : starting values are obtained from
192!                      the extrapolated collocation solution
193!                      (the default)
194!        ICNTRL(6)=1 : starting values are zero
195!
196!    ICNTRL(7)  -> method to solve the linear ADJ equations:
197!        ICNTRL(7)=0,1 : modified Newton re-using LU (the default)
198!                        with a fixed number of iterations
199!        ICNTRL(7)=2 :   direct solution (additional one LU factorization
200!                        of 3Nx3N matrix per step); good for debugging
201!        ICNTRL(7)=3 :   adaptive solution (if Newton does not converge
202!                        switch to direct)
203!
204!    ICNTRL(8)  -> checkpointing the LU factorization at each step:
205!        ICNTRL(8)=0 : do *not* save LU factorization (the default)
206!        ICNTRL(8)=1 : save LU factorization
207!        Note: if ICNTRL(7)=1 the LU factorization is *not* saved
208!
209!    ICNTRL(9) -> Type of adjoint algorithm
210!         = 0 : default is discrete adjoint ( of method ICNTRL(3) )
211!         = 1 : no adjoint       
212!         = 2 : discrete adjoint ( of method ICNTRL(3) )
213!         = 3 : fully adaptive continuous adjoint ( with method ICNTRL(6) )
214!         = 4 : simplified continuous adjoint ( with method ICNTRL(6) )
215!
216!    ICNTRL(10) -> switch for error estimation strategy
217!               ICNTRL(10) = 0: one additional stage at c=0,
218!                               see Hairer (default)
219!               ICNTRL(10) = 1: two additional stages at c=0
220!                               and SDIRK at c=1, stiffly accurate
221!
222!    ICNTRL(11) -> switch for step size strategy
223!              ICNTRL(11)=1:  mod. predictive controller (Gustafsson, default)
224!              ICNTRL(11)=2:  classical step size control
225!              the choice 1 seems to produce safer results;
226!              for simple problems, the choice 2 produces
227!              often slightly faster runs
228!
229!~~~>  Real input parameters:
230!
231!    RCNTRL(1)  -> Hmin, lower bound for the integration step size
232!                  (highly recommended to keep Hmin = ZERO, the default)
233!
234!    RCNTRL(2)  -> Hmax, upper bound for the integration step size
235!
236!    RCNTRL(3)  -> Hstart, the starting step size
237!
238!    RCNTRL(4)  -> FacMin, lower bound on step decrease factor (default=0.2)
239!
240!    RCNTRL(5)  -> FacMax, upper bound on step increase factor (default=6)
241!
242!    RCNTRL(6)  -> FacRej, step decrease factor after multiple rejections
243!                 (default=0.1)
244!
245!    RCNTRL(7)  -> FacSafe, by which the new step is slightly smaller
246!                  than the predicted value  (default=0.9)
247!
248!    RCNTRL(8)  -> ThetaMin. If Newton convergence rate smaller
249!                  than ThetaMin the Jacobian is not recomputed;
250!                  (default=0.001)
251!
252!    RCNTRL(9)  -> NewtonTol, stopping criterion for Newton's method
253!                  (default=0.03)
254!
255!    RCNTRL(10) -> Qmin
256!
257!    RCNTRL(11) -> Qmax. If Qmin < Hnew/Hold < Qmax, then the
258!                  step size is kept constant and the LU factorization
259!                  reused (default Qmin=1, Qmax=1.2)
260!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
261!
262!
263!    OUTPUT ARGUMENTS:
264!    -----------------
265!
266!    T           -> T value for which the solution has been computed
267!                     (after successful return T=Tend).
268!
269!    Y(N)        -> Numerical solution at T
270!
271!    IERR        -> Reports on successfulness upon return:
272!                    = 1 for success
273!                    < 0 for error (value equals error code)
274!
275!    ISTATUS(1)  -> No. of function calls
276!    ISTATUS(2)  -> No. of Jacobian calls
277!    ISTATUS(3)  -> No. of steps
278!    ISTATUS(4)  -> No. of accepted steps
279!    ISTATUS(5)  -> No. of rejected steps (except at very beginning)
280!    ISTATUS(6)  -> No. of LU decompositions
281!    ISTATUS(7)  -> No. of forward/backward substitutions
282!    ISTATUS(8)  -> No. of singular matrix decompositions
283!
284!    RSTATUS(1)  -> Texit, the time corresponding to the
285!                     computed Y upon return
286!    RSTATUS(2)  -> Hexit, last accepted step before exit
287!    RSTATUS(3)  -> Hnew, last predicted step (not yet taken)
288!                   For multiple restarts, use Hnew as Hstart
289!                     in the subsequent run
290!
291!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
292
293      IMPLICIT NONE
294     
295      INTEGER :: N
296      INTEGER, INTENT(IN)     :: NADJ
297      KPP_REAL, INTENT(INOUT) :: Lambda(N,NADJ)
298      KPP_REAL, INTENT(IN) :: AbsTol_adj(N,NADJ), RelTol_adj(N,NADJ)
299      KPP_REAL :: Y(N),AbsTol(N),RelTol(N),RCNTRL(20),RSTATUS(20)
300      INTEGER :: ICNTRL(20), ISTATUS(20)
301      LOGICAL :: StartNewton, Gustafsson, SdirkError, SaveLU
302      INTEGER :: IERR, ITOL
303      KPP_REAL, INTENT(IN):: Tstart,Tend
304      KPP_REAL :: Texit
305
306      !~~~> Control arguments
307      INTEGER :: Max_no_steps, NewtonMaxit, AdjointType, rkMethod
308      KPP_REAL :: Hmin,Hmax,Hstart,Qmin,Qmax
309      KPP_REAL :: Roundoff, ThetaMin, NewtonTol
310      KPP_REAL :: FacSafe,FacMin,FacMax,FacRej
311      ! Runge-Kutta method parameters
312      INTEGER, PARAMETER :: R2A=1, R1A=2, L3C=3, GAU=4, L3A=5
313      KPP_REAL :: rkT(3,3), rkTinv(3,3), rkTinvAinv(3,3), rkAinvT(3,3), &
314                       rkA(3,3), rkB(3),  rkC(3), rkD(0:3), rkE(0:3),      &
315                       rkBgam(0:4), rkBhat(0:4), rkTheta(0:3),             &
316                       rkGamma,  rkAlpha, rkBeta, rkELO
317      ! ADJ method parameters
318      INTEGER, PARAMETER :: KPP_ROOT_none = 1, KPP_ROOT_discrete = 2,              &
319                            KPP_ROOT_continuous = 3, KPP_ROOT_simple_continuous = 4
320      INTEGER :: AdjointSolve                     
321      INTEGER, PARAMETER :: Solve_direct = 1, Solve_fixed = 2,             &
322                            Solve_adaptive = 3
323      INTEGER, PARAMETER :: bufsize = 10000
324      INTEGER :: stack_ptr = 0 ! last written entry
325      KPP_REAL, DIMENSION(:), POINTER :: chk_H, chk_T
326      KPP_REAL, DIMENSION(:,:), POINTER :: chk_Y, chk_Z, chk_E1
327      INTEGER, DIMENSION(:), POINTER :: chk_NiT
328      COMPLEX(kind=dp), DIMENSION(:,:), POINTER :: chk_E2
329      KPP_REAL, DIMENSION(:,:), POINTER :: chk_dY, chk_d2Y
330      !~~~> Local variables
331      INTEGER :: i 
332      KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0
333   
334!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
335!        SETTING THE PARAMETERS
336!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
337      IERR = 0
338      ISTATUS(1:20) = 0
339      RSTATUS(1:20) = ZERO
340       
341!~~~> ICNTRL(1) - autonomous system - not used       
342!~~~> ITOL: 1 for vector and 0 for scalar AbsTol/RelTol
343      IF (ICNTRL(2) == 0) THEN
344         ITOL = 1
345      ELSE
346         ITOL = 0
347      END IF
348!~~~> Error control selection 
349      IF (ICNTRL(10) == 0) THEN
350         SdirkError = .FALSE.
351      ELSE
352         SdirkError = .TRUE.
353      END IF     
354!~~~> Method selection 
355      SELECT CASE (ICNTRL(3))     
356      CASE (0,1)
357         CALL Radau2A_Coefficients
358      CASE (2)
359         CALL Lobatto3C_Coefficients
360      CASE (3)
361         CALL Gauss_Coefficients
362      CASE (4)
363         CALL Radau1A_Coefficients
364      CASE DEFAULT
365         WRITE(6,*) 'ICNTRL(3)=',ICNTRL(3)
366         CALL RK_ErrorMsg(-13,Tstart,ZERO,IERR)
367      END SELECT
368!~~~> Max_no_steps: the maximal number of time steps
369      IF (ICNTRL(4) == 0) THEN
370         Max_no_steps = 200000
371      ELSE
372         Max_no_steps=ICNTRL(4)
373         IF (Max_no_steps <= 0) THEN
374            WRITE(6,*) 'ICNTRL(4)=',ICNTRL(4)
375            CALL RK_ErrorMsg(-1,Tstart,ZERO,IERR)
376         END IF
377      END IF
378!~~~> NewtonMaxit    maximal number of Newton iterations
379      IF (ICNTRL(5) == 0) THEN
380         NewtonMaxit = 8
381      ELSE
382         NewtonMaxit=ICNTRL(5)
383         IF (NewtonMaxit <= 0) THEN
384            WRITE(6,*) 'ICNTRL(5)=',ICNTRL(5)
385            CALL RK_ErrorMsg(-2,Tstart,ZERO,IERR)
386          END IF
387      END IF
388!~~~> StartNewton:  Use extrapolation for starting values of Newton iterations
389      IF (ICNTRL(6) == 0) THEN
390         StartNewton = .TRUE.
391      ELSE
392         StartNewton = .FALSE.
393      END IF     
394!~~~>  How to solve the linear adjoint system
395      SELECT CASE (ICNTRL(7))     
396      CASE (0,1)
397         AdjointSolve = Solve_fixed
398      CASE (2)
399         AdjointSolve = Solve_direct
400      CASE (3)
401         AdjointSolve = Solve_adaptive
402      CASE DEFAULT 
403         PRINT * , 'User-selected adjoint solution: ICNTRL(7)=', ICNTRL(7)
404         PRINT * , 'Implemented: =(0,1) (fixed), =2 (direct), =3 (adaptive)'
405         CALL rk_ErrorMsg(-9,Tstart,ZERO,IERR)
406         RETURN     
407      END SELECT
408!~~~>  Discrete or continuous adjoint formulation
409      SELECT CASE (ICNTRL(9))     
410      CASE (0,2)
411         AdjointType = KPP_ROOT_discrete
412      CASE (1)
413         AdjointType = KPP_ROOT_none
414      CASE DEFAULT 
415         PRINT * , 'User-selected adjoint type: ICNTRL(9)=', ICNTRL(9)
416         PRINT * , 'Implemented: =(0,2) (discrete adj) and =1 (no adj)'
417         CALL rk_ErrorMsg(-9,Tstart,ZERO,IERR)
418         RETURN     
419      END SELECT
420!~~~> Save or not the forward LU factorization
421      SaveLU = (ICNTRL(8) /= 0) 
422      IF (AdjointSolve == Solve_direct) SaveLU = .FALSE.
423!~~~> Gustafsson: step size controller
424      IF (ICNTRL(11) == 0) THEN
425         Gustafsson = .TRUE.
426      ELSE
427         Gustafsson = .FALSE.
428      END IF
429
430!~~~> Roundoff: smallest number s.t. 1.0 + Roundoff > 1.0
431      Roundoff = WLAMCH('E');
432
433!~~~> Hmin = minimal step size
434      IF (RCNTRL(1) == ZERO) THEN
435         Hmin = ZERO
436      ELSE
437         Hmin = MIN(ABS(RCNTRL(1)),ABS(Tend-Tstart))
438      END IF
439!~~~> Hmax = maximal step size
440      IF (RCNTRL(2) == ZERO) THEN
441         Hmax = ABS(Tend-Tstart)
442      ELSE
443         Hmax = MIN(ABS(RCNTRL(2)),ABS(Tend-Tstart))
444      END IF
445!~~~> Hstart = starting step size
446      IF (RCNTRL(3) == ZERO) THEN
447         Hstart = ZERO
448      ELSE
449         Hstart = MIN(ABS(RCNTRL(3)),ABS(Tend-Tstart))
450      END IF
451!~~~> FacMin: lower bound on step decrease factor
452      IF(RCNTRL(4) == ZERO)THEN
453         FacMin = 0.2d0
454      ELSE
455         FacMin = RCNTRL(4)
456      END IF
457!~~~> FacMax: upper bound on step increase factor
458      IF(RCNTRL(5) == ZERO)THEN
459         FacMax = 8.D0
460      ELSE
461         FacMax = RCNTRL(5)
462      END IF
463!~~~> FacRej: step decrease factor after 2 consecutive rejections
464      IF(RCNTRL(6) == ZERO)THEN
465         FacRej = 0.1d0
466      ELSE
467         FacRej = RCNTRL(6)
468      END IF
469!~~~> FacSafe:  by which the new step is slightly smaller
470!               than the predicted value
471      IF (RCNTRL(7) == ZERO) THEN
472         FacSafe=0.9d0
473      ELSE
474         FacSafe=RCNTRL(7)
475      END IF
476      IF ( (FacMax < ONE) .OR. (FacMin > ONE) .OR. &
477           (FacSafe <= 1.0d-3) .OR. (FacSafe >= ONE) ) THEN
478            WRITE(6,*)'RCNTRL(4:7)=',RCNTRL(4:7)
479            CALL RK_ErrorMsg(-4,Tstart,ZERO,IERR)
480      END IF
481
482!~~~> ThetaMin:  decides whether the Jacobian should be recomputed
483      IF (RCNTRL(8) == ZERO) THEN
484         ThetaMin = 1.0d-3
485      ELSE
486         ThetaMin=RCNTRL(8)
487         IF (ThetaMin <= 0.0d0 .OR. ThetaMin >= 1.0d0) THEN
488            WRITE(6,*) 'RCNTRL(8)=', RCNTRL(8)
489            CALL RK_ErrorMsg(-5,Tstart,ZERO,IERR)
490         END IF
491      END IF
492!~~~> NewtonTol:  stopping crierion for Newton's method
493      IF (RCNTRL(9) == ZERO) THEN
494         NewtonTol = 3.0d-2
495      ELSE
496         NewtonTol = RCNTRL(9)
497         IF (NewtonTol <= Roundoff) THEN
498            WRITE(6,*) 'RCNTRL(9)=',RCNTRL(9)
499            CALL RK_ErrorMsg(-6,Tstart,ZERO,IERR)
500         END IF
501      END IF
502!~~~> Qmin AND Qmax: IF Qmin < Hnew/Hold < Qmax then step size = const.
503      IF (RCNTRL(10) == ZERO) THEN
504         Qmin=1.D0
505      ELSE
506         Qmin=RCNTRL(10)
507      END IF
508      IF (RCNTRL(11) == ZERO) THEN
509         Qmax=1.2D0
510      ELSE
511         Qmax=RCNTRL(11)
512      END IF
513      IF (Qmin > ONE .OR. Qmax < ONE) THEN
514         WRITE(6,*) 'RCNTRL(10:11)=',Qmin,Qmax
515         CALL RK_ErrorMsg(-7,Tstart,ZERO,IERR)
516      END IF
517!~~~> Check if tolerances are reasonable
518      IF (ITOL == 0) THEN
519          IF (AbsTol(1) <= ZERO.OR.RelTol(1) <= 10.d0*Roundoff) THEN
520              WRITE (6,*) 'AbsTol/RelTol=',AbsTol,RelTol 
521              CALL RK_ErrorMsg(-8,Tstart,ZERO,IERR)
522          END IF
523      ELSE
524          DO i=1,N
525          IF (AbsTol(i) <= ZERO.OR.RelTol(i) <= 10.d0*Roundoff) THEN
526              WRITE (6,*) 'AbsTol/RelTol(',i,')=',AbsTol(i),RelTol(i)
527              CALL RK_ErrorMsg(-8,Tstart,ZERO,IERR)
528          END IF
529          END DO
530      END IF
531
532!~~~> Parameters are wrong
533      IF (IERR < 0) RETURN
534
535!~~~>  Allocate checkpoint space or open checkpoint files
536   IF (AdjointType == KPP_ROOT_discrete) THEN
537       CALL rk_AllocateDBuffers()
538   ELSEIF ( (AdjointType == KPP_ROOT_continuous).OR. &
539           (AdjointType == KPP_ROOT_simple_continuous) ) THEN
540       CALL rk_AllocateCBuffers
541   END IF
542
543!~~~> Call the core method
544      CALL RK_FwdIntegrator( N,Tstart,Tend,Y,AdjointType,IERR )
545!   PRINT*,'FORWARD STATISTICS'
546!   PRINT*,'Step=',Istatus(istp),' Acc=',Istatus(iacc),   &
547!        ' Rej=',Istatus(irej), ' Singular=',Istatus(isng)
548   Nstp = 0
549   Nacc = 0
550   Nrej = 0
551   Nsng = 0
552
553!~~~>  If Forward integration failed return   
554   IF (IERR<0) RETURN
555
556   SELECT CASE (AdjointType)   
557   CASE (KPP_ROOT_discrete)   
558     CALL rk_DadjInt (                          &
559        NADJ, Lambda,                           &
560        Tstart, Tend, Texit,                    &
561        IERR )
562   CASE (KPP_ROOT_continuous) 
563     CALL rk_CadjInt (                          &
564        NADJ, Lambda,                           &
565        Tend, Tstart, Texit,                    &
566        IERR )
567   CASE (KPP_ROOT_simple_continuous)
568     CALL rk_SimpleCadjInt (                    &
569        NADJ, Lambda,                           &
570        Tstart, Tend, Texit,                    &
571        IERR )
572   END SELECT ! AdjointType
573!   PRINT*,'ADJOINT STATISTICS'
574!   PRINT*,'Step=',Nstp,' Acc=',Nacc,             &
575!        ' Rej=',Nrej, ' Singular=',Nsng
576
577!~~~>  Free checkpoint space or close checkpoint files
578   IF (AdjointType == KPP_ROOT_discrete) THEN
579      CALL rk_FreeDBuffers
580   ELSEIF ( (AdjointType == KPP_ROOT_continuous) .OR. &
581           (AdjointType == KPP_ROOT_simple_continuous) ) THEN
582      CALL rk_FreeCBuffers
583   END IF
584
585
586!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
587CONTAINS ! Internal procedures to RungeKuttaADJ
588!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
589
590
591!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
592  SUBROUTINE rk_AllocateDBuffers()
593!~~~>  Allocate buffer space for discrete adjoint
594!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
595   INTEGER :: i
596   
597   ALLOCATE( chk_H(bufsize), STAT=i )
598   IF (i/=0) THEN
599      PRINT*,'Failed allocation of buffer H'; STOP
600   END IF   
601   ALLOCATE( chk_T(bufsize), STAT=i )
602   IF (i/=0) THEN
603      PRINT*,'Failed allocation of buffer T'; STOP
604   END IF   
605   ALLOCATE( chk_Y(NVAR,bufsize), STAT=i )
606   IF (i/=0) THEN
607      PRINT*,'Failed allocation of buffer Y'; STOP
608   END IF   
609   ALLOCATE( chk_Z(NVAR*3,bufsize), STAT=i )
610   IF (i/=0) THEN
611      PRINT*,'Failed allocation of buffer Z'; STOP
612   END IF   
613   ALLOCATE( chk_NiT(bufsize), STAT=i )
614   IF (i/=0) THEN
615      PRINT*,'Failed allocation of buffer NiT'; STOP
616   END IF   
617   IF (SaveLU) THEN 
618#ifdef FULL_ALGEBRA     
619     ALLOCATE( chk_E1(NVAR*NVAR,bufsize), STAT=i )
620     IF (i/=0) THEN
621        PRINT*,'Failed allocation of buffer E1'; STOP
622     END IF   
623     ALLOCATE( chk_E2(NVAR*NVAR,bufsize), STAT=i )
624     IF (i/=0) THEN
625        PRINT*,'Failed allocation of buffer E2'; STOP
626     END IF   
627#else     
628     ALLOCATE( chk_E1(LU_NONZERO,bufsize), STAT=i )
629     IF (i/=0) THEN
630        PRINT*,'Failed allocation of buffer E1'; STOP
631     END IF   
632     ALLOCATE( chk_E2(LU_NONZERO,bufsize), STAT=i )
633     IF (i/=0) THEN
634        PRINT*,'Failed allocation of buffer E2'; STOP
635     END IF   
636#endif
637   END IF   
638   
639 
640 END SUBROUTINE rk_AllocateDBuffers
641
642
643!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
644 SUBROUTINE rk_FreeDBuffers()
645!~~~>  Dallocate buffer space for discrete adjoint
646!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
647   INTEGER :: i
648   
649   DEALLOCATE( chk_H, STAT=i )
650   IF (i/=0) THEN
651      PRINT*,'Failed deallocation of buffer H'; STOP
652   END IF   
653   DEALLOCATE( chk_T, STAT=i )
654   IF (i/=0) THEN
655      PRINT*,'Failed deallocation of buffer T'; STOP
656   END IF   
657   DEALLOCATE( chk_Y, STAT=i )
658   IF (i/=0) THEN
659      PRINT*,'Failed deallocation of buffer Y'; STOP
660   END IF   
661   DEALLOCATE( chk_Z, STAT=i )
662   IF (i/=0) THEN
663      PRINT*,'Failed deallocation of buffer Z'; STOP
664   END IF   
665   DEALLOCATE( chk_NiT, STAT=i )
666   IF (i/=0) THEN
667      PRINT*,'Failed deallocation of buffer NiT'; STOP
668   END IF   
669   IF (SaveLU) THEN
670     DEALLOCATE( chk_E1, STAT=i )
671     IF (i/=0) THEN
672        PRINT*,'Failed allocation of buffer E1'; STOP
673     END IF   
674     DEALLOCATE( chk_E2, STAT=i )
675     IF (i/=0) THEN
676        PRINT*,'Failed allocation of buffer E2'; STOP
677     END IF   
678   END IF   
679
680 END SUBROUTINE rk_FreeDBuffers
681
682
683!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
684 SUBROUTINE rk_AllocateCBuffers()
685!~~~>  Allocate buffer space for continuous adjoint
686!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
687   INTEGER :: i
688   
689   ALLOCATE( chk_H(bufsize), STAT=i )
690   IF (i/=0) THEN
691      PRINT*,'Failed allocation of buffer H'; STOP
692   END IF   
693   ALLOCATE( chk_T(bufsize), STAT=i )
694   IF (i/=0) THEN
695      PRINT*,'Failed allocation of buffer T'; STOP
696   END IF   
697   ALLOCATE( chk_Y(NVAR,bufsize), STAT=i )
698   IF (i/=0) THEN
699      PRINT*,'Failed allocation of buffer Y'; STOP
700   END IF   
701   ALLOCATE( chk_dY(NVAR,bufsize), STAT=i )
702   IF (i/=0) THEN
703      PRINT*,'Failed allocation of buffer dY'; STOP
704   END IF   
705   ALLOCATE( chk_d2Y(NVAR,bufsize), STAT=i )
706   IF (i/=0) THEN
707      PRINT*,'Failed allocation of buffer d2Y'; STOP
708   END IF   
709 
710 END SUBROUTINE rk_AllocateCBuffers
711
712
713!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
714 SUBROUTINE rk_FreeCBuffers()
715!~~~>  Dallocate buffer space for continuous adjoint
716!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
717   INTEGER :: i
718   print*,'cbuffers deallocate???'
719   DEALLOCATE( chk_H, STAT=i )
720   IF (i/=0) THEN
721      PRINT*,'Failed deallocation of buffer H'; STOP
722   END IF   
723   DEALLOCATE( chk_T, STAT=i )
724   IF (i/=0) THEN
725      PRINT*,'Failed deallocation of buffer T'; STOP
726   END IF   
727   DEALLOCATE( chk_Y, STAT=i )
728   IF (i/=0) THEN
729      PRINT*,'Failed deallocation of buffer Y'; STOP
730   END IF   
731   DEALLOCATE( chk_dY, STAT=i )
732   IF (i/=0) THEN
733      PRINT*,'Failed deallocation of buffer dY'; STOP
734   END IF   
735   DEALLOCATE( chk_d2Y, STAT=i )
736   IF (i/=0) THEN
737      PRINT*,'Failed deallocation of buffer d2Y'; STOP
738   END IF   
739 
740 END SUBROUTINE rk_FreeCBuffers
741
742!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
743 SUBROUTINE rk_DPush( T, H, Y, Zstage, NewIt, E1, E2 )!, Jcb )
744!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
745!~~~> Saves the next trajectory snapshot for discrete adjoints
746!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
747   KPP_REAL :: T, H, Y(NVAR), Zstage(NVAR*3)
748   INTEGER :: NewIt
749#ifdef FULL_ALGEBRA     
750   KPP_REAL :: E1(NVAR,NVAR)
751   COMPLEX(kind=dp) :: E2(NVAR,NVAR)
752   INTEGER :: i, j   
753#else     
754   KPP_REAL :: E1(LU_NONZERO)
755   COMPLEX(kind=dp) :: E2(LU_NONZERO)   
756#endif     
757
758   stack_ptr = stack_ptr + 1
759   IF ( stack_ptr > bufsize ) THEN
760     PRINT*,'Push failed: buffer overflow'
761     STOP
762   END IF 
763   chk_H( stack_ptr ) = H
764   chk_T( stack_ptr ) = T
765   ! CALL WCOPY(NVAR,Y,1,chk_Y(1,stack_ptr),1)
766   ! CALL WCOPY(NVAR*3,Zstage,1,chk_Z(1,stack_ptr),1)
767   chk_Y(1:N,stack_ptr) = Y(1:N)
768   chk_Z(1:3*N,stack_ptr) = Zstage(1:3*N)
769   chk_NiT( stack_ptr ) = NewIt
770   IF (SaveLU) THEN 
771#ifdef FULL_ALGEBRA     
772     ! CALL WCOPY(NVAR*NVAR, E1,1,chk_E1(1,stack_ptr),1)
773     ! CALL WCOPYCmplx(NVAR*NVAR, E2,1,chk_E2(1,stack_ptr),1)
774     DO j=1,NVAR
775       DO i=1,NVAR
776         chk_E1(NVAR*(j-1)+i,stack_ptr) = E1(i,j)
777         chk_E2(NVAR*(j-1)+i,stack_ptr) = E2(i,j)
778       END DO
779     END DO
780#else     
781     ! CALL WCOPY(LU_NONZERO, E1,1,chk_E1(1,stack_ptr),1)
782     ! CALL WCOPYCmplx(LU_NONZERO, E2,1,chk_E2(1,stack_ptr),1)
783     chk_E1(1:LU_NONZERO,stack_ptr) = E1(1:LU_NONZERO)
784     chk_E2(1:LU_NONZERO,stack_ptr) = E2(1:LU_NONZERO)
785#endif     
786   END IF 
787   !CALL WCOPY(LU_NONZERO,Jcb,1,chk_J(1,stack_ptr),1)
788 
789  END SUBROUTINE rk_DPush
790 
791   
792!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
793  SUBROUTINE rk_DPop( T, H, Y, Zstage, NewIt, E1, E2 ) !, Jcb )
794!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
795!~~~> Retrieves the next trajectory snapshot for discrete adjoints
796!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
797
798   KPP_REAL :: T, H, Y(NVAR), Zstage(NVAR*3) ! , Jcb(LU_NONZERO)
799   INTEGER :: NewIt
800#ifdef FULL_ALGEBRA     
801   KPP_REAL :: E1(NVAR,NVAR)
802   COMPLEX(kind=dp) :: E2(NVAR,NVAR)   
803   INTEGER :: i, j   
804#else     
805   KPP_REAL :: E1(LU_NONZERO)
806   COMPLEX(kind=dp) :: E2(LU_NONZERO)   
807#endif     
808   
809   IF ( stack_ptr <= 0 ) THEN
810     PRINT*,'Pop failed: empty buffer'
811     STOP
812   END IF 
813   H = chk_H( stack_ptr )
814   T = chk_T( stack_ptr )
815   ! CALL WCOPY(NVAR,chk_Y(1,stack_ptr),1,Y,1)
816   Y(1:NVAR) = chk_Y(1:NVAR,stack_ptr)
817   ! CALL WCOPY(NVAR*3,chk_Z(1,stack_ptr),1,Zstage,1)
818   Zstage(1:3*NVAR) = chk_Z(1:3*NVAR,stack_ptr)
819   NewIt = chk_NiT( stack_ptr )
820   IF (SaveLU) THEN
821#ifdef FULL_ALGEBRA     
822     ! CALL WCOPY(NVAR*NVAR,chk_E1(1,stack_ptr),1, E1,1)
823     ! CALL WCOPYCmplx(NVAR*NVAR,chk_E2(1,stack_ptr),1, E2,1)
824     DO j=1,NVAR
825       DO i=1,NVAR
826         E1(i,j) = chk_E1(NVAR*(j-1)+i,stack_ptr)
827         E2(i,j) = chk_E2(NVAR*(j-1)+i,stack_ptr)
828       END DO
829     END DO
830#else     
831     ! CALL WCOPY(LU_NONZERO,chk_E1(1,stack_ptr),1, E1,1)
832     ! CALL WCOPYCmplx(LU_NONZERO,chk_E2(1,stack_ptr),1, E2,1)
833     E1(1:LU_NONZERO) = chk_E1(1:LU_NONZERO,stack_ptr)
834     E2(1:LU_NONZERO) = chk_E2(1:LU_NONZERO,stack_ptr)
835#endif     
836   END IF 
837   !CALL WCOPY(LU_NONZERO,chk_J(1,stack_ptr),1,Jcb,1)
838
839   stack_ptr = stack_ptr - 1
840 
841  END SUBROUTINE rk_DPop
842 
843!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
844  SUBROUTINE rk_CPush(T, H, Y, dY, d2Y )
845!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
846!~~~> Saves the next trajectory snapshot for discrete adjoints
847!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
848
849   KPP_REAL :: T, H, Y(NVAR), dY(NVAR), d2Y(NVAR)
850   
851   stack_ptr = stack_ptr + 1
852   IF ( stack_ptr > bufsize ) THEN
853     PRINT*,'Push failed: buffer overflow'
854     STOP
855   END IF 
856   chk_H( stack_ptr ) = H
857   chk_T( stack_ptr ) = T
858   ! CALL WCOPY(NVAR,Y,1,chk_Y(1,stack_ptr),1)
859   ! CALL WCOPY(NVAR,dY,1,chk_dY(1,stack_ptr),1)
860   ! CALL WCOPY(NVAR,d2Y,1,chk_d2Y(1,stack_ptr),1)
861   chk_Y(1:NVAR,stack_ptr)  =  Y(1:NVAR)
862   chk_dY(1:NVAR,stack_ptr) = dY(1:NVAR)
863   chk_d2Y(1:NVAR,stack_ptr)= d2Y(1:NVAR)
864 
865  END SUBROUTINE rk_CPush
866 
867   
868!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
869  SUBROUTINE rk_CPop( T, H, Y, dY, d2Y )
870!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
871!~~~> Retrieves the next trajectory snapshot for discrete adjoints
872!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
873
874   KPP_REAL :: T, H, Y(NVAR), dY(NVAR), d2Y(NVAR)
875   
876   IF ( stack_ptr <= 0 ) THEN
877     PRINT*,'Pop failed: empty buffer'
878     STOP
879   END IF 
880   H = chk_H( stack_ptr )
881   T = chk_T( stack_ptr )
882   ! CALL WCOPY(NVAR,chk_Y(1,stack_ptr),1,Y,1)
883   ! CALL WCOPY(NVAR,chk_dY(1,stack_ptr),1,dY,1)
884   ! CALL WCOPY(NVAR,chk_d2Y(1,stack_ptr),1,d2Y,1)
885     Y(1:NVAR) = chk_Y(1:NVAR,stack_ptr)
886    dY(1:NVAR) = chk_dY(1:NVAR,stack_ptr)
887   d2Y(1:NVAR) = chk_d2Y(1:NVAR,stack_ptr)
888
889   stack_ptr = stack_ptr - 1
890 
891  END SUBROUTINE rk_CPop
892
893
894!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
895   SUBROUTINE RK_FwdIntegrator( N,Tstart,Tend,Y,AdjointType,IERR )
896!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
897
898      IMPLICIT NONE
899!~~~> Arguments
900      INTEGER,  INTENT(IN)         :: N
901      KPP_REAL, INTENT(IN)    :: Tend, Tstart
902      KPP_REAL, INTENT(INOUT) :: Y(N)
903      INTEGER,  INTENT(OUT)        :: IERR
904      INTEGER,  INTENT(IN)         :: AdjointType
905
906!~~~> Local variables
907#ifdef FULL_ALGEBRA
908      KPP_REAL    :: FJAC(N,N), E1(N,N)
909      COMPLEX(kind=dp) :: E2(N,N)   
910#else
911      KPP_REAL    :: FJAC(LU_NONZERO), E1(LU_NONZERO)
912      COMPLEX(kind=dp) :: E2(LU_NONZERO)   
913#endif               
914      KPP_REAL, DIMENSION(N) :: Z1,Z2,Z3,Z4,SCAL,DZ1,DZ2,DZ3,DZ4,G,TMP,FCN0
915      KPP_REAL  :: CONT(N,3), Tdirection,  H, T, Hacc, Hnew, Hold, Fac,  &
916                 FacGus, Theta, Err, ErrOld, NewtonRate, NewtonIncrement,  &
917                 Hratio, Qnewton, NewtonPredictedErr,NewtonIncrementOld, ThetaSD
918      INTEGER :: IP1(N),IP2(N),NewtonIter, ISING, Nconsecutive, NewIt
919      LOGICAL :: Reject, FirstStep, SkipJac, NewtonDone, SkipLU
920     
921      KPP_REAL, DIMENSION(:), POINTER :: Zstage
922           
923!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
924!~~~>  Initial setting
925!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
926      IF (AdjointType == KPP_ROOT_discrete) THEN ! Save stage solution
927        ALLOCATE(Zstage(N*3), STAT=i)
928        IF (i/=0) THEN
929          PRINT*,'Allocation of Zstage failed'
930          STOP
931        END IF
932      END IF   
933      T=Tstart
934
935      Tdirection = SIGN(ONE,Tend-Tstart)
936      H = MIN( MAX(ABS(Hmin),ABS(Hstart)) , Hmax )
937      IF (ABS(H) <= 10.d0*Roundoff) H = 1.0d-6
938      H = SIGN(H,Tdirection)
939      Hold      = H
940      Reject    = .FALSE.
941      FirstStep = .TRUE.
942      SkipJac   = .FALSE.
943      SkipLU    = .FALSE.
944      IF ((T+H*1.0001D0-Tend)*Tdirection >= ZERO) THEN
945         H = Tend-T
946      END IF
947      Nconsecutive = 0
948      CALL RK_ErrorScale(N,ITOL,AbsTol,RelTol,Y,SCAL)
949!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
950!~~~>  Time loop begins
951!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
952Tloop: DO WHILE ( (Tend-T)*Tdirection - Roundoff > ZERO )
953
954      IF ( .NOT.SkipLU ) THEN ! This time around skip the Jac update and LU
955        !~~~> Compute the Jacobian matrix
956        IF ( .NOT.SkipJac ) THEN
957          CALL JAC_CHEM(T,Y,FJAC)
958          ISTATUS(ijac) = ISTATUS(ijac) + 1
959        END IF
960        !~~~> Compute the matrices E1 and E2 and their decompositions
961        CALL RK_Decomp(N,H,FJAC,E1,IP1,E2,IP2,ISING)
962        IF (ISING /= 0) THEN
963          ISTATUS(isng) = ISTATUS(isng) + 1; Nconsecutive = Nconsecutive + 1
964          IF (Nconsecutive >= 5) THEN
965            CALL RK_ErrorMsg(-12,T,H,IERR); RETURN
966          END IF
967          H=H*0.5d0; Reject=.TRUE.; SkipJac = .TRUE.;  SkipLU = .FALSE.
968          CYCLE Tloop
969        ELSE
970          Nconsecutive = 0   
971        END IF   
972      END IF ! SkipLU
973   
974      ISTATUS(istp) = ISTATUS(istp) + 1
975      IF (ISTATUS(istp) > Max_no_steps) THEN
976        PRINT*,'Max number of time steps is ',Max_no_steps
977        CALL RK_ErrorMsg(-9,T,H,IERR); RETURN
978      END IF
979      IF (0.1D0*ABS(H) <= ABS(T)*Roundoff)  THEN
980        CALL RK_ErrorMsg(-10,T,H,IERR); RETURN
981      END IF
982     
983!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
984!~~~>  Loop for the simplified Newton iterations
985!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
986           
987      !~~~>  Starting values for Newton iteration
988      IF ( FirstStep .OR. (.NOT.StartNewton) ) THEN
989         CALL Set2zero(N,Z1)
990         CALL Set2zero(N,Z2)
991         CALL Set2zero(N,Z3)
992      ELSE
993         ! Evaluate quadratic polynomial
994         CALL RK_Interpolate('eval',N,H,Hold,Z1,Z2,Z3,CONT)
995      END IF
996     
997      !~~~>  Initializations for Newton iteration
998      NewtonDone = .FALSE.
999      Fac = 0.5d0 ! Step reduction if too many iterations
1000     
1001NewtonLoop:DO  NewtonIter = 1, NewtonMaxit 
1002 
1003            !~~~> Prepare the right-hand side
1004            CALL RK_PrepareRHS(N,T,H,Y,Z1,Z2,Z3,DZ1,DZ2,DZ3)
1005           
1006            !~~~> Solve the linear systems
1007            CALL RK_Solve( N,H,E1,IP1,E2,IP2,DZ1,DZ2,DZ3,ISING )
1008
1009
1010            NewtonIncrement = SQRT( ( RK_ErrorNorm(N,SCAL,DZ1)**2 + &
1011                                RK_ErrorNorm(N,SCAL,DZ2)**2 + &
1012                                RK_ErrorNorm(N,SCAL,DZ3)**2 )/3.0d0 )
1013           
1014            IF ( NewtonIter == 1 ) THEN
1015                Theta      = ABS(ThetaMin)
1016                NewtonRate = 2.0d0 
1017            ELSE
1018                Theta = NewtonIncrement/NewtonIncrementOld
1019                IF (Theta < 0.99d0) THEN
1020                    NewtonRate = Theta/(ONE-Theta)
1021                ELSE ! Non-convergence of Newton: Theta too large
1022                    EXIT NewtonLoop
1023                END IF
1024                IF ( NewtonIter < NewtonMaxit ) THEN 
1025                  ! Predict error at the end of Newton process
1026                  NewtonPredictedErr = NewtonIncrement &
1027                               *Theta**(NewtonMaxit-NewtonIter)/(ONE-Theta)
1028                  IF (NewtonPredictedErr >= NewtonTol) THEN
1029                    ! Non-convergence of Newton: predicted error too large
1030                    Qnewton = MIN(10.0d0,NewtonPredictedErr/NewtonTol)
1031                    Fac=0.8d0*Qnewton**(-ONE/(1+NewtonMaxit-NewtonIter))
1032                    EXIT NewtonLoop
1033                  END IF
1034                END IF
1035            END IF
1036
1037            NewtonIncrementOld = MAX(NewtonIncrement,Roundoff) 
1038            ! Update solution
1039            CALL WAXPY(N,-ONE,DZ1,1,Z1,1) ! Z1 <- Z1 - DZ1
1040            CALL WAXPY(N,-ONE,DZ2,1,Z2,1) ! Z2 <- Z2 - DZ2
1041            CALL WAXPY(N,-ONE,DZ3,1,Z3,1) ! Z3 <- Z3 - DZ3
1042           
1043            ! Check error in Newton iterations
1044            NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol)
1045            IF (NewtonDone) THEN
1046              NewIt = NewtonIter
1047              EXIT NewtonLoop
1048            END IF 
1049            IF (NewtonIter == NewtonMaxit) THEN
1050                PRINT*, 'Slow or no convergence in Newton Iteration: Max no. of', &
1051                        'Newton iterations reached'
1052            END IF
1053
1054      END DO NewtonLoop
1055
1056           
1057      IF (.NOT.NewtonDone) THEN
1058           !CALL RK_ErrorMsg(-12,T,H,IERR);
1059           H = Fac*H; Reject=.TRUE.; SkipJac = .TRUE.;  SkipLU = .FALSE.
1060           ISTATUS(irej) = ISTATUS(irej) + 1
1061           CYCLE Tloop
1062      END IF
1063
1064!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1065!~~~> SDIRK Stage
1066!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1067      IF (SdirkError) THEN
1068
1069!~~~>  Starting values for Newton iterations
1070       Z4(1:N) = Z3(1:N)
1071       
1072!~~~>   Prepare the loop-independent part of the right-hand side
1073       CALL FUN_CHEM(T,Y,DZ4)
1074       ISTATUS(Nfun) = ISTATUS(Nfun) + 1
1075       
1076!       G = H*rkBgam(0)*DZ4 + rkTheta(1)*Z1 + rkTheta(2)*Z2 + rkTheta(3)*Z3
1077       CALL Set2Zero(N, G)
1078       CALL WAXPY(N,rkBgam(0)*H, DZ4,1,G,1) 
1079       CALL WAXPY(N,rkTheta(1),Z1,1,G,1)
1080       CALL WAXPY(N,rkTheta(2),Z2,1,G,1)
1081       CALL WAXPY(N,rkTheta(3),Z3,1,G,1)
1082
1083       !~~~>  Initializations for Newton iteration
1084       NewtonDone = .FALSE.
1085       Fac = 0.5d0 ! Step reduction factor if too many iterations
1086           
1087SDNewtonLoop:DO NewtonIter = 1, NewtonMaxit
1088
1089!~~~>   Prepare the loop-dependent part of the right-hand side
1090            CALL WADD(N,Y,Z4,TMP)         ! TMP <- Y + Z4
1091            CALL FUN_CHEM(T+H,TMP,DZ4)    ! DZ4 <- Fun(Y+Z4)         
1092            ISTATUS(Nfun) = ISTATUS(Nfun) + 1
1093!            DZ4(1:N) = (G(1:N)-Z4(1:N))*(rkGamma/H) + DZ4(1:N)
1094            CALL WAXPY (N, -ONE*rkGamma/H, Z4, 1, DZ4, 1)
1095            CALL WAXPY (N, rkGamma/H, G,1, DZ4,1)
1096
1097!~~~>   Solve the linear system
1098#ifdef FULL_ALGEBRA 
1099            CALL DGETRS( 'N', N, 1, E1, N, IP1, DZ4, N, ISING )
1100#else
1101            CALL KppSolve(E1, DZ4)
1102#endif
1103           
1104!~~~>   Check convergence of Newton iterations
1105            NewtonIncrement = RK_ErrorNorm(N,SCAL,DZ4)
1106            IF ( NewtonIter == 1 ) THEN
1107                ThetaSD      = ABS(ThetaMin)
1108                NewtonRate = 2.0d0 
1109            ELSE
1110                ThetaSD = NewtonIncrement/NewtonIncrementOld
1111                IF (ThetaSD < 0.99d0) THEN
1112                    NewtonRate = ThetaSD/(ONE-ThetaSD)
1113                    ! Predict error at the end of Newton process
1114                    NewtonPredictedErr = NewtonIncrement &
1115                               *ThetaSD**(NewtonMaxit-NewtonIter)/(ONE-ThetaSD)
1116                    IF (NewtonPredictedErr >= NewtonTol) THEN
1117                      ! Non-convergence of Newton: predicted error too large
1118                      !print*,'Error too large: ', NewtonPredictedErr
1119                      Qnewton = MIN(10.0d0,NewtonPredictedErr/NewtonTol)
1120                      Fac = 0.8d0*Qnewton**(-ONE/(1+NewtonMaxit-NewtonIter))
1121                      EXIT SDNewtonLoop
1122                    END IF
1123                ELSE ! Non-convergence of Newton: ThetaSD too large
1124                    !print*,'Theta too large: ',ThetaSD
1125                    EXIT SDNewtonLoop
1126                END IF
1127            END IF
1128            NewtonIncrementOld = NewtonIncrement
1129            ! Update solution: Z4 <-- Z4 + DZ4
1130            CALL WAXPY(N,ONE,DZ4,1,Z4,1) 
1131           
1132            ! Check error in Newton iterations
1133            NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol)
1134            IF (NewtonDone) EXIT SDNewtonLoop
1135           
1136            END DO SDNewtonLoop
1137           
1138            IF (.NOT.NewtonDone) THEN
1139                H = Fac*H; Reject=.TRUE.; SkipJac = .TRUE.;  SkipLU = .FALSE.
1140                CYCLE Tloop
1141            END IF
1142      END IF
1143!~~~>  End of implified SDIRK Newton iterations
1144
1145           
1146!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1147!~~~> Error estimation
1148!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1149      IF (SdirkError) THEN
1150!         DZ4(1:N) =  rkD(1)*Z1 + rkD(2)*Z2 + rkD(3)*Z3 - Z4   
1151         CALL Set2Zero(N, DZ4)
1152         IF (rkD(1) /= ZERO) CALL WAXPY(N, rkD(1), Z1, 1, DZ4, 1)
1153         IF (rkD(2) /= ZERO) CALL WAXPY(N, rkD(2), Z2, 1, DZ4, 1)
1154         IF (rkD(3) /= ZERO) CALL WAXPY(N, rkD(3), Z3, 1, DZ4, 1)
1155         CALL WAXPY(N, -ONE, Z4, 1, DZ4, 1)
1156         Err = RK_ErrorNorm(N,SCAL,DZ4)   
1157      ELSE
1158         CALL  RK_ErrorEstimate(N,H,Y,T, &
1159               E1,IP1,Z1,Z2,Z3,SCAL,Err,FirstStep,Reject)
1160      END IF
1161!~~~> Computation of new step size Hnew
1162      Fac  = Err**(-ONE/rkELO)*          &
1163             MIN(FacSafe,(ONE+2*NewtonMaxit)/(NewtonIter+2*NewtonMaxit))
1164      Fac  = MIN(FacMax,MAX(FacMin,Fac))
1165      Hnew = Fac*H
1166
1167!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1168!~~~> Accept/reject step
1169!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1170accept:IF (Err < ONE) THEN !~~~> STEP IS ACCEPTED
1171         IF (AdjointType == KPP_ROOT_discrete) THEN ! Save stage solution
1172            ! CALL WCOPY(N,Z1,1,Zstage(1),1)
1173            ! CALL WCOPY(N,Z2,1,Zstage(N+1),1)
1174            ! CALL WCOPY(N,Z3,1,Zstage(2*N+1),1)
1175            Zstage(1:N)       = Z1(1:N) 
1176            Zstage(N+1:2*N)   = Z2(1:N)
1177            Zstage(2*N+1:3*N) = Z3(1:N)
1178            ! Push old Y - Y at the beginning of the stage
1179            CALL rk_DPush(T, H, Y, Zstage, NewIt, E1, E2)
1180         END IF
1181
1182         FirstStep=.FALSE.
1183         ISTATUS(iacc) = ISTATUS(iacc) + 1
1184         IF (Gustafsson) THEN
1185            !~~~> Predictive controller of Gustafsson
1186            IF (ISTATUS(iacc) > 1) THEN
1187               FacGus=FacSafe*(H/Hacc)*(Err**2/ErrOld)**(-0.25d0)
1188               FacGus=MIN(FacMax,MAX(FacMin,FacGus))
1189               Fac=MIN(Fac,FacGus)
1190               Hnew = Fac*H
1191            END IF
1192            Hacc=H
1193            ErrOld=MAX(1.0d-2,Err)
1194         END IF
1195         Hold = H
1196         T = T+H 
1197         ! Update solution: Y <- Y + sum (d_i Z_i)
1198         IF (rkD(1) /= ZERO) CALL WAXPY(N,rkD(1),Z1,1,Y,1)
1199         IF (rkD(2) /= ZERO) CALL WAXPY(N,rkD(2),Z2,1,Y,1)
1200         IF (rkD(3) /= ZERO) CALL WAXPY(N,rkD(3),Z3,1,Y,1)
1201         ! Construct the solution quadratic interpolant Q(c_i) = Z_i, i=1:3
1202         IF (StartNewton) CALL RK_Interpolate('make',N,H,Hold,Z1,Z2,Z3,CONT)
1203         CALL RK_ErrorScale(N,ITOL,AbsTol,RelTol,Y,SCAL)
1204         RSTATUS(itexit) = T
1205         RSTATUS(ihnew)  = Hnew
1206         RSTATUS(ihacc)  = H
1207         Hnew = Tdirection*MIN( MAX(ABS(Hnew),Hmin) , Hmax )
1208         IF (Reject) Hnew = Tdirection*MIN(ABS(Hnew),ABS(H))
1209         Reject = .FALSE.
1210         IF ((T+Hnew/Qmin-Tend)*Tdirection >=  ZERO) THEN
1211            H = Tend-T
1212         ELSE
1213            Hratio=Hnew/H
1214            ! Reuse the LU decomposition
1215            SkipLU = (Theta<=ThetaMin) .AND. (Hratio>=Qmin) .AND. (Hratio<=Qmax)
1216            SkipLU = .false.
1217            IF (.NOT.SkipLU) H=Hnew
1218         END IF
1219         ! If convergence is fast enough, do not update Jacobian
1220         ! SkipJac = (Theta <= ThetaMin)
1221         SkipJac = .FALSE.
1222
1223      ELSE accept !~~~> Step is rejected
1224         IF (FirstStep .OR. Reject) THEN
1225             H = FacRej*H
1226         ELSE
1227             H = Hnew
1228         END IF
1229         Reject   = .TRUE.
1230         SkipJac  = .TRUE.
1231         SkipLU   = .FALSE. 
1232         IF (ISTATUS(iacc) >= 1) ISTATUS(irej) = ISTATUS(irej) + 1
1233      END IF accept
1234     
1235    END DO Tloop
1236
1237!~~~> Save last state: only needed for continuous adjoint
1238   IF ( (AdjointType == KPP_ROOT_continuous) .OR. &
1239       (AdjointType == KPP_ROOT_simple_continuous) ) THEN
1240       CALL Fun_Chem(T,Y,Fcn0)
1241       CALL rk_CPush( T, H, Y, Fcn0, DZ3)
1242!~~~> Deallocate stage buffer: only needed for discrete adjoint
1243   ELSEIF (AdjointType == KPP_ROOT_discrete) THEN
1244        DEALLOCATE(Zstage, STAT=i)
1245        IF (i/=0) THEN
1246          PRINT*,'Deallocation of Zstage failed'
1247          STOP
1248        END IF
1249   END IF   
1250   
1251    ! Successful exit
1252    IERR = 1 
1253
1254 END SUBROUTINE RK_FwdIntegrator
1255
1256
1257!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1258 SUBROUTINE rk_DadjInt( NADJ,Lambda,Tstart,Tend,T,IERR )
1259!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1260
1261      IMPLICIT NONE
1262!~~~> Arguments
1263!~~~> Input: the initial condition at Tstart; Output: the solution at T   
1264      INTEGER, INTENT(IN)     :: NADJ
1265!~~~> First order adjoint   
1266      KPP_REAL, INTENT(INOUT) :: Lambda(N,NADJ)
1267      KPP_REAL, INTENT(IN)    :: Tstart, Tend
1268      KPP_REAL, INTENT(INOUT) :: T
1269      INTEGER,  INTENT(OUT)        :: IERR
1270
1271!~~~> Local variables
1272#ifdef FULL_ALGEBRA
1273      KPP_REAL    :: Jac0(N,N), Jac1(N,N),Jac2(N,N),Jac3(N,N), E1(N,N)
1274      COMPLEX(kind=dp) :: E2(N,N)   
1275      KPP_REAL    :: Jbig(3*N,3*N), X(3*N)
1276      INTEGER          :: IPbig(3*N)
1277#else
1278      KPP_REAL, DIMENSION(LU_NONZERO) :: Jac0, Jac1, Jac2, Jac3, E1
1279      COMPLEX(kind=dp), DIMENSION(LU_NONZERO) :: E2   
1280      ! Next two lines for sparse big algebra:
1281      ! KPP_REAL :: Jbig(3,3,LU_NONZERO), X(3,N)
1282      ! INTEGER :: IPbig(3,N)
1283      KPP_REAL :: Jbig(3*N,3*N), X(3*N)
1284      INTEGER :: IPbig(3*N)
1285#endif               
1286      KPP_REAL, DIMENSION(N,NADJ) :: U1,U2,U3
1287      KPP_REAL, DIMENSION(N) :: SCAL,DU1,DU2,DU3
1288      KPP_REAL :: Y(N), Zstage(3*N)
1289      KPP_REAL  :: H, NewtonRate, NewtonIncrement, NewtonIncrementOld
1290      INTEGER :: IP1(N),IP2(N),NewtonIter, ISING, iadj, NewIt
1291      LOGICAL :: Reject, NewtonDone, NewtonConverge
1292      KPP_REAL :: T1, Theta
1293      KPP_REAL, DIMENSION(N) :: TMP, G1, G2, G3
1294      INTEGER :: info, i, j
1295           
1296!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1297!~~~>  Initial setting
1298!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1299      T1 = Tend     
1300!      Tdirection = SIGN(ONE,Tend-Tstart)
1301      NewtonConverge = .TRUE.
1302      Reject = .FALSE.
1303     
1304!~~~> Time loop begins below
1305TimeLoop:DO WHILE ( stack_ptr > 0 )
1306
1307   IF (.not.Reject) THEN
1308   
1309     !~~~>  Recover checkpoints for stage values and vectors
1310     CALL rk_DPop(T, H, Y, Zstage, NewIt, E1, E2)
1311
1312     !~~~>  Compute LU decomposition
1313     IF (.NOT.SaveLU) THEN
1314        !~~~> Compute the Jacobian matrix
1315        CALL JAC_CHEM(T,Y,Jac0)
1316        ISTATUS(ijac) = ISTATUS(ijac) + 1
1317        !~~~> Compute the matrices E1 and E2 and their decompositions
1318        CALL RK_Decomp(N,H,Jac0,E1,IP1,E2,IP2,ISING)
1319     END IF
1320   
1321     !~~~>   Jacobian values at stage vectors
1322     CALL WADD(N,Y,Zstage(1),TMP)       ! TMP  <- Y + Z1
1323     CALL JAC_CHEM(T+rkC(1)*H,TMP,Jac1) ! Jac1 <- Jac(Y+Z1) 
1324     CALL WADD(N,Y,Zstage(1+N),TMP)     ! TMP  <- Y + Z2
1325     CALL JAC_CHEM(T+rkC(2)*H,TMP,Jac2) ! Jac2 <- Jac(Y+Z2) 
1326     CALL WADD(N,Y,Zstage(1+2*N),TMP)   ! TMP  <- Y + Z3
1327     CALL JAC_CHEM(T+rkC(3)*H,TMP,Jac3) ! Jac3 <- Jac(Y+Z3)
1328     
1329   END IF ! .not.Reject
1330
1331111 CONTINUE
1332
1333   IF ( (AdjointSolve == Solve_adaptive .and. .not.NewtonConverge) &
1334         .or. (AdjointSolve == Solve_direct) ) THEN
1335#ifdef FULL_ALGEBRA 
1336       DO j=1,N
1337         DO i=1,N
1338          Jbig(i,j)         = -H*rkA(1,1)*Jac1(i,j)
1339          Jbig(N+i,j)       = -H*rkA(2,1)*Jac1(i,j)
1340          Jbig(2*N+i,j)     = -H*rkA(3,1)*Jac1(i,j)
1341          Jbig(i,N+j)       = -H*rkA(1,2)*Jac2(i,j)
1342          Jbig(N+i,N+j)     = -H*rkA(2,2)*Jac2(i,j)
1343          Jbig(2*N+i,N+j)   = -H*rkA(3,2)*Jac2(i,j)
1344          Jbig(i,2*N+j)     = -H*rkA(1,3)*Jac3(i,j)
1345          Jbig(N+i,2*N+j)   = -H*rkA(2,3)*Jac3(i,j)
1346          Jbig(2*N+i,2*N+j) = -H*rkA(3,3)*Jac3(i,j)
1347         END DO 
1348       END DO
1349       DO i=1,3*N
1350          Jbig(i,i) = Jbig(i,i) + ONE
1351       END DO
1352       CALL DGETRF(3*N,3*N,Jbig,3*N,IPbig,info) 
1353       ! CALL WGEFA(3*N,Jbig,IPbig,info)
1354       IF (info /= 0) THEN
1355         PRINT*,'Big guy is singular'; STOP
1356       END IF 
1357#else   
1358!  Commented lines for sparse big algebra:   
1359!      !~~~>  Construct the big Jacobian
1360!      DO j=1,LU_NONZERO
1361!        DO i=1,3
1362!          Jbig(i,1,j) = -H*rkA(i,1)*Jac1(j)
1363!          Jbig(i,2,j) = -H*rkA(i,2)*Jac2(j)
1364!          Jbig(i,3,j) = -H*rkA(i,3)*Jac3(j)
1365!        END DO
1366!      END DO
1367!      DO j=1,NVAR
1368!        DO i=1,3
1369!          Jbig(i,i,LU_DIAG(j)) = ONE + Jbig(i,i,LU_DIAG(j))
1370!        END DO
1371!      END DO
1372!      !~~~>  Solve the big system
1373!      CALL KppDecompBig( Jbig, IPbig, info )
1374!  Use full big algebra:   
1375      Jbig(1:3*N,1:3*N) = 0.0d0
1376      DO i=1,LU_NONZERO
1377          Jbig(LU_irow(i),LU_icol(i))         = -H*rkA(1,1)*Jac1(i)
1378          Jbig(LU_irow(i),N+LU_icol(i))       = -H*rkA(1,2)*Jac2(i)
1379          Jbig(LU_irow(i),2*N+LU_icol(i))     = -H*rkA(1,3)*Jac3(i)
1380          Jbig(N+LU_irow(i),LU_icol(i))       = -H*rkA(2,1)*Jac1(i)
1381          Jbig(N+LU_irow(i),N+LU_icol(i))     = -H*rkA(2,2)*Jac2(i)
1382          Jbig(N+LU_irow(i),2*N+LU_icol(i))   = -H*rkA(2,3)*Jac3(i)
1383          Jbig(2*N+LU_irow(i),LU_icol(i))     = -H*rkA(3,1)*Jac1(i)
1384          Jbig(2*N+LU_irow(i),N+LU_icol(i))   = -H*rkA(3,2)*Jac2(i)
1385          Jbig(2*N+LU_irow(i),2*N+LU_icol(i)) = -H*rkA(3,3)*Jac3(i)
1386      END DO
1387      DO i=1, 3*N
1388         Jbig(i,i) = ONE + Jbig(i,i)
1389      END DO
1390      ! CALL DGETRF(3*N,3*N,Jbig,3*N,IPbig,info)
1391      CALL WGEFA(3*N,Jbig,IPbig,info) 
1392      IF (info /= 0) THEN
1393         PRINT*,'Big guy is singular'; STOP
1394      END IF
1395#endif       
1396    END IF
1397     
1398!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1399!~~~>  Loop for the simplified Newton iterations
1400!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1401Adj:DO iadj = 1, NADJ     
1402      !~~~>  Starting values for Newton iteration
1403      ! CALL WCOPY(N,Lambda(1,iadj),1,U1(1,iadj),1)
1404      ! CALL WCOPY(N,Lambda(1,iadj),1,U2(1,iadj),1)
1405      ! CALL WCOPY(N,Lambda(1,iadj),1,U3(1,iadj),1)
1406      CALL Set2Zero(N,U1(1,iadj))
1407      CALL Set2Zero(N,U2(1,iadj))
1408      CALL Set2Zero(N,U3(1,iadj))
1409     
1410      !~~~>  Initializations for Newton iteration
1411      NewtonDone = .FALSE.
1412      !~~~>    Right Hand Side - part G for all Newton iterations
1413      CALL RK_PrepareRHS_G(N,H,Jac1,Jac2,Jac3,Lambda(1,iadj),         &
1414                        G1, G2, G3)
1415                       
1416      IF ( (AdjointSolve == Solve_adaptive .and. NewtonConverge)      &
1417            .or. (AdjointSolve == Solve_fixed) ) THEN
1418
1419NewtonLoopAdj:DO  NewtonIter = 1, NewtonMaxit 
1420
1421            !~~~> Prepare the right-hand side
1422            CALL RK_PrepareRHS_Adj(N,H,Jac1,Jac2,Jac3,Lambda(1,iadj), &
1423                        U1(1,iadj),U2(1,iadj),U3(1,iadj),             &
1424                        G1, G2, G3, DU1,DU2,DU3)
1425
1426            !~~~> Solve the linear systems
1427            CALL RK_SolveTR( N,H,E1,IP1,E2,IP2,DU1,DU2,DU3,ISING )
1428
1429!~~~> The following code performs an adaptive number of Newton
1430!     iterations for solving adjoint system
1431            IF (AdjointSolve == Solve_adaptive) THEN
1432
1433            CALL RK_ErrorScale(N,ITOL,                      &
1434                 AbsTol_adj(1:N,iadj),RelTol_adj(1:N,iadj), &
1435                 Lambda(1:N,iadj),SCAL)
1436
1437            ! SCAL(1:N) = 1.0d0
1438            NewtonIncrement = SQRT( ( RK_ErrorNorm(N,SCAL,DU1)**2 +   &
1439                                RK_ErrorNorm(N,SCAL,DU2)**2 +         &
1440                                RK_ErrorNorm(N,SCAL,DU3)**2 )/3.0d0 )
1441           
1442       
1443            IF ( NewtonIter == 1 ) THEN
1444                Theta      = ABS(ThetaMin)
1445                NewtonRate = 2.0d0 
1446            ELSE
1447                Theta = NewtonIncrement/NewtonIncrementOld
1448                IF (Theta < 0.99d0) THEN
1449                    NewtonRate = Theta/(ONE-Theta)
1450                ELSE ! Non-convergence of Newton: Theta too large
1451                    Reject = .TRUE.
1452                    NewtonDone = .FALSE.
1453                    EXIT NewtonLoopAdj
1454                END IF
1455
1456            END IF
1457         
1458            NewtonIncrementOld = MAX(NewtonIncrement,Roundoff) 
1459
1460            END IF ! (AdjointSolve == Solve_adaptive)
1461           
1462            ! Update solution
1463            CALL WAXPY(N,-ONE,DU1,1,U1(1,iadj),1) ! U1 <- U1 - DU1
1464            CALL WAXPY(N,-ONE,DU2,1,U2(1,iadj),1) ! U2 <- U2 - DU2
1465            CALL WAXPY(N,-ONE,DU3,1,U3(1,iadj),1) ! U3 <- U3 - DU3
1466
1467            IF (AdjointSolve == Solve_adaptive) THEN
1468            ! When performing an adaptive number of iterations
1469            !       check the error in Newton iterations
1470               NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol)
1471               IF ((NewtonDone).and.(NewtonIter>NewIt+1)) EXIT NewtonLoopAdj
1472            ELSE IF (AdjointSolve == Solve_fixed) THEN
1473               IF (NewtonIter>MAX(NewIt+1,6)) EXIT NewtonLoopAdj
1474            END IF   
1475           
1476      END DO NewtonLoopAdj
1477     
1478      IF ((AdjointSolve==Solve_adaptive).AND.(.NOT.NewtonDone)) THEN
1479         ! print*,'Newton iterations do not converge, switching to full system.'
1480         NewtonConverge = .FALSE.
1481         Reject = .TRUE.
1482         GOTO 111
1483      END IF
1484     
1485      ! Update adjoint solution: Y_adj <- Y_adj + sum (U_i)
1486      CALL WAXPY(N,ONE,U1(1,iadj),1,Lambda(1,iadj),1)
1487      CALL WAXPY(N,ONE,U2(1,iadj),1,Lambda(1,iadj),1)
1488      CALL WAXPY(N,ONE,U3(1,iadj),1,Lambda(1,iadj),1)
1489
1490      ELSE ! NewtonConverge = .false.
1491
1492#ifdef FULL_ALGEBRA 
1493        X(1:N)       = -G1(1:N)
1494        X(N+1:2*N)   = -G2(1:N)
1495        X(2*N+1:3*N) = -G3(1:N)
1496        CALL DGETRS('T',3*N,1,Jbig,3*N,IPbig,X,3*N,0) 
1497        ! CALL WGESL('T',3*N,Jbig,IPbig,X)
1498        Lambda(1:N,iadj) = Lambda(1:N,iadj)+X(1:N)+X(N+1:2*N)+X(2*N+1:3*N)
1499#else       
1500!   Commented lines for sparse big algebra:
1501!       X(1,1:N) = -G1(1:N)
1502!       X(2,1:N) = -G2(1:N)
1503!       X(3,1:N) = -G3(1:N)
1504!       CALL KppSolveBigTR( Jbig, IPbig, X )               
1505!       Lambda(1:N,iadj) = Lambda(1:N,iadj)+X(1,1:N)+X(2,1:N)+X(3,1:N)
1506!   Use fill big algebra:
1507        X(1:N)       = -G1(1:N)
1508        X(N+1:2*N)   = -G2(1:N)
1509        X(2*N+1:3*N) = -G3(1:N)
1510        ! CALL DGETRS('T',3*N,1,Jbig,3*N,IPbig,X,3*N,0)
1511        CALL WGESL('T',3*N,Jbig,IPbig,X)
1512        Lambda(1:N,iadj) = Lambda(1:N,iadj)+X(1:N)+X(N+1:2*N)+X(2*N+1:3*N)
1513#endif       
1514        IF ((AdjointSolve==Solve_adaptive).AND.(iadj>=NADJ)) THEN
1515          NewtonConverge = .TRUE.
1516          Reject = .FALSE.
1517        END IF
1518     
1519     END IF ! NewtonConverge
1520 
1521   END DO Adj
1522   
1523   T1 = T1-H
1524     
1525   END DO TimeLoop
1526   
1527    ! Successful exit
1528    IERR = 1 
1529
1530 END SUBROUTINE RK_DadjInt
1531
1532
1533!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1534 SUBROUTINE rk_CadjInt (                        &
1535        NADJ, Y,                                &
1536        Tstart, Tend, T, IERR)
1537!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1538      IMPLICIT NONE
1539!~~~> Arguments
1540!~~~> Input: the initial condition at Tstart; Output: the solution at T   
1541      INTEGER, INTENT(IN)     :: NADJ
1542!~~~> First order adjoint   
1543      KPP_REAL, INTENT(INOUT) :: Y(N,NADJ)
1544      KPP_REAL, INTENT(IN)    :: Tstart, Tend
1545      KPP_REAL, INTENT(INOUT) :: T
1546      INTEGER,  INTENT(OUT)        :: IERR
1547
1548  END SUBROUTINE rk_CadjInt
1549
1550!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1551 SUBROUTINE rk_SimpleCadjInt (                  &
1552        NADJ, Y,                                &
1553        Tstart, Tend, T,                        &
1554        IERR )
1555!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1556      IMPLICIT NONE
1557!~~~> Arguments
1558!~~~> Input: the initial condition at Tstart; Output: the solution at T   
1559      INTEGER, INTENT(IN)     :: NADJ
1560!~~~> First order adjoint   
1561      KPP_REAL, INTENT(INOUT) :: Y(N,NADJ)
1562      KPP_REAL, INTENT(IN)    :: Tstart, Tend
1563      KPP_REAL, INTENT(INOUT) :: T
1564      INTEGER,  INTENT(OUT)        :: IERR
1565
1566  END SUBROUTINE rk_SimpleCadjInt
1567
1568
1569!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1570 SUBROUTINE RK_ErrorMsg(Code,T,H,IERR)
1571!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1572!    Handles all error messages
1573!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1574
1575   IMPLICIT NONE
1576   KPP_REAL, INTENT(IN) :: T, H
1577   INTEGER, INTENT(IN)  :: Code
1578   INTEGER, INTENT(OUT) :: IERR
1579
1580   IERR = Code
1581   PRINT * , &
1582     'Forced exit from RungeKutta due to the following error:'
1583
1584
1585   SELECT CASE (Code)
1586    CASE (-1)
1587      PRINT * , '--> Improper value for maximal no of steps'
1588    CASE (-2)
1589      PRINT * , '--> Improper value for maximal no of Newton iterations'
1590    CASE (-3)
1591      PRINT * , '--> Hmin/Hmax/Hstart must be positive'
1592    CASE (-4)
1593      PRINT * , '--> Improper values for FacMin/FacMax/FacSafe/FacRej'
1594    CASE (-5)
1595      PRINT * , '--> Improper value for ThetaMin'
1596    CASE (-6)
1597      PRINT * , '--> Newton stopping tolerance too small'
1598    CASE (-7)
1599      PRINT * , '--> Improper values for Qmin, Qmax'
1600    CASE (-8)
1601      PRINT * , '--> Tolerances are too small'
1602    CASE (-9)
1603      PRINT * , '--> No of steps exceeds maximum bound'
1604    CASE (-10)
1605      PRINT * , '--> Step size too small: T + 10*H = T', &
1606            ' or H < Roundoff'
1607    CASE (-11)
1608      PRINT * , '--> Matrix is repeatedly singular'
1609    CASE (-12)
1610      PRINT * , '--> Non-convergence of Newton iterations'
1611    CASE (-13)
1612      PRINT * , '--> Requested RK method not implemented'
1613    CASE DEFAULT
1614      PRINT *, 'Unknown Error code: ', Code
1615   END SELECT
1616
1617   WRITE(6,FMT="(5X,'T=',E12.5,'  H=',E12.5)") T, H 
1618
1619 END SUBROUTINE RK_ErrorMsg
1620
1621
1622!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1623 SUBROUTINE RK_ErrorScale(N,ITOL,AbsTol,RelTol,Y,SCAL)
1624!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1625!    Handles all error messages
1626!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1627   IMPLICIT NONE
1628   INTEGER, INTENT(IN)  :: N, ITOL
1629   KPP_REAL, INTENT(IN) :: AbsTol(*), RelTol(*), Y(N)
1630   KPP_REAL, INTENT(OUT) :: SCAL(N)
1631   INTEGER :: i
1632   
1633   IF (ITOL==0) THEN
1634       DO i=1,N
1635          SCAL(i)= ONE/(AbsTol(1)+RelTol(1)*ABS(Y(i)))
1636       END DO
1637   ELSE
1638       DO i=1,N
1639          SCAL(i)=ONE/(AbsTol(i)+RelTol(i)*ABS(Y(i)))
1640       END DO
1641   END IF
1642     
1643 END SUBROUTINE RK_ErrorScale
1644
1645
1646!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1647!  SUBROUTINE RK_Transform(N,Tr,Z1,Z2,Z3,W1,W2,W3)
1648!!~~~>                 W <-- Tr x Z
1649!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1650!      IMPLICIT NONE
1651!      INTEGER :: N, i
1652!      KPP_REAL :: Tr(3,3),Z1(N),Z2(N),Z3(N),W1(N),W2(N),W3(N)
1653!      KPP_REAL :: x1, x2, x3
1654!      DO i=1,N
1655!          x1 = Z1(i); x2 = Z2(i); x3 = Z3(i)
1656!          W1(i) = Tr(1,1)*x1 + Tr(1,2)*x2 + Tr(1,3)*x3
1657!          W2(i) = Tr(2,1)*x1 + Tr(2,2)*x2 + Tr(2,3)*x3
1658!          W3(i) = Tr(3,1)*x1 + Tr(3,2)*x2 + Tr(3,3)*x3
1659!      END DO
1660!  END SUBROUTINE RK_Transform
1661 
1662!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1663  SUBROUTINE RK_Interpolate(action,N,H,Hold,Z1,Z2,Z3,CONT)
1664!~~~>   Constructs or evaluates a quadratic polynomial
1665!         that interpolates the Z solution at current step
1666!         and provides starting values for the next step
1667!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1668      INTEGER :: N, i
1669      KPP_REAL :: H,Hold,Z1(N),Z2(N),Z3(N),CONT(N,3)
1670      KPP_REAL :: r, x1, x2, x3, den
1671      CHARACTER(LEN=4) :: action
1672       
1673      SELECT CASE (action) 
1674      CASE ('make')
1675         ! Construct the solution quadratic interpolant Q(c_i) = Z_i, i=1:3
1676         den = (rkC(3)-rkC(2))*(rkC(2)-rkC(1))*(rkC(1)-rkC(3))
1677         DO i=1,N
1678             CONT(i,1)=(-rkC(3)**2*rkC(2)*Z1(i)+Z3(i)*rkC(2)*rkC(1)**2 &
1679                        +rkC(2)**2*rkC(3)*Z1(i)-rkC(2)**2*rkC(1)*Z3(i) &
1680                        +rkC(3)**2*rkC(1)*Z2(i)-Z2(i)*rkC(3)*rkC(1)**2)&
1681                        /den-Z3(i)
1682             CONT(i,2)= -( rkC(1)**2*(Z3(i)-Z2(i)) + rkC(2)**2*(Z1(i)  &
1683                          -Z3(i)) +rkC(3)**2*(Z2(i)-Z1(i)) )/den
1684             CONT(i,3)= ( rkC(1)*(Z3(i)-Z2(i)) + rkC(2)*(Z1(i)-Z3(i))  &
1685                           +rkC(3)*(Z2(i)-Z1(i)) )/den
1686         END DO
1687      CASE ('eval')
1688          ! Evaluate quadratic polynomial
1689          r = H/Hold
1690         x1 = ONE + rkC(1)*r
1691         x2 = ONE + rkC(2)*r
1692         x3 = ONE + rkC(3)*r
1693         DO i=1,N
1694            Z1(i) = CONT(i,1)+x1*(CONT(i,2)+x1*CONT(i,3))
1695            Z2(i) = CONT(i,1)+x2*(CONT(i,2)+x2*CONT(i,3))
1696            Z3(i) = CONT(i,1)+x3*(CONT(i,2)+x3*CONT(i,3))
1697         END DO
1698       END SELECT   
1699  END SUBROUTINE RK_Interpolate
1700
1701
1702!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1703   SUBROUTINE RK_PrepareRHS(N,T,H,Y,Z1,Z2,Z3,R1,R2,R3)
1704!~~~> Prepare the right-hand side for Newton iterations
1705!     R = Z - hA x F
1706!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1707      IMPLICIT NONE
1708     
1709      INTEGER :: N
1710      KPP_REAL :: T,H
1711      KPP_REAL, DIMENSION(N) :: Y,Z1,Z2,Z3,F,R1,R2,R3,TMP
1712
1713      CALL WCOPY(N,Z1,1,R1,1) ! R1 <- Z1
1714      CALL WCOPY(N,Z2,1,R2,1) ! R2 <- Z2
1715      CALL WCOPY(N,Z3,1,R3,1) ! R3 <- Z3
1716
1717      CALL WADD(N,Y,Z1,TMP)              ! TMP <- Y + Z1
1718      CALL FUN_CHEM(T+rkC(1)*H,TMP,F)    ! F1 <- Fun(Y+Z1)         
1719      CALL WAXPY(N,-H*rkA(1,1),F,1,R1,1) ! R1 <- R1 - h*A_11*F1
1720      CALL WAXPY(N,-H*rkA(2,1),F,1,R2,1) ! R2 <- R2 - h*A_21*F1
1721      CALL WAXPY(N,-H*rkA(3,1),F,1,R3,1) ! R3 <- R3 - h*A_31*F1
1722
1723      CALL WADD(N,Y,Z2,TMP)              ! TMP <- Y + Z2
1724      CALL FUN_CHEM(T+rkC(2)*H,TMP,F)    ! F2 <- Fun(Y+Z2)       
1725      CALL WAXPY(N,-H*rkA(1,2),F,1,R1,1) ! R1 <- R1 - h*A_12*F2
1726      CALL WAXPY(N,-H*rkA(2,2),F,1,R2,1) ! R2 <- R2 - h*A_22*F2
1727      CALL WAXPY(N,-H*rkA(3,2),F,1,R3,1) ! R3 <- R3 - h*A_32*F2
1728
1729      CALL WADD(N,Y,Z3,TMP)              ! TMP <- Y + Z3
1730      CALL FUN_CHEM(T+rkC(3)*H,TMP,F)    ! F3 <- Fun(Y+Z3)     
1731      CALL WAXPY(N,-H*rkA(1,3),F,1,R1,1) ! R1 <- R1 - h*A_13*F3
1732      CALL WAXPY(N,-H*rkA(2,3),F,1,R2,1) ! R2 <- R2 - h*A_23*F3
1733      CALL WAXPY(N,-H*rkA(3,3),F,1,R3,1) ! R3 <- R3 - h*A_33*F3
1734           
1735  END SUBROUTINE RK_PrepareRHS
1736 
1737!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1738   SUBROUTINE RK_PrepareRHS_Adj(N,H,Jac1,Jac2,Jac3,Lambda, &
1739                      U1,U2,U3,G1,G2,G3,R1,R2,R3)
1740!~~~> Prepare the right-hand side for Newton iterations
1741!     R = Z_adj - hA x Jac*Z_adj - h J^t b lambda
1742!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1743      IMPLICIT NONE
1744     
1745      INTEGER, INTENT(IN) :: N
1746      KPP_REAL, INTENT(IN) :: H
1747      KPP_REAL, DIMENSION(N), INTENT(IN) :: U1,U2,U3,Lambda,G1,G2,G3
1748      KPP_REAL, DIMENSION(N), INTENT(OUT) :: R1,R2,R3
1749#ifdef FULL_ALGEBRA     
1750      KPP_REAL, DIMENSION(N,N), INTENT(IN) :: Jac1, Jac2, Jac3
1751#else     
1752      KPP_REAL, DIMENSION(LU_NONZERO),INTENT(IN) :: Jac1, Jac2, Jac3
1753#endif
1754      KPP_REAL, DIMENSION(N) ::  F,TMP
1755
1756
1757      CALL WCOPY(N,G1,1,R1,1) ! R1 <- G1
1758      CALL WCOPY(N,G2,1,R2,1) ! R2 <- G2
1759      CALL WCOPY(N,G3,1,R3,1) ! R3 <- G3
1760
1761      CALL SET2ZERO(N,F)
1762      CALL WAXPY(N,-H*rkA(1,1),U1,1,F,1) ! F1 <- -h*A_11*U1
1763      CALL WAXPY(N,-H*rkA(2,1),U2,1,F,1) ! F1 <- F1 - h*A_21*U2
1764      CALL WAXPY(N,-H*rkA(3,1),U3,1,F,1) ! F1 <- F1 - h*A_31*U3
1765#ifdef FULL_ALGEBRA 
1766      TMP = MATMUL(TRANSPOSE(Jac1),F)   
1767#else     
1768      CALL JacTR_SP_Vec ( Jac1, F, TMP )   ! R1 <- -Jac(Y+Z1)^t*h*sum(A_j1*U_j) 
1769#endif     
1770      CALL WAXPY(N,ONE,U1,1,TMP,1) ! R1 <- U1 -Jac(Y+Z1)^t*h*sum(A_j1*U_j)
1771      CALL WAXPY(N,ONE,TMP,1,R1,1) ! R1 <- U1 -Jac(Y+Z1)^t*h*sum(A_j1*U_j)
1772
1773      CALL SET2ZERO(N,F)
1774      CALL WAXPY(N,-H*rkA(1,2),U1,1,F,1) ! F2 <- -h*A_11*U1
1775      CALL WAXPY(N,-H*rkA(2,2),U2,1,F,1) ! F2 <- F2 - h*A_21*U2
1776      CALL WAXPY(N,-H*rkA(3,2),U3,1,F,1) ! F2 <- F2 - h*A_31*U3
1777#ifdef FULL_ALGEBRA 
1778      TMP = MATMUL(TRANSPOSE(Jac2),F)   
1779#else     
1780      CALL JacTR_SP_Vec ( Jac2, F, TMP )   ! R2 <- -Jac(Y+Z2)^t*h*sum(A_j2*U_j) 
1781#endif     
1782      CALL WAXPY(N,ONE,U2,1,TMP,1) ! R2 <- U2 -Jac(Y+Z2)^t*h*sum(A_j2*U_j)
1783      CALL WAXPY(N,ONE,TMP,1,R2,1) ! R2 <- U2 -Jac(Y+Z2)^t*h*sum(A_j2*U_j)
1784
1785      CALL SET2ZERO(N,F)
1786      CALL WAXPY(N,-H*rkA(1,3),U1,1,F,1) ! F3 <- -h*A_11*U1
1787      CALL WAXPY(N,-H*rkA(2,3),U2,1,F,1) ! F3 <- F3 - h*A_21*U2
1788      CALL WAXPY(N,-H*rkA(3,3),U3,1,F,1) ! F3 <- F3 - h*A_31*U3
1789#ifdef FULL_ALGEBRA 
1790      TMP = MATMUL(TRANSPOSE(Jac3),F)   
1791#else     
1792      CALL JacTR_SP_Vec ( Jac3, F, TMP )   ! R3 <- -Jac(Y+Z3)^t*h*sum(A_j3*U_j) 
1793#endif     
1794      CALL WAXPY(N,ONE,U3,1,TMP,1) ! R3 <- U3 -Jac(Y+Z3)^t*h*sum(A_j3*U_j)
1795      CALL WAXPY(N,ONE,TMP,1,R3,1) ! R3 <- U3 -Jac(Y+Z3)^t*h*sum(A_j3*U_j)
1796
1797
1798  END SUBROUTINE RK_PrepareRHS_Adj
1799
1800!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1801   SUBROUTINE RK_PrepareRHS_G(N,H,Jac1,Jac2,Jac3,Lambda, &
1802                      G1,G2,G3)
1803!~~~> Prepare the right-hand side for Newton iterations
1804!     R = Z_adj - hA x Jac*Z_adj - h J^t b lambda
1805!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1806      IMPLICIT NONE
1807     
1808      INTEGER, INTENT(IN) :: N
1809      KPP_REAL, INTENT(IN) :: H
1810      KPP_REAL, DIMENSION(N), INTENT(IN) :: Lambda
1811      KPP_REAL, DIMENSION(N), INTENT(OUT) :: G1,G2,G3
1812#ifdef FULL_ALGEBRA     
1813      KPP_REAL, DIMENSION(N,N), INTENT(IN) :: Jac1, Jac2, Jac3
1814#else     
1815      KPP_REAL, DIMENSION(LU_NONZERO),INTENT(IN) :: Jac1, Jac2, Jac3
1816#endif
1817      KPP_REAL, DIMENSION(N) ::  F 
1818
1819      CALL SET2ZERO(N,G1)
1820      CALL SET2ZERO(N,G2)
1821      CALL SET2ZERO(N,G3)
1822#ifdef FULL_ALGEBRA 
1823      F = MATMUL(TRANSPOSE(Jac1),Lambda)   
1824#else     
1825      CALL JacTR_SP_Vec ( Jac1, Lambda, F )   ! F1 <- Jac(Y+Z1)^t*Lambda 
1826#endif     
1827      CALL WAXPY(N,-H*rkB(1),F,1,G1,1) ! R1 <- R1 - h*B_1*F1
1828
1829#ifdef FULL_ALGEBRA     
1830      F = MATMUL(TRANSPOSE(Jac2),Lambda)   
1831#else     
1832      CALL JacTR_SP_Vec ( Jac2, Lambda, F )   ! F2 <- Jac(Y+Z2)^t*Lambda 
1833#endif     
1834      CALL WAXPY(N,-H*rkB(2),F,1,G2,1) ! R2 <- R2 - h*B_2*F2
1835
1836#ifdef FULL_ALGEBRA     
1837      F = MATMUL(TRANSPOSE(Jac3),Lambda)   
1838#else     
1839      CALL JacTR_SP_Vec ( Jac3, Lambda, F )   ! F3 <- Jac(Y+Z3)^t*Lambda 
1840#endif     
1841      CALL WAXPY(N,-H*rkB(3),F,1,G3,1) ! R3 <- R3 - h*B_3*F3
1842
1843           
1844  END SUBROUTINE RK_PrepareRHS_G
1845 
1846
1847!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1848   SUBROUTINE RK_Decomp(N,H,FJAC,E1,IP1,E2,IP2,ISING)
1849   !~~~> Compute the matrices E1 and E2 and their decompositions
1850!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1851      IMPLICIT NONE
1852     
1853      INTEGER :: N, ISING
1854      KPP_REAL    :: H, Alpha, Beta, Gamma
1855#ifdef FULL_ALGEBRA     
1856      KPP_REAL    :: FJAC(N,N),E1(N,N)
1857      COMPLEX(kind=dp) :: E2(N,N)
1858#else     
1859      KPP_REAL    :: FJAC(LU_NONZERO),E1(LU_NONZERO)
1860      COMPLEX(kind=dp) :: E2(LU_NONZERO)
1861#endif     
1862      INTEGER :: IP1(N), IP2(N), i, j
1863     
1864      Gamma = rkGamma/H
1865      Alpha = rkAlpha/H
1866      Beta  = rkBeta /H
1867
1868#ifdef FULL_ALGEBRA     
1869      DO j=1,N
1870         DO  i=1,N
1871            E1(i,j)=-FJAC(i,j)
1872         END DO
1873         E1(j,j)=E1(j,j)+Gamma
1874      END DO
1875      CALL DGETRF(N,N,E1,N,IP1,ISING) 
1876#else     
1877      DO i=1,LU_NONZERO
1878         E1(i)=-FJAC(i)
1879      END DO
1880      DO i=1,N
1881         j=LU_DIAG(i); E1(j)=E1(j)+Gamma
1882      END DO
1883      CALL KppDecomp(E1,ISING)
1884#endif     
1885     
1886      IF (ISING /= 0) THEN
1887         ISTATUS(idec) = ISTATUS(idec) + 1
1888         RETURN
1889      END IF
1890     
1891#ifdef FULL_ALGEBRA     
1892      DO j=1,N
1893        DO i=1,N
1894          E2(i,j) = DCMPLX( -FJAC(i,j), ZERO )
1895        END DO
1896        E2(j,j) = E2(j,j) + CMPLX( Alpha, Beta )
1897      END DO
1898      CALL ZGETRF(N,N,E2,N,IP2,ISING)     
1899#else 
1900      DO i=1,LU_NONZERO
1901         E2(i) = DCMPLX( -FJAC(i), ZERO )
1902      END DO
1903      DO i=1,N
1904         j = LU_DIAG(i); E2(j)=E2(j) + CMPLX( Alpha, Beta )
1905      END DO
1906      CALL KppDecompCmplx(E2,ISING)   
1907#endif     
1908      ISTATUS(idec) = ISTATUS(idec) + 1
1909     
1910   END SUBROUTINE RK_Decomp
1911
1912
1913!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1914   SUBROUTINE RK_Solve(N,H,E1,IP1,E2,IP2,R1,R2,R3,ISING)
1915!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1916      IMPLICIT NONE
1917      INTEGER :: N,IP1(N),IP2(N),ISING
1918#ifdef FULL_ALGEBRA     
1919      KPP_REAL    :: E1(N,N)
1920      COMPLEX(kind=dp) :: E2(N,N)
1921#else     
1922      KPP_REAL    :: E1(LU_NONZERO)
1923      COMPLEX(kind=dp) :: E2(LU_NONZERO)
1924#endif     
1925      KPP_REAL    :: R1(N),R2(N),R3(N)
1926      KPP_REAL    :: H, x1, x2, x3
1927      COMPLEX(kind=dp) :: BC(N)
1928      INTEGER :: i
1929!     
1930     ! Z <- h^{-1) T^{-1) A^{-1) x Z
1931      DO i=1,N
1932          x1 = R1(i)/H; x2 = R2(i)/H; x3 = R3(i)/H
1933          R1(i) = rkTinvAinv(1,1)*x1 + rkTinvAinv(1,2)*x2 + rkTinvAinv(1,3)*x3
1934          R2(i) = rkTinvAinv(2,1)*x1 + rkTinvAinv(2,2)*x2 + rkTinvAinv(2,3)*x3
1935          R3(i) = rkTinvAinv(3,1)*x1 + rkTinvAinv(3,2)*x2 + rkTinvAinv(3,3)*x3
1936      END DO
1937
1938#ifdef FULL_ALGEBRA     
1939      CALL DGETRS ('N',N,1,E1,N,IP1,R1,N,0) 
1940#else     
1941      CALL KppSolve (E1,R1)
1942#endif     
1943!     
1944      DO i=1,N
1945        BC(i) = DCMPLX(R2(i),R3(i))
1946      END DO
1947#ifdef FULL_ALGEBRA     
1948      CALL ZGETRS ('N',N,1,E2,N,IP2,BC,N,0) 
1949#else     
1950      CALL KppSolveCmplx (E2,BC)
1951#endif     
1952      DO i=1,N
1953        R2(i) = DBLE( BC(i) )
1954        R3(i) = AIMAG( BC(i) )
1955      END DO
1956
1957      ! Z <- T x Z
1958      DO i=1,N
1959          x1 = R1(i); x2 = R2(i); x3 = R3(i)
1960          R1(i) = rkT(1,1)*x1 + rkT(1,2)*x2 + rkT(1,3)*x3
1961          R2(i) = rkT(2,1)*x1 + rkT(2,2)*x2 + rkT(2,3)*x3
1962          R3(i) = rkT(3,1)*x1 + rkT(3,2)*x2 + rkT(3,3)*x3
1963      END DO
1964
1965      ISTATUS(isol) = ISTATUS(isol) + 1
1966
1967   END SUBROUTINE RK_Solve
1968
1969
1970!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1971   SUBROUTINE RK_SolveTR(N,H,E1,IP1,E2,IP2,R1,R2,R3,ISING)
1972!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1973      IMPLICIT NONE
1974      INTEGER, INTENT(IN)  :: N,IP1(N),IP2(N)
1975      INTEGER, INTENT(OUT) :: ISING
1976#ifdef FULL_ALGEBRA     
1977      KPP_REAL,    INTENT(IN) :: E1(N,N)
1978      COMPLEX(kind=dp), INTENT(IN) :: E2(N,N)
1979#else     
1980      KPP_REAL,    INTENT(IN) :: E1(LU_NONZERO)
1981      COMPLEX(kind=dp), INTENT(IN) :: E2(LU_NONZERO)
1982#endif     
1983      KPP_REAL, INTENT(INOUT) :: R1(N),R2(N),R3(N)
1984      KPP_REAL    :: H, x1, x2, x3
1985      COMPLEX(kind=dp) :: BC(N)
1986      INTEGER :: i
1987!     
1988     ! RHS <- h^{-1) (A^{-1) T^{-1))^t x RHS
1989      DO i=1,N
1990          x1 = R1(i)/H; x2 = R2(i)/H; x3 = R3(i)/H
1991          R1(i) = rkAinvT(1,1)*x1 + rkAinvT(2,1)*x2 + rkAinvT(3,1)*x3
1992          R2(i) = rkAinvT(1,2)*x1 + rkAinvT(2,2)*x2 + rkAinvT(3,2)*x3
1993          R3(i) = rkAinvT(1,3)*x1 + rkAinvT(2,3)*x2 + rkAinvT(3,3)*x3
1994      END DO
1995
1996#ifdef FULL_ALGEBRA     
1997      CALL DGETRS ('T',N,1,E1,N,IP1,R1,N,0) 
1998#else     
1999      CALL KppSolveTR (E1,R1,R1)
2000#endif     
2001!     
2002      DO i=1,N
2003        BC(i) = DCMPLX(R2(i),-R3(i))
2004      END DO
2005#ifdef FULL_ALGEBRA     
2006      CALL ZGETRS ('T',N,1,E2,N,IP2,BC,N,0) 
2007#else     
2008      CALL KppSolveTRCmplx (E2,BC)
2009#endif     
2010      DO i=1,N
2011        R2(i) = DBLE( BC(i) )
2012        R3(i) = -AIMAG( BC(i) )
2013      END DO
2014
2015      ! X <- (T^{-1})^t x X
2016      DO i=1,N
2017          x1 = R1(i); x2 = R2(i); x3 = R3(i)
2018          R1(i) = rkTinv(1,1)*x1 + rkTinv(2,1)*x2 + rkTinv(3,1)*x3
2019          R2(i) = rkTinv(1,2)*x1 + rkTinv(2,2)*x2 + rkTinv(3,2)*x3
2020          R3(i) = rkTinv(1,3)*x1 + rkTinv(2,3)*x2 + rkTinv(3,3)*x3
2021      END DO
2022
2023      ISTATUS(isol) = ISTATUS(isol) + 1
2024
2025   END SUBROUTINE RK_SolveTR
2026
2027!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2028   SUBROUTINE RK_ErrorEstimate(N,H,Y,T, &
2029               E1,IP1,Z1,Z2,Z3,SCAL,Err,     &
2030               FirstStep,Reject)
2031!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2032      IMPLICIT NONE
2033     
2034      INTEGER :: N
2035#ifdef FULL_ALGEBRA     
2036      KPP_REAL :: E1(N,N)
2037#else     
2038      KPP_REAL :: E1(LU_NONZERO)
2039#endif     
2040      KPP_REAL :: SCAL(N),Z1(N),Z2(N),Z3(N),F1(N),F2(N), &
2041                        F0(N),Y(N),TMP(N),T,H
2042      INTEGER :: IP1(N), i
2043      LOGICAL FirstStep,Reject
2044      KPP_REAL :: HEE1,HEE2,HEE3,Err
2045
2046      HEE1  = rkE(1)/H
2047      HEE2  = rkE(2)/H
2048      HEE3  = rkE(3)/H
2049
2050      CALL FUN_CHEM(T,Y,F0)
2051      ISTATUS(ifun) = ISTATUS(ifun) + 1
2052
2053      DO  i=1,N
2054         F2(i)  = HEE1*Z1(i)+HEE2*Z2(i)+HEE3*Z3(i)
2055         TMP(i) = rkE(0)*F0(i) + F2(i)
2056      END DO
2057
2058#ifdef FULL_ALGEBRA     
2059      CALL DGETRS ('N',N,1,E1,N,IP1,TMP,N,0) 
2060      IF ((rkMethod==R1A).OR.(rkMethod==GAU).OR.(rkMethod==L3A)) CALL DGETRS ('N',N,1,E1,N,IP1,TMP,N,0)
2061      IF (rkMethod==GAU) CALL DGETRS ('N',N,1,E1,N,IP1,TMP,N,0)
2062#else     
2063      CALL KppSolve (E1, TMP)
2064      IF ((rkMethod==R1A).OR.(rkMethod==GAU).OR.(rkMethod==L3A)) CALL KppSolve (E1,TMP)
2065      IF (rkMethod==GAU) CALL KppSolve (E1,TMP)
2066#endif     
2067
2068      Err = RK_ErrorNorm(N,SCAL,TMP)
2069!
2070      IF (Err < ONE) RETURN
2071firej:IF (FirstStep.OR.Reject) THEN
2072          DO i=1,N
2073             TMP(i)=Y(i)+TMP(i)
2074          END DO
2075          CALL FUN_CHEM(T,TMP,F1)   
2076          ISTATUS(ifun) = ISTATUS(ifun) + 1
2077          DO i=1,N
2078             TMP(i)=F1(i)+F2(i)
2079          END DO
2080
2081#ifdef FULL_ALGEBRA     
2082          CALL DGETRS ('N',N,1,E1,N,IP1,TMP,N,0) 
2083#else     
2084          CALL KppSolve (E1, TMP)
2085#endif     
2086          Err = RK_ErrorNorm(N,SCAL,TMP)
2087       END IF firej
2088 
2089   END SUBROUTINE RK_ErrorEstimate
2090
2091!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2092   KPP_REAL FUNCTION RK_ErrorNorm(N,SCAL,DY)
2093!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2094      IMPLICIT NONE
2095     
2096      INTEGER :: N
2097      KPP_REAL :: SCAL(N),DY(N)
2098      INTEGER :: i
2099
2100      RK_ErrorNorm = ZERO
2101      DO i=1,N
2102          RK_ErrorNorm = RK_ErrorNorm + (DY(i)*SCAL(i))**2
2103      END DO
2104      RK_ErrorNorm = MAX( SQRT(RK_ErrorNorm/N), 1.0d-10 )
2105 
2106   END FUNCTION RK_ErrorNorm
2107 
2108!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2109  SUBROUTINE Radau2A_Coefficients
2110!    The coefficients of the 3-stage Radau-2A method
2111!    (given to ~30 accurate digits)
2112!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
2113      IMPLICIT NONE
2114! The coefficients of the Radau2A method
2115      KPP_REAL :: b0
2116
2117!      b0 = 1.0d0
2118      IF (SdirkError) THEN
2119        b0 = 0.2d-1
2120      ELSE
2121        b0 = 0.5d-1
2122      END IF
2123
2124! The coefficients of the Radau2A method
2125      rkMethod = R2A
2126
2127      rkA(1,1) =  1.968154772236604258683861429918299d-1
2128      rkA(1,2) = -6.55354258501983881085227825696087d-2
2129      rkA(1,3) =  2.377097434822015242040823210718965d-2
2130      rkA(2,1) =  3.944243147390872769974116714584975d-1
2131      rkA(2,2) =  2.920734116652284630205027458970589d-1
2132      rkA(2,3) = -4.154875212599793019818600988496743d-2
2133      rkA(3,1) =  3.764030627004672750500754423692808d-1
2134      rkA(3,2) =  5.124858261884216138388134465196080d-1
2135      rkA(3,3) =  1.111111111111111111111111111111111d-1
2136
2137      rkB(1) = 3.764030627004672750500754423692808d-1
2138      rkB(2) = 5.124858261884216138388134465196080d-1
2139      rkB(3) = 1.111111111111111111111111111111111d-1
2140
2141      rkC(1) = 1.550510257216821901802715925294109d-1
2142      rkC(2) = 6.449489742783178098197284074705891d-1
2143      rkC(3) = 1.0d0
2144     
2145      ! New solution: H* Sum B_j*f(Z_j) = Sum D_j*Z_j
2146      rkD(1) = 0.0d0
2147      rkD(2) = 0.0d0
2148      rkD(3) = 1.0d0
2149
2150      ! Classical error estimator:
2151      ! H* Sum (B_j-Bhat_j)*f(Z_j) = H*E(0)*f(0) + Sum E_j*Z_j
2152      rkE(0) = 1.0d0*b0
2153      rkE(1) = -10.04880939982741556246032950764708d0*b0
2154      rkE(2) = 1.382142733160748895793662840980412d0*b0
2155      rkE(3) = -.3333333333333333333333333333333333d0*b0
2156
2157      ! Sdirk error estimator
2158      rkBgam(0) = b0
2159      rkBgam(1) = .3764030627004672750500754423692807d0-1.558078204724922382431975370686279d0*b0
2160      rkBgam(2) = .8914115380582557157653087040196118d0*b0+.5124858261884216138388134465196077d0
2161      rkBgam(3) = -.1637777184845662566367174924883037d0-.3333333333333333333333333333333333d0*b0
2162      rkBgam(4) = .2748888295956773677478286035994148d0
2163
2164      ! H* Sum Bgam_j*f(Z_j) = H*Bgam(0)*f(0) + Sum Theta_j*Z_j
2165      rkTheta(1) = -1.520677486405081647234271944611547d0-10.04880939982741556246032950764708d0*b0
2166      rkTheta(2) = 2.070455145596436382729929151810376d0+1.382142733160748895793662840980413d0*b0
2167      rkTheta(3) = -.3333333333333333333333333333333333d0*b0-.3744441479783868387391430179970741d0
2168
2169      ! Local order of error estimator
2170      IF (b0==0.0d0) THEN
2171        rkELO  = 6.0d0
2172      ELSE     
2173        rkELO  = 4.0d0
2174      END IF   
2175
2176      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2177      !~~~> Diagonalize the RK matrix:               
2178      ! rkTinv * inv(rkA) * rkT =         
2179      !           |  rkGamma      0           0     |
2180      !           |      0      rkAlpha   -rkBeta   |
2181      !           |      0      rkBeta     rkAlpha  |
2182      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2183
2184      rkGamma = 3.637834252744495732208418513577775d0
2185      rkAlpha = 2.681082873627752133895790743211112d0
2186      rkBeta  = 3.050430199247410569426377624787569d0
2187
2188      rkT(1,1) =  9.443876248897524148749007950641664d-2
2189      rkT(1,2) = -1.412552950209542084279903838077973d-1
2190      rkT(1,3) = -3.00291941051474244918611170890539d-2
2191      rkT(2,1) =  2.502131229653333113765090675125018d-1
2192      rkT(2,2) =  2.041293522937999319959908102983381d-1
2193      rkT(2,3) =  3.829421127572619377954382335998733d-1
2194      rkT(3,1) =  1.0d0
2195      rkT(3,2) =  1.0d0
2196      rkT(3,3) =  0.0d0
2197
2198      rkTinv(1,1) =  4.178718591551904727346462658512057d0
2199      rkTinv(1,2) =  3.27682820761062387082533272429617d-1
2200      rkTinv(1,3) =  5.233764454994495480399309159089876d-1
2201      rkTinv(2,1) = -4.178718591551904727346462658512057d0
2202      rkTinv(2,2) = -3.27682820761062387082533272429617d-1
2203      rkTinv(2,3) =  4.766235545005504519600690840910124d-1
2204      rkTinv(3,1) = -5.02872634945786875951247343139544d-1
2205      rkTinv(3,2) =  2.571926949855605429186785353601676d0
2206      rkTinv(3,3) = -5.960392048282249249688219110993024d-1
2207
2208      rkTinvAinv(1,1) =  1.520148562492775501049204957366528d+1
2209      rkTinvAinv(1,2) =  1.192055789400527921212348994770778d0
2210      rkTinvAinv(1,3) =  1.903956760517560343018332287285119d0
2211      rkTinvAinv(2,1) = -9.669512977505946748632625374449567d0
2212      rkTinvAinv(2,2) = -8.724028436822336183071773193986487d0
2213      rkTinvAinv(2,3) =  3.096043239482439656981667712714881d0
2214      rkTinvAinv(3,1) = -1.409513259499574544876303981551774d+1
2215      rkTinvAinv(3,2) =  5.895975725255405108079130152868952d0
2216      rkTinvAinv(3,3) = -1.441236197545344702389881889085515d-1
2217
2218      rkAinvT(1,1) = .3435525649691961614912493915818282d0
2219      rkAinvT(1,2) = -.4703191128473198422370558694426832d0
2220      rkAinvT(1,3) = .3503786597113668965366406634269080d0
2221      rkAinvT(2,1) = .9102338692094599309122768354288852d0
2222      rkAinvT(2,2) = 1.715425895757991796035292755937326d0
2223      rkAinvT(2,3) = .4040171993145015239277111187301784d0
2224      rkAinvT(3,1) = 3.637834252744495732208418513577775d0
2225      rkAinvT(3,2) = 2.681082873627752133895790743211112d0
2226      rkAinvT(3,3) = -3.050430199247410569426377624787569d0
2227
2228  END SUBROUTINE Radau2A_Coefficients
2229
2230   
2231
2232   
2233!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2234  SUBROUTINE Lobatto3C_Coefficients
2235!    The coefficients of the 3-stage Lobatto-3C method
2236!    (given to ~30 accurate digits)
2237!    The parameter b0 can be chosen to tune the error estimator
2238!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
2239      IMPLICIT NONE
2240      KPP_REAL :: b0
2241
2242      rkMethod = L3C
2243
2244!      b0 = 1.0d0
2245      IF (SdirkError) THEN
2246        b0 = 0.2d0
2247      ELSE
2248        b0 = 0.5d0
2249      END IF
2250! The coefficients of the Lobatto3C method
2251
2252      rkA(1,1) =  .1666666666666666666666666666666667d0
2253      rkA(1,2) = -.3333333333333333333333333333333333d0
2254      rkA(1,3) =  .1666666666666666666666666666666667d0
2255      rkA(2,1) =  .1666666666666666666666666666666667d0
2256      rkA(2,2) =  .4166666666666666666666666666666667d0
2257      rkA(2,3) = -.8333333333333333333333333333333333d-1
2258      rkA(3,1) =  .1666666666666666666666666666666667d0
2259      rkA(3,2) =  .6666666666666666666666666666666667d0
2260      rkA(3,3) =  .1666666666666666666666666666666667d0
2261
2262      rkB(1) = .1666666666666666666666666666666667d0
2263      rkB(2) = .6666666666666666666666666666666667d0
2264      rkB(3) = .1666666666666666666666666666666667d0
2265
2266      rkC(1) = 0.0d0
2267      rkC(2) = 0.5d0
2268      rkC(3) = 1.0d0
2269
2270      ! Classical error estimator, embedded solution:
2271      rkBhat(0) = b0
2272      rkBhat(1) = .16666666666666666666666666666666667d0-b0
2273      rkBhat(2) = .66666666666666666666666666666666667d0
2274      rkBhat(3) = .16666666666666666666666666666666667d0
2275     
2276      ! New solution: h Sum_j b_j f(Z_j) = sum d_j Z_j
2277      rkD(1) = 0.0d0
2278      rkD(2) = 0.0d0
2279      rkD(3) = 1.0d0
2280
2281      ! Classical error estimator:
2282      !   H* Sum (B_j-Bhat_j)*f(Z_j) = H*E(0)*f(0) + Sum E_j*Z_j
2283      rkE(0) =   .3808338772072650364017425226487022*b0
2284      rkE(1) = -1.142501631621795109205227567946107*b0
2285      rkE(2) = -1.523335508829060145606970090594809*b0
2286      rkE(3) =   .3808338772072650364017425226487022*b0
2287
2288      ! Sdirk error estimator
2289      rkBgam(0) = b0
2290      rkBgam(1) = .1666666666666666666666666666666667d0-1.d0*b0
2291      rkBgam(2) = .6666666666666666666666666666666667d0
2292      rkBgam(3) = -.2141672105405983697350758559820354d0
2293      rkBgam(4) = .3808338772072650364017425226487021d0
2294
2295      ! H* Sum Bgam_j*f(Z_j) = H*Bgam(0)*f(0) + Sum Theta_j*Z_j
2296      rkTheta(1) = -3.d0*b0-.3808338772072650364017425226487021d0
2297      rkTheta(2) = -4.d0*b0+1.523335508829060145606970090594808d0
2298      rkTheta(3) = -.142501631621795109205227567946106d0+b0
2299
2300      ! Local order of error estimator
2301      IF (b0==0.0d0) THEN
2302        rkELO  = 5.0d0
2303      ELSE     
2304        rkELO  = 4.0d0
2305      END IF   
2306
2307      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2308      !~~~> Diagonalize the RK matrix:               
2309      ! rkTinv * inv(rkA) * rkT =         
2310      !           |  rkGamma      0           0     |
2311      !           |      0      rkAlpha   -rkBeta   |
2312      !           |      0      rkBeta     rkAlpha  |
2313      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2314
2315      rkGamma = 2.625816818958466716011888933765284d0
2316      rkAlpha = 1.687091590520766641994055533117359d0
2317      rkBeta  = 2.508731754924880510838743672432351d0
2318
2319      rkT(1,1) = 1.d0
2320      rkT(1,2) = 1.d0
2321      rkT(1,3) = 0.d0
2322      rkT(2,1) = .4554100411010284672111720348287483d0
2323      rkT(2,2) = -.6027050205505142336055860174143743d0
2324      rkT(2,3) = -.4309321229203225731070721341350346d0
2325      rkT(3,1) = 2.195823345445647152832799205549709d0
2326      rkT(3,2) = -1.097911672722823576416399602774855d0
2327      rkT(3,3) = .7850032632435902184104551358922130d0
2328
2329      rkTinv(1,1) = .4205559181381766909344950150991349d0
2330      rkTinv(1,2) = .3488903392193734304046467270632057d0
2331      rkTinv(1,3) = .1915253879645878102698098373933487d0
2332      rkTinv(2,1) = .5794440818618233090655049849008650d0
2333      rkTinv(2,2) = -.3488903392193734304046467270632057d0
2334      rkTinv(2,3) = -.1915253879645878102698098373933487d0
2335      rkTinv(3,1) = -.3659705575742745254721332009249516d0
2336      rkTinv(3,2) = -1.463882230297098101888532803699806d0
2337      rkTinv(3,3) = .4702733607340189781407813565524989d0
2338
2339      rkTinvAinv(1,1) = 1.104302803159744452668648155627548d0
2340      rkTinvAinv(1,2) = .916122120694355522658740710823143d0
2341      rkTinvAinv(1,3) = .5029105849749601702795812241441172d0
2342      rkTinvAinv(2,1) = 1.895697196840255547331351844372453d0
2343      rkTinvAinv(2,2) = 3.083877879305644477341259289176857d0
2344      rkTinvAinv(2,3) = -1.502910584974960170279581224144117d0
2345      rkTinvAinv(3,1) = .8362439183082935036129145574774502d0
2346      rkTinvAinv(3,2) = -3.344975673233174014451658229909802d0
2347      rkTinvAinv(3,3) = .312908409479233358005944466882642d0
2348
2349      rkAinvT(1,1) = 2.625816818958466716011888933765282d0
2350      rkAinvT(1,2) = 1.687091590520766641994055533117358d0
2351      rkAinvT(1,3) = -2.508731754924880510838743672432351d0
2352      rkAinvT(2,1) = 1.195823345445647152832799205549710d0
2353      rkAinvT(2,2) = -2.097911672722823576416399602774855d0
2354      rkAinvT(2,3) = .7850032632435902184104551358922130d0
2355      rkAinvT(3,1) = 5.765829871932827589653709477334136d0
2356      rkAinvT(3,2) = .1170850640335862051731452613329320d0
2357      rkAinvT(3,3) = 4.078738281412060947659653944216779d0
2358
2359  END SUBROUTINE Lobatto3C_Coefficients
2360
2361!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2362  SUBROUTINE Gauss_Coefficients
2363!    The coefficients of the 3-stage Gauss method
2364!    (given to ~30 accurate digits)
2365!    The parameter b3 can be chosen by the user
2366!    to tune the error estimator
2367!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
2368      IMPLICIT NONE
2369      KPP_REAL :: b0
2370! The coefficients of the Gauss method
2371
2372
2373      rkMethod = GAU
2374     
2375!      b0 = 4.0d0
2376      b0 = 0.1d0
2377     
2378! The coefficients of the Gauss method
2379
2380      rkA(1,1) =  .1388888888888888888888888888888889d0
2381      rkA(1,2) = -.359766675249389034563954710966045d-1
2382      rkA(1,3) =  .97894440153083260495800422294756d-2
2383      rkA(2,1) =  .3002631949808645924380249472131556d0
2384      rkA(2,2) =  .2222222222222222222222222222222222d0
2385      rkA(2,3) = -.224854172030868146602471694353778d-1
2386      rkA(3,1) =  .2679883337624694517281977355483022d0
2387      rkA(3,2) =  .4804211119693833479008399155410489d0
2388      rkA(3,3) =  .1388888888888888888888888888888889d0
2389
2390      rkB(1) = .2777777777777777777777777777777778d0
2391      rkB(2) = .4444444444444444444444444444444444d0
2392      rkB(3) = .2777777777777777777777777777777778d0
2393
2394      rkC(1) = .1127016653792583114820734600217600d0
2395      rkC(2) = .5000000000000000000000000000000000d0
2396      rkC(3) = .8872983346207416885179265399782400d0
2397
2398      ! Classical error estimator, embedded solution:
2399      rkBhat(0) = b0
2400      rkBhat(1) =-1.4788305577012361475298775666303999d0*b0 &
2401                  +.27777777777777777777777777777777778d0
2402      rkBhat(2) =  .44444444444444444444444444444444444d0 &
2403                  +.66666666666666666666666666666666667d0*b0
2404      rkBhat(3) = -.18783610896543051913678910003626672d0*b0 &
2405                  +.27777777777777777777777777777777778d0
2406
2407      ! New solution: h Sum_j b_j f(Z_j) = sum d_j Z_j
2408      rkD(1) = .1666666666666666666666666666666667d1
2409      rkD(2) = -.1333333333333333333333333333333333d1
2410      rkD(3) = .1666666666666666666666666666666667d1
2411
2412      ! Classical error estimator:
2413      !   H* Sum (B_j-Bhat_j)*f(Z_j) = H*E(0)*f(0) + Sum E_j*Z_j
2414      rkE(0) = .2153144231161121782447335303806954d0*b0
2415      rkE(1) = -2.825278112319014084275808340593191d0*b0
2416      rkE(2) = .2870858974881495709929780405075939d0*b0
2417      rkE(3) = -.4558086256248162565397206448274867d-1*b0
2418
2419      ! Sdirk error estimator
2420      rkBgam(0) = 0.d0
2421      rkBgam(1) = .2373339543355109188382583162660537d0
2422      rkBgam(2) = .5879873931885192299409334646982414d0
2423      rkBgam(3) = -.4063577064014232702392531134499046d-1
2424      rkBgam(4) = .2153144231161121782447335303806955d0
2425
2426      ! H* Sum Bgam_j*f(Z_j) = H*Bgam(0)*f(0) + Sum Theta_j*Z_j
2427      rkTheta(1) = -2.594040933093095272574031876464493d0
2428      rkTheta(2) = 1.824611539036311947589425112250199d0
2429      rkTheta(3) = .1856563166634371860478043996459493d0
2430
2431      ! ELO = local order of classical error estimator
2432      rkELO = 4.0d0
2433
2434      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2435      !~~~> Diagonalize the RK matrix:               
2436      ! rkTinv * inv(rkA) * rkT =         
2437      !           |  rkGamma      0           0     |
2438      !           |      0      rkAlpha   -rkBeta   |
2439      !           |      0      rkBeta     rkAlpha  |
2440      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2441
2442      rkGamma = 4.644370709252171185822941421408064d0
2443      rkAlpha = 3.677814645373914407088529289295970d0
2444      rkBeta  = 3.508761919567443321903661209182446d0
2445
2446      rkT(1,1) =  .7215185205520017032081769924397664d-1
2447      rkT(1,2) = -.8224123057363067064866206597516454d-1
2448      rkT(1,3) = -.6012073861930850173085948921439054d-1
2449      rkT(2,1) =  .1188325787412778070708888193730294d0
2450      rkT(2,2) =  .5306509074206139504614411373957448d-1
2451      rkT(2,3) =  .3162050511322915732224862926182701d0
2452      rkT(3,1) = 1.d0
2453      rkT(3,2) = 1.d0
2454      rkT(3,3) = 0.d0
2455
2456      rkTinv(1,1) =  5.991698084937800775649580743981285d0
2457      rkTinv(1,2) =  1.139214295155735444567002236934009d0
2458      rkTinv(1,3) =   .4323121137838583855696375901180497d0
2459      rkTinv(2,1) = -5.991698084937800775649580743981285d0
2460      rkTinv(2,2) = -1.139214295155735444567002236934009d0
2461      rkTinv(2,3) =   .5676878862161416144303624098819503d0
2462      rkTinv(3,1) = -1.246213273586231410815571640493082d0
2463      rkTinv(3,2) =  2.925559646192313662599230367054972d0
2464      rkTinv(3,3) =  -.2577352012734324923468722836888244d0
2465
2466      rkTinvAinv(1,1) =  27.82766708436744962047620566703329d0
2467      rkTinvAinv(1,2) =   5.290933503982655311815946575100597d0
2468      rkTinvAinv(1,3) =   2.007817718512643701322151051660114d0
2469      rkTinvAinv(2,1) = -17.66368928942422710690385180065675d0
2470      rkTinvAinv(2,2) = -14.45491129892587782538830044147713d0
2471      rkTinvAinv(2,3) =   2.992182281487356298677848948339886d0
2472      rkTinvAinv(3,1) = -25.60678350282974256072419392007303d0
2473      rkTinvAinv(3,2) =   6.762434375611708328910623303779923d0
2474      rkTinvAinv(3,3) =   1.043979339483109825041215970036771d0
2475     
2476      rkAinvT(1,1) = .3350999483034677402618981153470483d0
2477      rkAinvT(1,2) = -.5134173605009692329246186488441294d0
2478      rkAinvT(1,3) = .6745196507033116204327635673208923d-1
2479      rkAinvT(2,1) = .5519025480108928886873752035738885d0
2480      rkAinvT(2,2) = 1.304651810077110066076640761092008d0
2481      rkAinvT(2,3) = .9767507983414134987545585703726984d0
2482      rkAinvT(3,1) = 4.644370709252171185822941421408064d0
2483      rkAinvT(3,2) = 3.677814645373914407088529289295970d0
2484      rkAinvT(3,3) = -3.508761919567443321903661209182446d0
2485     
2486  END SUBROUTINE Gauss_Coefficients
2487
2488
2489
2490!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2491  SUBROUTINE Radau1A_Coefficients
2492!    The coefficients of the 3-stage Gauss method
2493!    (given to ~30 accurate digits)
2494!    The parameter b3 can be chosen by the user
2495!    to tune the error estimator
2496!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
2497      IMPLICIT NONE
2498!      KPP_REAL :: b0 = 0.3d0
2499      KPP_REAL :: b0 = 0.1d0
2500
2501! The coefficients of the Radau1A method
2502
2503      rkMethod = R1A
2504
2505      rkA(1,1) =  .1111111111111111111111111111111111d0
2506      rkA(1,2) = -.1916383190435098943442935597058829d0
2507      rkA(1,3) =  .8052720793239878323318244859477174d-1
2508      rkA(2,1) =  .1111111111111111111111111111111111d0
2509      rkA(2,2) =  .2920734116652284630205027458970589d0
2510      rkA(2,3) = -.481334970546573839513422644787591d-1
2511      rkA(3,1) =  .1111111111111111111111111111111111d0
2512      rkA(3,2) =  .5370223859435462728402311533676479d0
2513      rkA(3,3) =  .1968154772236604258683861429918299d0
2514
2515      rkB(1) = .1111111111111111111111111111111111d0
2516      rkB(2) = .5124858261884216138388134465196080d0
2517      rkB(3) = .3764030627004672750500754423692808d0
2518
2519      rkC(1) = 0.d0
2520      rkC(2) = .3550510257216821901802715925294109d0
2521      rkC(3) = .8449489742783178098197284074705891d0
2522
2523      ! Classical error estimator, embedded solution:
2524      rkBhat(0) = b0
2525      rkBhat(1) = .11111111111111111111111111111111111d0-b0
2526      rkBhat(2) = .51248582618842161383881344651960810d0
2527      rkBhat(3) = .37640306270046727505007544236928079d0
2528
2529      ! New solution: H* Sum B_j*f(Z_j) = Sum D_j*Z_j
2530      rkD(1) = .3333333333333333333333333333333333d0
2531      rkD(2) = -.8914115380582557157653087040196127d0
2532      rkD(3) = .1558078204724922382431975370686279d1
2533
2534      ! Classical error estimator:
2535      ! H* Sum (b_j-bhat_j) f(Z_j) = H*E(0)*F(0) + Sum E_j Z_j
2536      rkE(0) =   .2748888295956773677478286035994148d0*b0
2537      rkE(1) = -1.374444147978386838739143017997074d0*b0
2538      rkE(2) = -1.335337922441686804550326197041126d0*b0
2539      rkE(3) =   .235782604058977333559011782643466d0*b0
2540
2541      ! Sdirk error estimator
2542      rkBgam(0) = 0.0d0
2543      rkBgam(1) = .1948150124588532186183490991130616d-1
2544      rkBgam(2) = .7575249005733381398986810981093584d0
2545      rkBgam(3) = -.518952314149008295083446116200793d-1
2546      rkBgam(4) = .2748888295956773677478286035994148d0
2547
2548      ! H* Sum Bgam_j*f(Z_j) = H*Bgam(0)*f(0) + Sum Theta_j*Z_j
2549      rkTheta(1) = -1.224370034375505083904362087063351d0
2550      rkTheta(2) = .9340045331532641409047527962010133d0
2551      rkTheta(3) = .4656990124352088397561234800640929d0
2552
2553      ! ELO = local order of classical error estimator
2554      rkELO = 4.0d0
2555
2556      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2557      !~~~> Diagonalize the RK matrix:               
2558      ! rkTinv * inv(rkA) * rkT =         
2559      !           |  rkGamma      0           0     |
2560      !           |      0      rkAlpha   -rkBeta   |
2561      !           |      0      rkBeta     rkAlpha  |
2562      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2563
2564      rkGamma = 3.637834252744495732208418513577775d0
2565      rkAlpha = 2.681082873627752133895790743211112d0
2566      rkBeta  = 3.050430199247410569426377624787569d0
2567
2568      rkT(1,1) =  .424293819848497965354371036408369d0
2569      rkT(1,2) = -.3235571519651980681202894497035503d0
2570      rkT(1,3) = -.522137786846287839586599927945048d0
2571      rkT(2,1) =  .57594609499806128896291585429339d-1
2572      rkT(2,2) =  .3148663231849760131614374283783d-2
2573      rkT(2,3) =  .452429247674359778577728510381731d0
2574      rkT(3,1) = 1.d0
2575      rkT(3,2) = 1.d0
2576      rkT(3,3) = 0.d0
2577
2578      rkTinv(1,1) = 1.233523612685027760114769983066164d0
2579      rkTinv(1,2) = 1.423580134265707095505388133369554d0
2580      rkTinv(1,3) = .3946330125758354736049045150429624d0
2581      rkTinv(2,1) = -1.233523612685027760114769983066164d0
2582      rkTinv(2,2) = -1.423580134265707095505388133369554d0
2583      rkTinv(2,3) = .6053669874241645263950954849570376d0
2584      rkTinv(3,1) = -.1484438963257383124456490049673414d0
2585      rkTinv(3,2) = 2.038974794939896109682070471785315d0
2586      rkTinv(3,3) = -.544501292892686735299355831692542d-1
2587
2588      rkTinvAinv(1,1) =  4.487354449794728738538663081025420d0
2589      rkTinvAinv(1,2) =  5.178748573958397475446442544234494d0
2590      rkTinvAinv(1,3) =  1.435609490412123627047824222335563d0
2591      rkTinvAinv(2,1) = -2.854361287939276673073807031221493d0
2592      rkTinvAinv(2,2) = -1.003648660720543859000994063139137d+1
2593      rkTinvAinv(2,3) =  1.789135380979465422050817815017383d0
2594      rkTinvAinv(3,1) = -4.160768067752685525282947313530352d0
2595      rkTinvAinv(3,2) =  1.124128569859216916690209918405860d0
2596      rkTinvAinv(3,3) =  1.700644430961823796581896350418417d0
2597
2598      rkAinvT(1,1) = 1.543510591072668287198054583233180d0
2599      rkAinvT(1,2) = -2.460228411937788329157493833295004d0
2600      rkAinvT(1,3) = -.412906170450356277003910443520499d0
2601      rkAinvT(2,1) = .209519643211838264029272585946993d0
2602      rkAinvT(2,2) = 1.388545667194387164417459732995766d0
2603      rkAinvT(2,3) = 1.20339553005832004974976023130002d0
2604      rkAinvT(3,1) = 3.637834252744495732208418513577775d0
2605      rkAinvT(3,2) = 2.681082873627752133895790743211112d0
2606      rkAinvT(3,3) = -3.050430199247410569426377624787569d0
2607
2608  END SUBROUTINE Radau1A_Coefficients
2609
2610 
2611!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2612  END SUBROUTINE RungeKuttaADJ ! and all its internal procedures
2613!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2614
2615!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2616  SUBROUTINE FUN_CHEM(T, V, FCT)
2617!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2618
2619    USE KPP_ROOT_Parameters
2620    USE KPP_ROOT_Global
2621    USE KPP_ROOT_Function, ONLY: Fun
2622    USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST, Update_PHOTO
2623
2624    IMPLICIT NONE
2625
2626    KPP_REAL :: V(NVAR), FCT(NVAR)
2627    KPP_REAL :: T, Told
2628
2629    Told = TIME
2630    TIME = T
2631    CALL Update_SUN()
2632    CALL Update_RCONST()
2633    CALL Update_PHOTO()
2634    TIME = Told
2635   
2636    CALL Fun(V, FIX, RCONST, FCT)
2637   
2638  END SUBROUTINE FUN_CHEM
2639
2640
2641!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2642  SUBROUTINE JAC_CHEM (T, V, JF)
2643!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2644
2645    USE KPP_ROOT_Parameters
2646    USE KPP_ROOT_Global
2647    USE KPP_ROOT_JacobianSP
2648    USE KPP_ROOT_Jacobian, ONLY: Jac_SP
2649    USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST, Update_PHOTO
2650
2651    IMPLICIT NONE
2652
2653    KPP_REAL :: V(NVAR), T , Told
2654#ifdef FULL_ALGEBRA   
2655    KPP_REAL :: JV(LU_NONZERO), JF(NVAR,NVAR)
2656    INTEGER       :: i, j 
2657#else
2658    KPP_REAL :: JF(LU_NONZERO)
2659#endif   
2660
2661    Told = TIME
2662    TIME = T
2663    CALL Update_SUN()
2664    CALL Update_RCONST()
2665    CALL Update_PHOTO()
2666    TIME = Told
2667   
2668#ifdef FULL_ALGEBRA   
2669    CALL Jac_SP(V, FIX, RCONST, JV)
2670    DO j=1,NVAR
2671      DO i=1,NVAR
2672         JF(i,j) = 0.0d0
2673      END DO
2674    END DO
2675    DO i=1,LU_NONZERO
2676       JF(LU_IROW(i),LU_ICOL(i)) = JV(i)
2677    END DO
2678#else
2679    CALL Jac_SP(V, FIX, RCONST, JF) 
2680#endif   
2681
2682  END SUBROUTINE JAC_CHEM
2683
2684
2685!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2686
2687END MODULE KPP_ROOT_Integrator
Note: See TracBrowser for help on using the repository browser.