source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_passive1/chem_gasphase_mod.f90 @ 3487

Last change on this file since 3487 was 3487, checked in by maronga, 3 years ago

processed comments from a previous revision, adjustments of palmrungui and palm.f90 to release 6.0

File size: 78.4 KB
Line 
1MODULE chem_gasphase_mod
2 
3!   Mechanism: passive1
4!
5!------------------------------------------------------------------------------!
6!
7! ******Module chem_gasphase_mod is automatically generated by kpp4palm ******
8!
9!   *********Please do NOT change this Code,it will be ovewritten *********
10!
11!------------------------------------------------------------------------------!
12! This file was created by KPP (http://people.cs.vt.edu/asandu/Software/Kpp/)
13! and kpp4palm (created by Klaus Ketelsen). kpp4palm is an adapted version
14! of KP4 (Jöckel,P.,Kerkweg,A.,Pozzer,A.,Sander,R.,Tost,H.,Riede,
15! H.,Baumgaertner,A.,Gromov,S.,and Kern,B.,2010: Development cycle 2 of
16! the Modular Earth Submodel System (MESSy2),Geosci. Model Dev.,3,717-752,
17! https://doi.org/10.5194/gmd-3-717-2010). KP4 is part of the Modular Earth
18! Submodel System (MESSy),which is is available under the  GNU General Public
19! License (GPL).
20!
21! KPP is free software; you can redistribute it and/or modify it under the terms
22! of the General Public Licence as published by the Free Software Foundation;
23! either version 2 of the License,or (at your option) any later version.
24! KPP is distributed in the hope that it will be useful,but WITHOUT ANY WARRANTY;
25! without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
26! PURPOSE. See the GNU General Public Licence for more details.
27!
28!------------------------------------------------------------------------------!
29! This file is part of the PALM model system.
30!
31! PALM is free software: you can redistribute it and/or modify it under the
32! terms of the GNU General Public License as published by the Free Software
33! Foundation,either version 3 of the License,or (at your option) any later
34! version.
35!
36! PALM is distributed in the hope that it will be useful,but WITHOUT ANY
37! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
38! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
39!
40! You should have received a copy of the GNU General Public License along with
41! PALM. If not,see <http://www.gnu.org/licenses/>.
42!
43! Copyright 1997-2018 Leibniz Universitaet Hannover
44!--------------------------------------------------------------------------------!
45!
46! Current revisions:
47! ------------------
48!
49!
50! Former revisions:
51! -----------------
52! $Id: module_header 2460 2017-09-13 14:47:48Z forkel $
53! forkel June 2018: qvap,fakt added
54! forkel June 2018: reset case in  Initialize,Integrate,Update_rconst
55!
56!
57! 2460 2017-09-13 14:47:48Z forkel
58!
59! forkel Sept. 2017: Variables for photolyis added
60!
61!
62! Nov. 2016: Intial version (Klaus Ketelsen)
63!
64!------------------------------------------------------------------------------!
65!
66
67
68! Set kpp Double Precision to PALM Default Precision
69
70  USE kinds,           ONLY: dp=>wp
71
72  USE pegrid,          ONLY: myid, threads_per_task
73
74  IMPLICIT        NONE
75  PRIVATE
76  !SAVE  ! note: occurs again in automatically generated code ...
77
78!  PUBLIC :: IERR_NAMES
79 
80! PUBLIC :: SPC_NAMES,EQN_NAMES,EQN_TAGS,REQ_HET,REQ_AEROSOL,REQ_PHOTRAT &
81!         ,REQ_MCFCT,IP_MAX,jname
82
83  PUBLIC :: eqn_names, phot_names, spc_names
84  PUBLIC :: nmaxfixsteps
85  PUBLIC :: atol, rtol
86  PUBLIC :: nspec, nreact
87  PUBLIC :: temp
88  PUBLIC :: qvap
89  PUBLIC :: fakt
90  PUBLIC :: phot
91  PUBLIC :: rconst
92  PUBLIC :: nvar
93  PUBLIC :: nphot
94  PUBLIC :: vl_dim                     ! PUBLIC to ebable other MODULEs to distiguish between scalar and vec
95 
96  PUBLIC :: initialize, integrate, update_rconst
97  PUBLIC :: chem_gasphase_integrate
98  PUBLIC :: initialize_kpp_ctrl
99
100! END OF MODULE HEADER TEMPLATE
101                                                                 
102! Variables used for vector mode                                 
103                                                                 
104  LOGICAL, PARAMETER          :: l_vector = .FALSE.           
105  INTEGER, PARAMETER          :: i_lu_di = 2
106  INTEGER, PARAMETER          :: vl_dim = 1
107  INTEGER                     :: vl                             
108                                                                 
109  INTEGER                     :: vl_glo                         
110  INTEGER                     :: is, ie                           
111                                                                 
112                                                                 
113  INTEGER, DIMENSION(vl_dim)  :: kacc, krej                       
114  INTEGER, DIMENSION(vl_dim)  :: ierrv                           
115  LOGICAL                     :: data_loaded = .FALSE.             
116! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
117!
118! Parameter Module File
119!
120! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
121!       (http://www.cs.vt.edu/~asandu/Software/KPP)
122! KPP is distributed under GPL,the general public licence
123!       (http://www.gnu.org/copyleft/gpl.html)
124! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
125! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
126!     With important contributions from:
127!        M. Damian,Villanova University,USA
128!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
129!
130! File                 : chem_gasphase_mod_Parameters.f90
131! Time                 : Tue Sep 25 18:34:54 2018
132! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
133! Equation file        : chem_gasphase_mod.kpp
134! Output root filename : chem_gasphase_mod
135!
136! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
137
138
139
140
141
142
143! NSPEC - Number of chemical species
144  INTEGER, PARAMETER :: nspec = 1 
145! NVAR - Number of Variable species
146  INTEGER, PARAMETER :: nvar = 1 
147! NVARACT - Number of Active species
148  INTEGER, PARAMETER :: nvaract = 1 
149! NFIX - Number of Fixed species
150  INTEGER, PARAMETER :: nfix = 1 
151! NREACT - Number of reactions
152  INTEGER, PARAMETER :: nreact = 1 
153! NVARST - Starting of variables in conc. vect.
154  INTEGER, PARAMETER :: nvarst = 1 
155! NFIXST - Starting of fixed in conc. vect.
156  INTEGER, PARAMETER :: nfixst = 2 
157! NONZERO - Number of nonzero entries in Jacobian
158  INTEGER, PARAMETER :: nonzero = 1 
159! LU_NONZERO - Number of nonzero entries in LU factoriz. of Jacobian
160  INTEGER, PARAMETER :: lu_nonzero = 1 
161! CNVAR - (NVAR+1) Number of elements in compressed row format
162  INTEGER, PARAMETER :: cnvar = 2 
163! CNEQN - (NREACT+1) Number stoicm elements in compressed col format
164  INTEGER, PARAMETER :: cneqn = 2 
165! NHESS - Length of Sparse Hessian
166  INTEGER, PARAMETER :: nhess = 1 
167! NMASS - Number of atoms to check mass balance
168  INTEGER, PARAMETER :: nmass = 1 
169
170! Index declaration for variable species in C and VAR
171!   VAR(ind_spc) = C(ind_spc)
172
173  INTEGER, PARAMETER, PUBLIC :: ind_pm10 = 1 
174
175! Index declaration for fixed species in C
176!   C(ind_spc)
177
178
179! Index declaration for fixed species in FIX
180!    FIX(indf_spc) = C(ind_spc) = C(NVAR+indf_spc)
181
182
183! NJVRP - Length of sparse Jacobian JVRP
184  INTEGER, PARAMETER :: njvrp = 1 
185
186! NSTOICM - Length of Sparse Stoichiometric Matrix
187  INTEGER, PARAMETER :: nstoicm = 1 
188
189
190! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
191!
192! Global Data Module File
193!
194! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
195!       (http://www.cs.vt.edu/~asandu/Software/KPP)
196! KPP is distributed under GPL,the general public licence
197!       (http://www.gnu.org/copyleft/gpl.html)
198! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
199! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
200!     With important contributions from:
201!        M. Damian,Villanova University,USA
202!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
203!
204! File                 : chem_gasphase_mod_Global.f90
205! Time                 : Tue Sep 25 18:34:54 2018
206! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
207! Equation file        : chem_gasphase_mod.kpp
208! Output root filename : chem_gasphase_mod
209!
210! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
211
212
213
214
215
216
217! Declaration of global variables
218
219! C - Concentration of all species
220  REAL(kind=dp):: c(nspec)
221! VAR - Concentrations of variable species (global)
222  REAL(kind=dp):: var(nvar)
223! FIX - Concentrations of fixed species (global)
224  REAL(kind=dp):: fix(nfix)
225! VAR,FIX are chunks of array C
226      EQUIVALENCE( c(1), var(1))
227! RCONST - Rate constants (global)
228  REAL(kind=dp):: rconst(nreact)
229! TIME - Current integration time
230  REAL(kind=dp):: time
231! TEMP - Temperature
232  REAL(kind=dp):: temp
233! TSTART - Integration start time
234  REAL(kind=dp):: tstart
235! ATOL - Absolute tolerance
236  REAL(kind=dp):: atol(nvar)
237! RTOL - Relative tolerance
238  REAL(kind=dp):: rtol(nvar)
239! STEPMIN - Lower bound for integration step
240  REAL(kind=dp):: stepmin
241! CFACTOR - Conversion factor for concentration units
242  REAL(kind=dp):: cfactor
243
244! INLINED global variable declarations
245
246! QVAP - Water vapor
247  REAL(kind=dp):: qvap
248! FAKT - Conversion factor
249  REAL(kind=dp):: fakt
250
251
252! INLINED global variable declarations
253
254
255
256! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
257!
258! Sparse Jacobian Data Structures File
259!
260! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
261!       (http://www.cs.vt.edu/~asandu/Software/KPP)
262! KPP is distributed under GPL,the general public licence
263!       (http://www.gnu.org/copyleft/gpl.html)
264! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
265! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
266!     With important contributions from:
267!        M. Damian,Villanova University,USA
268!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
269!
270! File                 : chem_gasphase_mod_JacobianSP.f90
271! Time                 : Tue Sep 25 18:34:54 2018
272! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
273! Equation file        : chem_gasphase_mod.kpp
274! Output root filename : chem_gasphase_mod
275!
276! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
277
278
279
280
281
282
283! Sparse Jacobian Data
284
285
286  INTEGER, PARAMETER, DIMENSION(1):: lu_irow =  (/ &
287       1 /) 
288
289  INTEGER, PARAMETER, DIMENSION(1):: lu_icol =  (/ &
290       1 /) 
291
292  INTEGER, PARAMETER, DIMENSION(2):: lu_crow =  (/ &
293       1, 2 /) 
294
295  INTEGER, PARAMETER, DIMENSION(2):: lu_diag =  (/ &
296       1, 2 /) 
297
298
299
300! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
301!
302! Utility Data Module File
303!
304! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
305!       (http://www.cs.vt.edu/~asandu/Software/KPP)
306! KPP is distributed under GPL,the general public licence
307!       (http://www.gnu.org/copyleft/gpl.html)
308! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
309! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
310!     With important contributions from:
311!        M. Damian,Villanova University,USA
312!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
313!
314! File                 : chem_gasphase_mod_Monitor.f90
315! Time                 : Tue Sep 25 18:34:54 2018
316! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
317! Equation file        : chem_gasphase_mod.kpp
318! Output root filename : chem_gasphase_mod
319!
320! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
321
322
323
324
325
326  CHARACTER(len=15), PARAMETER, DIMENSION(1):: spc_names =  (/ &
327     'PM10           ' /)
328
329  CHARACTER(len=100), PARAMETER, DIMENSION(1):: eqn_names =  (/ &
330     'PM10 --> PM10                                                                                       ' /)
331
332! INLINED global variables
333
334  !   inline f90_data: declaration of global variables for photolysis
335  !   REAL(kind=dp):: phot(nphot)must eventually be moved to global later for
336  INTEGER, PARAMETER :: nphot = 1
337  !   phot photolysis frequencies
338  REAL(kind=dp):: phot(nphot)
339
340  INTEGER, PARAMETER, PUBLIC :: j_no2 = 1
341
342  CHARACTER(len=15), PARAMETER, DIMENSION(nphot):: phot_names =   (/ &
343     'J_NO2          '/)
344
345! End INLINED global variables
346
347
348! Automatic generated PUBLIC Statements for ip_ and ihs_ variables
349 
350! Automatic generated PUBLIC Statements for ip_ and ihs_ variables
351 
352! Automatic generated PUBLIC Statements for ip_ and ihs_ variables
353 
354! Automatic generated PUBLIC Statements for ip_ and ihs_ variables
355 
356 
357!  variable definations from  individual module headers
358 
359! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
360!
361! Initialization File
362!
363! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
364!       (http://www.cs.vt.edu/~asandu/Software/KPP)
365! KPP is distributed under GPL,the general public licence
366!       (http://www.gnu.org/copyleft/gpl.html)
367! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
368! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
369!     With important contributions from:
370!        M. Damian,Villanova University,USA
371!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
372!
373! File                 : chem_gasphase_mod_Initialize.f90
374! Time                 : Tue Sep 25 18:34:54 2018
375! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
376! Equation file        : chem_gasphase_mod.kpp
377! Output root filename : chem_gasphase_mod
378!
379! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
380
381
382
383
384
385! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
386!
387! Numerical Integrator (Time-Stepping) File
388!
389! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
390!       (http://www.cs.vt.edu/~asandu/Software/KPP)
391! KPP is distributed under GPL,the general public licence
392!       (http://www.gnu.org/copyleft/gpl.html)
393! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
394! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
395!     With important contributions from:
396!        M. Damian,Villanova University,USA
397!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
398!
399! File                 : chem_gasphase_mod_Integrator.f90
400! Time                 : Tue Sep 25 18:34:54 2018
401! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
402! Equation file        : chem_gasphase_mod.kpp
403! Output root filename : chem_gasphase_mod
404!
405! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
406
407
408
409! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
410!
411! INTEGRATE - Integrator routine
412!   Arguments :
413!      TIN       - Start Time for Integration
414!      TOUT      - End Time for Integration
415!
416! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
417
418!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
419!  Rosenbrock - Implementation of several Rosenbrock methods:             !
420!               *Ros2                                                    !
421!               *Ros3                                                    !
422!               *Ros4                                                    !
423!               *Rodas3                                                  !
424!               *Rodas4                                                  !
425!  By default the code employs the KPP sparse linear algebra routines     !
426!  Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK)        !
427!                                                                         !
428!    (C)  Adrian Sandu,August 2004                                       !
429!    Virginia Polytechnic Institute and State University                  !
430!    Contact: sandu@cs.vt.edu                                             !
431!    Revised by Philipp Miehe and Adrian Sandu,May 2006                  !                               !
432!    This implementation is part of KPP - the Kinetic PreProcessor        !
433!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
434
435
436  SAVE
437 
438!~~~>  statistics on the work performed by the rosenbrock method
439  INTEGER, PARAMETER :: nfun=1, njac=2, nstp=3, nacc=4, &
440                        nrej=5, ndec=6, nsol=7, nsng=8, &
441                        ntexit=1, nhexit=2, nhnew = 3
442
443! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
444!
445! Linear Algebra Data and Routines File
446!
447! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
448!       (http://www.cs.vt.edu/~asandu/Software/KPP)
449! KPP is distributed under GPL,the general public licence
450!       (http://www.gnu.org/copyleft/gpl.html)
451! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
452! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
453!     With important contributions from:
454!        M. Damian,Villanova University,USA
455!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
456!
457! File                 : chem_gasphase_mod_LinearAlgebra.f90
458! Time                 : Tue Sep 25 18:34:54 2018
459! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
460! Equation file        : chem_gasphase_mod.kpp
461! Output root filename : chem_gasphase_mod
462!
463! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
464
465
466
467
468
469
470! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
471!
472! The ODE Jacobian of Chemical Model File
473!
474! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
475!       (http://www.cs.vt.edu/~asandu/Software/KPP)
476! KPP is distributed under GPL,the general public licence
477!       (http://www.gnu.org/copyleft/gpl.html)
478! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
479! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
480!     With important contributions from:
481!        M. Damian,Villanova University,USA
482!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
483!
484! File                 : chem_gasphase_mod_Jacobian.f90
485! Time                 : Tue Sep 25 18:34:54 2018
486! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
487! Equation file        : chem_gasphase_mod.kpp
488! Output root filename : chem_gasphase_mod
489!
490! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
491
492
493
494
495
496
497! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
498!
499! The ODE Function of Chemical Model File
500!
501! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
502!       (http://www.cs.vt.edu/~asandu/Software/KPP)
503! KPP is distributed under GPL,the general public licence
504!       (http://www.gnu.org/copyleft/gpl.html)
505! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
506! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
507!     With important contributions from:
508!        M. Damian,Villanova University,USA
509!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
510!
511! File                 : chem_gasphase_mod_Function.f90
512! Time                 : Tue Sep 25 18:34:54 2018
513! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
514! Equation file        : chem_gasphase_mod.kpp
515! Output root filename : chem_gasphase_mod
516!
517! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
518
519
520
521
522
523! A - Rate for each equation
524  REAL(kind=dp):: a(nreact)
525
526! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
527!
528! The Reaction Rates File
529!
530! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
531!       (http://www.cs.vt.edu/~asandu/Software/KPP)
532! KPP is distributed under GPL,the general public licence
533!       (http://www.gnu.org/copyleft/gpl.html)
534! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
535! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
536!     With important contributions from:
537!        M. Damian,Villanova University,USA
538!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
539!
540! File                 : chem_gasphase_mod_Rates.f90
541! Time                 : Tue Sep 25 18:34:54 2018
542! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
543! Equation file        : chem_gasphase_mod.kpp
544! Output root filename : chem_gasphase_mod
545!
546! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
547
548
549
550
551
552! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
553!
554! Auxiliary Routines File
555!
556! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
557!       (http://www.cs.vt.edu/~asandu/Software/KPP)
558! KPP is distributed under GPL,the general public licence
559!       (http://www.gnu.org/copyleft/gpl.html)
560! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
561! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
562!     With important contributions from:
563!        M. Damian,Villanova University,USA
564!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
565!
566! File                 : chem_gasphase_mod_Util.f90
567! Time                 : Tue Sep 25 18:34:54 2018
568! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
569! Equation file        : chem_gasphase_mod.kpp
570! Output root filename : chem_gasphase_mod
571!
572! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
573
574
575
576
577
578
579  ! header MODULE initialize_kpp_ctrl_template
580
581  ! notes:
582  ! - l_vector is automatically defined by kp4
583  ! - vl_dim is automatically defined by kp4
584  ! - i_lu_di is automatically defined by kp4
585  ! - wanted is automatically defined by xmecca
586  ! - icntrl rcntrl are automatically defined by kpp
587  ! - "USE messy_main_tools" is in MODULE_header of messy_mecca_kpp.f90
588  ! - SAVE will be automatically added by kp4
589
590  !SAVE
591
592  ! for fixed time step control
593  ! ... max. number of fixed time steps (sum must be 1)
594  INTEGER, PARAMETER         :: nmaxfixsteps = 50
595  ! ... switch for fixed time stepping
596  LOGICAL, PUBLIC            :: l_fixed_step = .FALSE.
597  INTEGER, PUBLIC            :: nfsteps = 1
598  ! ... number of kpp control PARAMETERs
599  INTEGER, PARAMETER, PUBLIC :: nkppctrl = 20
600  !
601  INTEGER, DIMENSION(nkppctrl), PUBLIC     :: icntrl = 0
602  REAL(dp), DIMENSION(nkppctrl), PUBLIC     :: rcntrl = 0.0_dp
603  REAL(dp), DIMENSION(nmaxfixsteps), PUBLIC :: t_steps = 0.0_dp
604
605  ! END header MODULE initialize_kpp_ctrl_template
606
607 
608! Interface Block
609 
610  INTERFACE            initialize
611    MODULE PROCEDURE   initialize
612  END INTERFACE        initialize
613 
614  INTERFACE            integrate
615    MODULE PROCEDURE   integrate
616  END INTERFACE        integrate
617 
618  INTERFACE            fun
619    MODULE PROCEDURE   fun
620  END INTERFACE        fun
621 
622  INTERFACE            kppsolve
623    MODULE PROCEDURE   kppsolve
624  END INTERFACE        kppsolve
625 
626  INTERFACE            jac_sp
627    MODULE PROCEDURE   jac_sp
628  END INTERFACE        jac_sp
629 
630  INTERFACE            k_arr
631    MODULE PROCEDURE   k_arr
632  END INTERFACE        k_arr
633 
634  INTERFACE            update_rconst
635    MODULE PROCEDURE   update_rconst
636  END INTERFACE        update_rconst
637 
638  INTERFACE            arr2
639    MODULE PROCEDURE   arr2
640  END INTERFACE        arr2
641 
642  INTERFACE            initialize_kpp_ctrl
643    MODULE PROCEDURE   initialize_kpp_ctrl
644  END INTERFACE        initialize_kpp_ctrl
645 
646  INTERFACE            error_output
647    MODULE PROCEDURE   error_output
648  END INTERFACE        error_output
649 
650  INTERFACE            wscal
651    MODULE PROCEDURE   wscal
652  END INTERFACE        wscal
653 
654!INTERFACE not working  INTERFACE            waxpy
655!INTERFACE not working    MODULE PROCEDURE   waxpy
656!INTERFACE not working  END INTERFACE        waxpy
657 
658  INTERFACE            rosenbrock
659    MODULE PROCEDURE   rosenbrock
660  END INTERFACE        rosenbrock
661 
662  INTERFACE            funtemplate
663    MODULE PROCEDURE   funtemplate
664  END INTERFACE        funtemplate
665 
666  INTERFACE            jactemplate
667    MODULE PROCEDURE   jactemplate
668  END INTERFACE        jactemplate
669 
670  INTERFACE            kppdecomp
671    MODULE PROCEDURE   kppdecomp
672  END INTERFACE        kppdecomp
673 
674  INTERFACE            chem_gasphase_integrate
675    MODULE PROCEDURE   chem_gasphase_integrate
676  END INTERFACE        chem_gasphase_integrate
677 
678 
679 CONTAINS
680 
681SUBROUTINE initialize()
682
683
684  INTEGER         :: j, k
685
686  INTEGER :: i
687  REAL(kind=dp):: x
688  k = is
689  cfactor = 1.000000e+00_dp
690
691  x = (0.) * cfactor
692  DO i = 1 , nvar
693  ENDDO
694
695  x = (0.) * cfactor
696  DO i = 1 , nfix
697    fix(i) = x
698  ENDDO
699
700! constant rate coefficients
701! END constant rate coefficients
702
703! INLINED initializations
704
705! End INLINED initializations
706
707     
708END SUBROUTINE initialize
709 
710SUBROUTINE integrate( tin, tout, &
711  icntrl_u, rcntrl_u, istatus_u, rstatus_u, ierr_u)
712
713
714   REAL(kind=dp), INTENT(IN):: tin  ! start time
715   REAL(kind=dp), INTENT(IN):: tout ! END time
716   ! OPTIONAL input PARAMETERs and statistics
717   INTEGER,      INTENT(IN), OPTIONAL :: icntrl_u(20)
718   REAL(kind=dp), INTENT(IN), OPTIONAL :: rcntrl_u(20)
719   INTEGER,      INTENT(OUT), OPTIONAL :: istatus_u(20)
720   REAL(kind=dp), INTENT(OUT), OPTIONAL :: rstatus_u(20)
721   INTEGER,      INTENT(OUT), OPTIONAL :: ierr_u
722
723   REAL(kind=dp):: rcntrl(20), rstatus(20)
724   INTEGER       :: icntrl(20), istatus(20), ierr
725
726   INTEGER, SAVE :: ntotal = 0
727
728   icntrl(:) = 0
729   rcntrl(:) = 0.0_dp
730   istatus(:) = 0
731   rstatus(:) = 0.0_dp
732
733    !~~~> fine-tune the integrator:
734   icntrl(1) = 0      ! 0 - non- autonomous, 1 - autonomous
735   icntrl(2) = 0      ! 0 - vector tolerances, 1 - scalars
736
737   ! IF OPTIONAL PARAMETERs are given, and IF they are >0,
738   ! THEN they overwrite default settings.
739   IF (PRESENT(icntrl_u))THEN
740     WHERE(icntrl_u(:)> 0)icntrl(:) = icntrl_u(:)
741   ENDIF
742   IF (PRESENT(rcntrl_u))THEN
743     WHERE(rcntrl_u(:)> 0)rcntrl(:) = rcntrl_u(:)
744   ENDIF
745
746
747   CALL rosenbrock(nvar, var, tin, tout,  &
748         atol, rtol,               &
749         rcntrl, icntrl, rstatus, istatus, ierr)
750
751   !~~~> debug option: show no of steps
752   ! ntotal = ntotal + istatus(nstp)
753   ! PRINT*,'NSTEPS=',ISTATUS(Nstp),' (',Ntotal,')','  O3=',VAR(ind_O3)
754
755   stepmin = rstatus(nhexit)
756   ! IF OPTIONAL PARAMETERs are given for output they
757   ! are updated with the RETURN information
758   IF (PRESENT(istatus_u))istatus_u(:) = istatus(:)
759   IF (PRESENT(rstatus_u))rstatus_u(:) = rstatus(:)
760   IF (PRESENT(ierr_u))  ierr_u       = ierr
761
762END SUBROUTINE integrate
763 
764SUBROUTINE fun(v, f, rct, vdot)
765
766! V - Concentrations of variable species (local)
767  REAL(kind=dp):: v(nvar)
768! F - Concentrations of fixed species (local)
769  REAL(kind=dp):: f(nfix)
770! RCT - Rate constants (local)
771  REAL(kind=dp):: rct(nreact)
772! Vdot - Time derivative of variable species concentrations
773  REAL(kind=dp):: vdot(nvar)
774
775
776! Computation of equation rates
777
778! Aggregate function
779  vdot(1) = 0
780     
781END SUBROUTINE fun
782 
783SUBROUTINE kppsolve(jvs, x)
784
785! JVS - sparse Jacobian of variables
786  REAL(kind=dp):: jvs(lu_nonzero)
787! X - Vector for variables
788  REAL(kind=dp):: x(nvar)
789
790  x(1) = x(1) / jvs(1)
791     
792END SUBROUTINE kppsolve
793 
794SUBROUTINE jac_sp(v, f, rct, jvs)
795
796! V - Concentrations of variable species (local)
797  REAL(kind=dp):: v(nvar)
798! F - Concentrations of fixed species (local)
799  REAL(kind=dp):: f(nfix)
800! RCT - Rate constants (local)
801  REAL(kind=dp):: rct(nreact)
802! JVS - sparse Jacobian of variables
803  REAL(kind=dp):: jvs(lu_nonzero)
804
805
806! Local variables
807! B - Temporary array
808  REAL(kind=dp):: b(1)
809
810! B(1) = dA(1)/dV(1)
811  b(1) = rct(1)
812
813! Construct the Jacobian terms from B's
814! JVS(1) = Jac_FULL(1,1)
815  jvs(1) = 0
816     
817END SUBROUTINE jac_sp
818 
819  elemental REAL(kind=dp)FUNCTION k_arr (k_298, tdep, temp)
820    ! arrhenius FUNCTION
821
822    REAL,    INTENT(IN):: k_298 ! k at t = 298.15k
823    REAL,    INTENT(IN):: tdep  ! temperature dependence
824    REAL(kind=dp), INTENT(IN):: temp  ! temperature
825
826    intrinsic exp
827
828    k_arr = k_298 * exp(tdep* (1._dp/temp- 3.3540e-3_dp))! 1/298.15=3.3540e-3
829
830  END FUNCTION k_arr
831 
832SUBROUTINE update_rconst()
833 INTEGER         :: k
834
835  k = is
836
837! Begin INLINED RCONST
838
839
840! End INLINED RCONST
841
842  rconst(1) = (1.0_dp)
843     
844END SUBROUTINE update_rconst
845 
846!  END FUNCTION ARR2
847REAL(kind=dp)FUNCTION arr2( a0, b0, temp)
848    REAL(kind=dp):: temp
849    REAL(kind=dp):: a0, b0
850    arr2 = a0 * exp( - b0 / temp)
851END FUNCTION arr2
852 
853SUBROUTINE initialize_kpp_ctrl(status)
854
855
856  ! i/o
857  INTEGER,         INTENT(OUT):: status
858
859  ! local
860  REAL(dp):: tsum
861  INTEGER  :: i
862
863  ! check fixed time steps
864  tsum = 0.0_dp
865  DO i=1, nmaxfixsteps
866     IF (t_steps(i)< tiny(0.0_dp))exit
867     tsum = tsum + t_steps(i)
868  ENDDO
869
870  nfsteps = i- 1
871
872  l_fixed_step = (nfsteps > 0).and.((tsum - 1.0)< tiny(0.0_dp))
873
874  IF (l_vector)THEN
875     WRITE(*,*) ' MODE             : VECTOR (LENGTH=',VL_DIM,')'
876  ELSE
877     WRITE(*,*) ' MODE             : SCALAR'
878  ENDIF
879  !
880  WRITE(*,*) ' DE-INDEXING MODE :',I_LU_DI
881  !
882  WRITE(*,*) ' ICNTRL           : ',icntrl
883  WRITE(*,*) ' RCNTRL           : ',rcntrl
884  !
885  ! note: this is ONLY meaningful for vectorized (kp4)rosenbrock- methods
886  IF (l_vector)THEN
887     IF (l_fixed_step)THEN
888        WRITE(*,*) ' TIME STEPS       : FIXED (',t_steps(1:nfsteps),')'
889     ELSE
890        WRITE(*,*) ' TIME STEPS       : AUTOMATIC'
891     ENDIF
892  ELSE
893     WRITE(*,*) ' TIME STEPS       : AUTOMATIC '//&
894          &'(t_steps (CTRL_KPP) ignored in SCALAR MODE)'
895  ENDIF
896  ! mz_pj_20070531-
897
898  status = 0
899
900
901END SUBROUTINE initialize_kpp_ctrl
902 
903SUBROUTINE error_output(c, ierr, pe)
904
905
906  INTEGER, INTENT(IN):: ierr
907  INTEGER, INTENT(IN):: pe
908  REAL(dp), DIMENSION(:), INTENT(IN):: c
909
910  write(6,*) 'ERROR in chem_gasphase_mod ',ierr,C(1)
911
912
913END SUBROUTINE error_output
914 
915      SUBROUTINE wscal(n, alpha, x, incx)
916!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
917!     constant times a vector: x(1:N) <- Alpha*x(1:N)
918!     only for incX=incY=1
919!     after BLAS
920!     replace this by the function from the optimized BLAS implementation:
921!         CALL SSCAL(N,Alpha,X,1) or  CALL DSCAL(N,Alpha,X,1)
922!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
923
924      INTEGER  :: i, incx, m, mp1, n
925      REAL(kind=dp) :: x(n), alpha
926      REAL(kind=dp), PARAMETER  :: zero=0.0_dp, one=1.0_dp
927
928      IF (alpha .eq. one)RETURN
929      IF (n .le. 0)RETURN
930
931      m = mod(n, 5)
932      IF ( m .ne. 0)THEN
933        IF (alpha .eq. (- one))THEN
934          DO i = 1, m
935            x(i) = - x(i)
936          ENDDO
937        ELSEIF (alpha .eq. zero)THEN
938          DO i = 1, m
939            x(i) = zero
940          ENDDO
941        ELSE
942          DO i = 1, m
943            x(i) = alpha* x(i)
944          ENDDO
945        ENDIF
946        IF ( n .lt. 5)RETURN
947      ENDIF
948      mp1 = m + 1
949      IF (alpha .eq. (- one))THEN
950        DO i = mp1, n, 5
951          x(i)   = - x(i)
952          x(i + 1) = - x(i + 1)
953          x(i + 2) = - x(i + 2)
954          x(i + 3) = - x(i + 3)
955          x(i + 4) = - x(i + 4)
956        ENDDO
957      ELSEIF (alpha .eq. zero)THEN
958        DO i = mp1, n, 5
959          x(i)   = zero
960          x(i + 1) = zero
961          x(i + 2) = zero
962          x(i + 3) = zero
963          x(i + 4) = zero
964        ENDDO
965      ELSE
966        DO i = mp1, n, 5
967          x(i)   = alpha* x(i)
968          x(i + 1) = alpha* x(i + 1)
969          x(i + 2) = alpha* x(i + 2)
970          x(i + 3) = alpha* x(i + 3)
971          x(i + 4) = alpha* x(i + 4)
972        ENDDO
973      ENDIF
974
975      END SUBROUTINE wscal
976 
977      SUBROUTINE waxpy(n, alpha, x, incx, y, incy)
978!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
979!     constant times a vector plus a vector: y <- y + Alpha*x
980!     only for incX=incY=1
981!     after BLAS
982!     replace this by the function from the optimized BLAS implementation:
983!         CALL SAXPY(N,Alpha,X,1,Y,1) or  CALL DAXPY(N,Alpha,X,1,Y,1)
984!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
985
986      INTEGER  :: i, incx, incy, m, mp1, n
987      REAL(kind=dp):: x(n), y(n), alpha
988      REAL(kind=dp), PARAMETER :: zero = 0.0_dp
989
990      IF (alpha .eq. zero)RETURN
991      IF (n .le. 0)RETURN
992
993      m = mod(n, 4)
994      IF ( m .ne. 0)THEN
995        DO i = 1, m
996          y(i) = y(i) + alpha* x(i)
997        ENDDO
998        IF ( n .lt. 4)RETURN
999      ENDIF
1000      mp1 = m + 1
1001      DO i = mp1, n, 4
1002        y(i) = y(i) + alpha* x(i)
1003        y(i + 1) = y(i + 1) + alpha* x(i + 1)
1004        y(i + 2) = y(i + 2) + alpha* x(i + 2)
1005        y(i + 3) = y(i + 3) + alpha* x(i + 3)
1006      ENDDO
1007     
1008      END SUBROUTINE waxpy
1009 
1010SUBROUTINE rosenbrock(n, y, tstart, tend, &
1011           abstol, reltol,             &
1012           rcntrl, icntrl, rstatus, istatus, ierr)
1013!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1014!
1015!    Solves the system y'=F(t,y) using a Rosenbrock method defined by:
1016!
1017!     G = 1/(H*gamma(1)) - Jac(t0,Y0)
1018!     T_i = t0 + Alpha(i)*H
1019!     Y_i = Y0 + \sum_{j=1}^{i-1} A(i,j)*K_j
1020!     G *K_i = Fun( T_i,Y_i)+ \sum_{j=1}^S C(i,j)/H *K_j +
1021!         gamma(i)*dF/dT(t0,Y0)
1022!     Y1 = Y0 + \sum_{j=1}^S M(j)*K_j
1023!
1024!    For details on Rosenbrock methods and their implementation consult:
1025!      E. Hairer and G. Wanner
1026!      "Solving ODEs II. Stiff and differential-algebraic problems".
1027!      Springer series in computational mathematics,Springer-Verlag,1996.
1028!    The codes contained in the book inspired this implementation.
1029!
1030!    (C)  Adrian Sandu,August 2004
1031!    Virginia Polytechnic Institute and State University
1032!    Contact: sandu@cs.vt.edu
1033!    Revised by Philipp Miehe and Adrian Sandu,May 2006                 
1034!    This implementation is part of KPP - the Kinetic PreProcessor
1035!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1036!
1037!~~~>   input arguments:
1038!
1039!-     y(n)  = vector of initial conditions (at t=tstart)
1040!-    [tstart, tend]  = time range of integration
1041!     (if Tstart>Tend the integration is performed backwards in time)
1042!-    reltol, abstol = user precribed accuracy
1043!- SUBROUTINE  fun( t, y, ydot) = ode FUNCTION,
1044!                       returns Ydot = Y' = F(T,Y)
1045!- SUBROUTINE  jac( t, y, jcb) = jacobian of the ode FUNCTION,
1046!                       returns Jcb = dFun/dY
1047!-    icntrl(1:20)  = INTEGER inputs PARAMETERs
1048!-    rcntrl(1:20)  = REAL inputs PARAMETERs
1049!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1050!
1051!~~~>     output arguments:
1052!
1053!-    y(n)  - > vector of final states (at t- >tend)
1054!-    istatus(1:20) - > INTEGER output PARAMETERs
1055!-    rstatus(1:20) - > REAL output PARAMETERs
1056!-    ierr            - > job status upon RETURN
1057!                        success (positive value) or
1058!                        failure (negative value)
1059!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1060!
1061!~~~>     input PARAMETERs:
1062!
1063!    Note: For input parameters equal to zero the default values of the
1064!       corresponding variables are used.
1065!
1066!    ICNTRL(1) = 1: F = F(y)   Independent of T (AUTONOMOUS)
1067!              = 0: F = F(t,y) Depends on T (NON-AUTONOMOUS)
1068!
1069!    ICNTRL(2) = 0: AbsTol,RelTol are N-dimensional vectors
1070!              = 1: AbsTol,RelTol are scalars
1071!
1072!    ICNTRL(3)  -> selection of a particular Rosenbrock method
1073!        = 0 :    Rodas3 (default)
1074!        = 1 :    Ros2
1075!        = 2 :    Ros3
1076!        = 3 :    Ros4
1077!        = 4 :    Rodas3
1078!        = 5 :    Rodas4
1079!
1080!    ICNTRL(4)  -> maximum number of integration steps
1081!        For ICNTRL(4) =0) the default value of 100000 is used
1082!
1083!    RCNTRL(1)  -> Hmin,lower bound for the integration step size
1084!          It is strongly recommended to keep Hmin = ZERO
1085!    RCNTRL(2)  -> Hmax,upper bound for the integration step size
1086!    RCNTRL(3)  -> Hstart,starting value for the integration step size
1087!
1088!    RCNTRL(4)  -> FacMin,lower bound on step decrease factor (default=0.2)
1089!    RCNTRL(5)  -> FacMax,upper bound on step increase factor (default=6)
1090!    RCNTRL(6)  -> FacRej,step decrease factor after multiple rejections
1091!                          (default=0.1)
1092!    RCNTRL(7)  -> FacSafe,by which the new step is slightly smaller
1093!         than the predicted value  (default=0.9)
1094!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1095!
1096!
1097!    OUTPUT ARGUMENTS:
1098!    -----------------
1099!
1100!    T           -> T value for which the solution has been computed
1101!                     (after successful return T=Tend).
1102!
1103!    Y(N)        -> Numerical solution at T
1104!
1105!    IDID        -> Reports on successfulness upon return:
1106!                    = 1 for success
1107!                    < 0 for error (value equals error code)
1108!
1109!    ISTATUS(1)  -> No. of function calls
1110!    ISTATUS(2)  -> No. of jacobian calls
1111!    ISTATUS(3)  -> No. of steps
1112!    ISTATUS(4)  -> No. of accepted steps
1113!    ISTATUS(5)  -> No. of rejected steps (except at very beginning)
1114!    ISTATUS(6)  -> No. of LU decompositions
1115!    ISTATUS(7)  -> No. of forward/backward substitutions
1116!    ISTATUS(8)  -> No. of singular matrix decompositions
1117!
1118!    RSTATUS(1)  -> Texit,the time corresponding to the
1119!                     computed Y upon return
1120!    RSTATUS(2)  -> Hexit,last accepted step before exit
1121!    RSTATUS(3)  -> Hnew,last predicted step (not yet taken)
1122!                   For multiple restarts,use Hnew as Hstart
1123!                     in the subsequent run
1124!
1125!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1126
1127
1128!~~~>  arguments
1129   INTEGER,      INTENT(IN)  :: n
1130   REAL(kind=dp), INTENT(INOUT):: y(n)
1131   REAL(kind=dp), INTENT(IN)  :: tstart, tend
1132   REAL(kind=dp), INTENT(IN)  :: abstol(n), reltol(n)
1133   INTEGER,      INTENT(IN)  :: icntrl(20)
1134   REAL(kind=dp), INTENT(IN)  :: rcntrl(20)
1135   INTEGER,      INTENT(INOUT):: istatus(20)
1136   REAL(kind=dp), INTENT(INOUT):: rstatus(20)
1137   INTEGER, INTENT(OUT) :: ierr
1138!~~~>  PARAMETERs of the rosenbrock method, up to 6 stages
1139   INTEGER ::  ros_s, rosmethod
1140   INTEGER, PARAMETER :: rs2=1, rs3=2, rs4=3, rd3=4, rd4=5, rg3=6
1141   REAL(kind=dp):: ros_a(15), ros_c(15), ros_m(6), ros_e(6), &
1142                    ros_alpha(6), ros_gamma(6), ros_elo
1143   LOGICAL :: ros_newf(6)
1144   CHARACTER(len=12):: ros_name
1145!~~~>  local variables
1146   REAL(kind=dp):: roundoff, facmin, facmax, facrej, facsafe
1147   REAL(kind=dp):: hmin, hmax, hstart
1148   REAL(kind=dp):: texit
1149   INTEGER       :: i, uplimtol, max_no_steps
1150   LOGICAL       :: autonomous, vectortol
1151!~~~>   PARAMETERs
1152   REAL(kind=dp), PARAMETER :: zero = 0.0_dp, one  = 1.0_dp
1153   REAL(kind=dp), PARAMETER :: deltamin = 1.0e-5_dp
1154
1155!~~~>  initialize statistics
1156   istatus(1:8) = 0
1157   rstatus(1:3) = zero
1158
1159!~~~>  autonomous or time dependent ode. default is time dependent.
1160   autonomous = .not.(icntrl(1) == 0)
1161
1162!~~~>  for scalar tolerances (icntrl(2).ne.0) the code uses abstol(1)and reltol(1)
1163!   For Vector tolerances (ICNTRL(2) == 0) the code uses AbsTol(1:N) and RelTol(1:N)
1164   IF (icntrl(2) == 0)THEN
1165      vectortol = .TRUE.
1166      uplimtol  = n
1167   ELSE
1168      vectortol = .FALSE.
1169      uplimtol  = 1
1170   ENDIF
1171
1172!~~~>   initialize the particular rosenbrock method selected
1173   select CASE (icntrl(3))
1174     CASE (1)
1175       CALL ros2
1176     CASE (2)
1177       CALL ros3
1178     CASE (3)
1179       CALL ros4
1180     CASE (0, 4)
1181       CALL rodas3
1182     CASE (5)
1183       CALL rodas4
1184     CASE (6)
1185       CALL rang3
1186     CASE default
1187       PRINT *,'Unknown Rosenbrock method: ICNTRL(3) =',ICNTRL(3) 
1188       CALL ros_errormsg(- 2, tstart, zero, ierr)
1189       RETURN
1190   END select
1191
1192!~~~>   the maximum number of steps admitted
1193   IF (icntrl(4) == 0)THEN
1194      max_no_steps = 200000
1195   ELSEIF (icntrl(4)> 0)THEN
1196      max_no_steps=icntrl(4)
1197   ELSE
1198      PRINT *,'User-selected max no. of steps: ICNTRL(4) =',ICNTRL(4)
1199      CALL ros_errormsg(- 1, tstart, zero, ierr)
1200      RETURN
1201   ENDIF
1202
1203!~~~>  unit roundoff (1+ roundoff>1)
1204   roundoff = epsilon(one)
1205
1206!~~~>  lower bound on the step size: (positive value)
1207   IF (rcntrl(1) == zero)THEN
1208      hmin = zero
1209   ELSEIF (rcntrl(1)> zero)THEN
1210      hmin = rcntrl(1)
1211   ELSE
1212      PRINT *,'User-selected Hmin: RCNTRL(1) =',RCNTRL(1)
1213      CALL ros_errormsg(- 3, tstart, zero, ierr)
1214      RETURN
1215   ENDIF
1216!~~~>  upper bound on the step size: (positive value)
1217   IF (rcntrl(2) == zero)THEN
1218      hmax = abs(tend-tstart)
1219   ELSEIF (rcntrl(2)> zero)THEN
1220      hmax = min(abs(rcntrl(2)), abs(tend-tstart))
1221   ELSE
1222      PRINT *,'User-selected Hmax: RCNTRL(2) =',RCNTRL(2)
1223      CALL ros_errormsg(- 3, tstart, zero, ierr)
1224      RETURN
1225   ENDIF
1226!~~~>  starting step size: (positive value)
1227   IF (rcntrl(3) == zero)THEN
1228      hstart = max(hmin, deltamin)
1229   ELSEIF (rcntrl(3)> zero)THEN
1230      hstart = min(abs(rcntrl(3)), abs(tend-tstart))
1231   ELSE
1232      PRINT *,'User-selected Hstart: RCNTRL(3) =',RCNTRL(3)
1233      CALL ros_errormsg(- 3, tstart, zero, ierr)
1234      RETURN
1235   ENDIF
1236!~~~>  step size can be changed s.t.  facmin < hnew/hold < facmax
1237   IF (rcntrl(4) == zero)THEN
1238      facmin = 0.2_dp
1239   ELSEIF (rcntrl(4)> zero)THEN
1240      facmin = rcntrl(4)
1241   ELSE
1242      PRINT *,'User-selected FacMin: RCNTRL(4) =',RCNTRL(4)
1243      CALL ros_errormsg(- 4, tstart, zero, ierr)
1244      RETURN
1245   ENDIF
1246   IF (rcntrl(5) == zero)THEN
1247      facmax = 6.0_dp
1248   ELSEIF (rcntrl(5)> zero)THEN
1249      facmax = rcntrl(5)
1250   ELSE
1251      PRINT *,'User-selected FacMax: RCNTRL(5) =',RCNTRL(5)
1252      CALL ros_errormsg(- 4, tstart, zero, ierr)
1253      RETURN
1254   ENDIF
1255!~~~>   facrej: factor to decrease step after 2 succesive rejections
1256   IF (rcntrl(6) == zero)THEN
1257      facrej = 0.1_dp
1258   ELSEIF (rcntrl(6)> zero)THEN
1259      facrej = rcntrl(6)
1260   ELSE
1261      PRINT *,'User-selected FacRej: RCNTRL(6) =',RCNTRL(6)
1262      CALL ros_errormsg(- 4, tstart, zero, ierr)
1263      RETURN
1264   ENDIF
1265!~~~>   facsafe: safety factor in the computation of new step size
1266   IF (rcntrl(7) == zero)THEN
1267      facsafe = 0.9_dp
1268   ELSEIF (rcntrl(7)> zero)THEN
1269      facsafe = rcntrl(7)
1270   ELSE
1271      PRINT *,'User-selected FacSafe: RCNTRL(7) =',RCNTRL(7)
1272      CALL ros_errormsg(- 4, tstart, zero, ierr)
1273      RETURN
1274   ENDIF
1275!~~~>  check IF tolerances are reasonable
1276    DO i=1, uplimtol
1277      IF ((abstol(i)<= zero).or. (reltol(i)<= 10.0_dp* roundoff)&
1278         .or. (reltol(i)>= 1.0_dp))THEN
1279        PRINT *,' AbsTol(',i,') = ',AbsTol(i)
1280        PRINT *,' RelTol(',i,') = ',RelTol(i)
1281        CALL ros_errormsg(- 5, tstart, zero, ierr)
1282        RETURN
1283      ENDIF
1284    ENDDO
1285
1286
1287!~~~>  CALL rosenbrock method
1288   CALL ros_integrator(y, tstart, tend, texit,  &
1289        abstol, reltol,                         &
1290!  Integration parameters
1291        autonomous, vectortol, max_no_steps,    &
1292        roundoff, hmin, hmax, hstart,           &
1293        facmin, facmax, facrej, facsafe,        &
1294!  Error indicator
1295        ierr)
1296
1297!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1298CONTAINS !  SUBROUTINEs internal to rosenbrock
1299!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1300   
1301!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
1302 SUBROUTINE ros_errormsg(code, t, h, ierr)
1303!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
1304!    Handles all error messages
1305!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
1306 
1307   REAL(kind=dp), INTENT(IN):: t, h
1308   INTEGER, INTENT(IN) :: code
1309   INTEGER, INTENT(OUT):: ierr
1310   
1311   ierr = code
1312   print * , &
1313     'Forced exit from Rosenbrock due to the following error:' 
1314     
1315   select CASE (code)
1316    CASE (- 1) 
1317      PRINT *,'--> Improper value for maximal no of steps'
1318    CASE (- 2) 
1319      PRINT *,'--> Selected Rosenbrock method not implemented'
1320    CASE (- 3) 
1321      PRINT *,'--> Hmin/Hmax/Hstart must be positive'
1322    CASE (- 4) 
1323      PRINT *,'--> FacMin/FacMax/FacRej must be positive'
1324    CASE (- 5)
1325      PRINT *,'--> Improper tolerance values'
1326    CASE (- 6)
1327      PRINT *,'--> No of steps exceeds maximum bound'
1328    CASE (- 7)
1329      PRINT *,'--> Step size too small: T + 10*H = T',&
1330            ' or H < Roundoff'
1331    CASE (- 8) 
1332      PRINT *,'--> Matrix is repeatedly singular'
1333    CASE default
1334      PRINT *,'Unknown Error code: ',Code
1335   END select
1336   
1337   print * , "t=", t, "and h=", h
1338     
1339 END SUBROUTINE ros_errormsg
1340
1341!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1342 SUBROUTINE ros_integrator (y, tstart, tend, t, &
1343        abstol, reltol,                         &
1344!~~~> integration PARAMETERs
1345        autonomous, vectortol, max_no_steps,    &
1346        roundoff, hmin, hmax, hstart,           &
1347        facmin, facmax, facrej, facsafe,        &
1348!~~~> error indicator
1349        ierr)
1350!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1351!   Template for the implementation of a generic Rosenbrock method
1352!      defined by ros_S (no of stages)
1353!      and its coefficients ros_{A,C,M,E,Alpha,Gamma}
1354!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1355
1356
1357!~~~> input: the initial condition at tstart; output: the solution at t
1358   REAL(kind=dp), INTENT(INOUT):: y(n)
1359!~~~> input: integration interval
1360   REAL(kind=dp), INTENT(IN):: tstart, tend
1361!~~~> output: time at which the solution is RETURNed (t=tendIF success)
1362   REAL(kind=dp), INTENT(OUT)::  t
1363!~~~> input: tolerances
1364   REAL(kind=dp), INTENT(IN)::  abstol(n), reltol(n)
1365!~~~> input: integration PARAMETERs
1366   LOGICAL, INTENT(IN):: autonomous, vectortol
1367   REAL(kind=dp), INTENT(IN):: hstart, hmin, hmax
1368   INTEGER, INTENT(IN):: max_no_steps
1369   REAL(kind=dp), INTENT(IN):: roundoff, facmin, facmax, facrej, facsafe
1370!~~~> output: error indicator
1371   INTEGER, INTENT(OUT):: ierr
1372! ~~~~ Local variables
1373   REAL(kind=dp):: ynew(n), fcn0(n), fcn(n)
1374   REAL(kind=dp):: k(n* ros_s), dfdt(n)
1375#ifdef full_algebra   
1376   REAL(kind=dp):: jac0(n, n), ghimj(n, n)
1377#else
1378   REAL(kind=dp):: jac0(lu_nonzero), ghimj(lu_nonzero)
1379#endif
1380   REAL(kind=dp):: h, hnew, hc, hg, fac, tau
1381   REAL(kind=dp):: err, yerr(n)
1382   INTEGER :: pivot(n), direction, ioffset, j, istage
1383   LOGICAL :: rejectlasth, rejectmoreh, singular
1384!~~~>  local PARAMETERs
1385   REAL(kind=dp), PARAMETER :: zero = 0.0_dp, one  = 1.0_dp
1386   REAL(kind=dp), PARAMETER :: deltamin = 1.0e-5_dp
1387!~~~>  locally called FUNCTIONs
1388!    REAL(kind=dp) WLAMCH
1389!    EXTERNAL WLAMCH
1390!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1391
1392
1393!~~~>  initial preparations
1394   t = tstart
1395   rstatus(nhexit) = zero
1396   h = min( max(abs(hmin), abs(hstart)), abs(hmax))
1397   IF (abs(h)<= 10.0_dp* roundoff)h = deltamin
1398
1399   IF (tend  >=  tstart)THEN
1400     direction = + 1
1401   ELSE
1402     direction = - 1
1403   ENDIF
1404   h = direction* h
1405
1406   rejectlasth=.FALSE.
1407   rejectmoreh=.FALSE.
1408
1409!~~~> time loop begins below
1410
1411timeloop: DO WHILE((direction > 0).and.((t- tend) + roundoff <= zero)&
1412       .or. (direction < 0).and.((tend-t) + roundoff <= zero))
1413
1414   IF (istatus(nstp)> max_no_steps)THEN  ! too many steps
1415      CALL ros_errormsg(- 6, t, h, ierr)
1416      RETURN
1417   ENDIF
1418   IF (((t+ 0.1_dp* h) == t).or.(h <= roundoff))THEN  ! step size too small
1419      CALL ros_errormsg(- 7, t, h, ierr)
1420      RETURN
1421   ENDIF
1422
1423!~~~>  limit h IF necessary to avoid going beyond tend
1424   h = min(h, abs(tend-t))
1425
1426!~~~>   compute the FUNCTION at current time
1427   CALL funtemplate(t, y, fcn0)
1428   istatus(nfun) = istatus(nfun) + 1
1429
1430!~~~>  compute the FUNCTION derivative with respect to t
1431   IF (.not.autonomous)THEN
1432      CALL ros_funtimederivative(t, roundoff, y, &
1433                fcn0, dfdt)
1434   ENDIF
1435
1436!~~~>   compute the jacobian at current time
1437   CALL jactemplate(t, y, jac0)
1438   istatus(njac) = istatus(njac) + 1
1439
1440!~~~>  repeat step calculation until current step accepted
1441untilaccepted: do
1442
1443   CALL ros_preparematrix(h, direction, ros_gamma(1), &
1444          jac0, ghimj, pivot, singular)
1445   IF (singular)THEN ! more than 5 consecutive failed decompositions
1446       CALL ros_errormsg(- 8, t, h, ierr)
1447       RETURN
1448   ENDIF
1449
1450!~~~>   compute the stages
1451stage: DO istage = 1, ros_s
1452
1453      ! current istage offset. current istage vector is k(ioffset+ 1:ioffset+ n)
1454       ioffset = n* (istage-1)
1455
1456      ! for the 1st istage the FUNCTION has been computed previously
1457       IF (istage == 1)THEN
1458         !slim: CALL wcopy(n, fcn0, 1, fcn, 1)
1459       fcn(1:n) = fcn0(1:n)
1460      ! istage>1 and a new FUNCTION evaluation is needed at the current istage
1461       ELSEIF(ros_newf(istage))THEN
1462         !slim: CALL wcopy(n, y, 1, ynew, 1)
1463       ynew(1:n) = y(1:n)
1464         DO j = 1, istage-1
1465           CALL waxpy(n, ros_a((istage-1) * (istage-2) /2+ j), &
1466            k(n* (j- 1) + 1), 1, ynew, 1)
1467         ENDDO
1468         tau = t + ros_alpha(istage) * direction* h
1469         CALL funtemplate(tau, ynew, fcn)
1470         istatus(nfun) = istatus(nfun) + 1
1471       ENDIF ! IF istage == 1 ELSEIF ros_newf(istage)
1472       !slim: CALL wcopy(n, fcn, 1, k(ioffset+ 1), 1)
1473       k(ioffset+ 1:ioffset+ n) = fcn(1:n)
1474       DO j = 1, istage-1
1475         hc = ros_c((istage-1) * (istage-2) /2+ j) /(direction* h)
1476         CALL waxpy(n, hc, k(n* (j- 1) + 1), 1, k(ioffset+ 1), 1)
1477       ENDDO
1478       IF ((.not. autonomous).and.(ros_gamma(istage).ne.zero))THEN
1479         hg = direction* h* ros_gamma(istage)
1480         CALL waxpy(n, hg, dfdt, 1, k(ioffset+ 1), 1)
1481       ENDIF
1482       CALL ros_solve(ghimj, pivot, k(ioffset+ 1))
1483
1484   END DO stage
1485
1486
1487!~~~>  compute the new solution
1488   !slim: CALL wcopy(n, y, 1, ynew, 1)
1489   ynew(1:n) = y(1:n)
1490   DO j=1, ros_s
1491         CALL waxpy(n, ros_m(j), k(n* (j- 1) + 1), 1, ynew, 1)
1492   ENDDO
1493
1494!~~~>  compute the error estimation
1495   !slim: CALL wscal(n, zero, yerr, 1)
1496   yerr(1:n) = zero
1497   DO j=1, ros_s
1498        CALL waxpy(n, ros_e(j), k(n* (j- 1) + 1), 1, yerr, 1)
1499   ENDDO
1500   err = ros_errornorm(y, ynew, yerr, abstol, reltol, vectortol)
1501
1502!~~~> new step size is bounded by facmin <= hnew/h <= facmax
1503   fac  = min(facmax, max(facmin, facsafe/err** (one/ros_elo)))
1504   hnew = h* fac
1505
1506!~~~>  check the error magnitude and adjust step size
1507   istatus(nstp) = istatus(nstp) + 1
1508   IF ((err <= one).or.(h <= hmin))THEN  !~~~> accept step
1509      istatus(nacc) = istatus(nacc) + 1
1510      !slim: CALL wcopy(n, ynew, 1, y, 1)
1511      y(1:n) = ynew(1:n)
1512      t = t + direction* h
1513      hnew = max(hmin, min(hnew, hmax))
1514      IF (rejectlasth)THEN  ! no step size increase after a rejected step
1515         hnew = min(hnew, h)
1516      ENDIF
1517      rstatus(nhexit) = h
1518      rstatus(nhnew) = hnew
1519      rstatus(ntexit) = t
1520      rejectlasth = .FALSE.
1521      rejectmoreh = .FALSE.
1522      h = hnew
1523      exit untilaccepted ! exit the loop: WHILE step not accepted
1524   ELSE           !~~~> reject step
1525      IF (rejectmoreh)THEN
1526         hnew = h* facrej
1527      ENDIF
1528      rejectmoreh = rejectlasth
1529      rejectlasth = .TRUE.
1530      h = hnew
1531      IF (istatus(nacc)>= 1) istatus(nrej) = istatus(nrej) + 1
1532   ENDIF ! err <= 1
1533
1534   END DO untilaccepted
1535
1536   END DO timeloop
1537
1538!~~~> succesful exit
1539   ierr = 1  !~~~> the integration was successful
1540
1541  END SUBROUTINE ros_integrator
1542
1543
1544!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1545  REAL(kind=dp)FUNCTION ros_errornorm(y, ynew, yerr, &
1546                               abstol, reltol, vectortol)
1547!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1548!~~~> computes the "scaled norm" of the error vector yerr
1549!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1550
1551! Input arguments
1552   REAL(kind=dp), INTENT(IN):: y(n), ynew(n), &
1553          yerr(n), abstol(n), reltol(n)
1554   LOGICAL, INTENT(IN)::  vectortol
1555! Local variables
1556   REAL(kind=dp):: err, scale, ymax
1557   INTEGER  :: i
1558   REAL(kind=dp), PARAMETER :: zero = 0.0_dp
1559
1560   err = zero
1561   DO i=1, n
1562     ymax = max(abs(y(i)), abs(ynew(i)))
1563     IF (vectortol)THEN
1564       scale = abstol(i) + reltol(i) * ymax
1565     ELSE
1566       scale = abstol(1) + reltol(1) * ymax
1567     ENDIF
1568     err = err+ (yerr(i) /scale) ** 2
1569   ENDDO
1570   err  = sqrt(err/n)
1571
1572   ros_errornorm = max(err, 1.0d-10)
1573
1574  END FUNCTION ros_errornorm
1575
1576
1577!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1578  SUBROUTINE ros_funtimederivative(t, roundoff, y, &
1579                fcn0, dfdt)
1580!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1581!~~~> the time partial derivative of the FUNCTION by finite differences
1582!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1583
1584!~~~> input arguments
1585   REAL(kind=dp), INTENT(IN):: t, roundoff, y(n), fcn0(n)
1586!~~~> output arguments
1587   REAL(kind=dp), INTENT(OUT):: dfdt(n)
1588!~~~> local variables
1589   REAL(kind=dp):: delta
1590   REAL(kind=dp), PARAMETER :: one = 1.0_dp, deltamin = 1.0e-6_dp
1591
1592   delta = sqrt(roundoff) * max(deltamin, abs(t))
1593   CALL funtemplate(t+ delta, y, dfdt)
1594   istatus(nfun) = istatus(nfun) + 1
1595   CALL waxpy(n, (- one), fcn0, 1, dfdt, 1)
1596   CALL wscal(n, (one/delta), dfdt, 1)
1597
1598  END SUBROUTINE ros_funtimederivative
1599
1600
1601!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1602  SUBROUTINE ros_preparematrix(h, direction, gam, &
1603             jac0, ghimj, pivot, singular)
1604! --- --- --- --- --- --- --- --- --- --- --- --- ---
1605!  Prepares the LHS matrix for stage calculations
1606!  1.  Construct Ghimj = 1/(H*ham) - Jac0
1607!      "(Gamma H) Inverse Minus Jacobian"
1608!  2.  Repeat LU decomposition of Ghimj until successful.
1609!       -half the step size if LU decomposition fails and retry
1610!       -exit after 5 consecutive fails
1611! --- --- --- --- --- --- --- --- --- --- --- --- ---
1612
1613!~~~> input arguments
1614#ifdef full_algebra   
1615   REAL(kind=dp), INTENT(IN)::  jac0(n, n)
1616#else
1617   REAL(kind=dp), INTENT(IN)::  jac0(lu_nonzero)
1618#endif   
1619   REAL(kind=dp), INTENT(IN)::  gam
1620   INTEGER, INTENT(IN)::  direction
1621!~~~> output arguments
1622#ifdef full_algebra   
1623   REAL(kind=dp), INTENT(OUT):: ghimj(n, n)
1624#else
1625   REAL(kind=dp), INTENT(OUT):: ghimj(lu_nonzero)
1626#endif   
1627   LOGICAL, INTENT(OUT)::  singular
1628   INTEGER, INTENT(OUT)::  pivot(n)
1629!~~~> inout arguments
1630   REAL(kind=dp), INTENT(INOUT):: h   ! step size is decreased when lu fails
1631!~~~> local variables
1632   INTEGER  :: i, ising, nconsecutive
1633   REAL(kind=dp):: ghinv
1634   REAL(kind=dp), PARAMETER :: one  = 1.0_dp, half = 0.5_dp
1635
1636   nconsecutive = 0
1637   singular = .TRUE.
1638
1639   DO WHILE (singular)
1640
1641!~~~>    construct ghimj = 1/(h* gam) - jac0
1642#ifdef full_algebra   
1643     !slim: CALL wcopy(n* n, jac0, 1, ghimj, 1)
1644     !slim: CALL wscal(n* n, (- one), ghimj, 1)
1645     ghimj = - jac0
1646     ghinv = one/(direction* h* gam)
1647     DO i=1, n
1648       ghimj(i, i) = ghimj(i, i) + ghinv
1649     ENDDO
1650#else
1651     !slim: CALL wcopy(lu_nonzero, jac0, 1, ghimj, 1)
1652     !slim: CALL wscal(lu_nonzero, (- one), ghimj, 1)
1653     ghimj(1:lu_nonzero) = - jac0(1:lu_nonzero)
1654     ghinv = one/(direction* h* gam)
1655     DO i=1, n
1656       ghimj(lu_diag(i)) = ghimj(lu_diag(i)) + ghinv
1657     ENDDO
1658#endif   
1659!~~~>    compute lu decomposition
1660     CALL ros_decomp( ghimj, pivot, ising)
1661     IF (ising == 0)THEN
1662!~~~>    IF successful done
1663        singular = .FALSE.
1664     ELSE ! ising .ne. 0
1665!~~~>    IF unsuccessful half the step size; IF 5 consecutive fails THEN RETURN
1666        istatus(nsng) = istatus(nsng) + 1
1667        nconsecutive = nconsecutive+1
1668        singular = .TRUE.
1669        PRINT*,'Warning: LU Decomposition returned ISING = ',ISING
1670        IF (nconsecutive <= 5)THEN ! less than 5 consecutive failed decompositions
1671           h = h* half
1672        ELSE  ! more than 5 consecutive failed decompositions
1673           RETURN
1674        ENDIF  ! nconsecutive
1675      ENDIF    ! ising
1676
1677   END DO ! WHILE singular
1678
1679  END SUBROUTINE ros_preparematrix
1680
1681
1682!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1683  SUBROUTINE ros_decomp( a, pivot, ising)
1684!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1685!  Template for the LU decomposition
1686!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1687!~~~> inout variables
1688#ifdef full_algebra   
1689   REAL(kind=dp), INTENT(INOUT):: a(n, n)
1690#else   
1691   REAL(kind=dp), INTENT(INOUT):: a(lu_nonzero)
1692#endif
1693!~~~> output variables
1694   INTEGER, INTENT(OUT):: pivot(n), ising
1695
1696#ifdef full_algebra   
1697   CALL  dgetrf( n, n, a, n, pivot, ising)
1698#else   
1699   CALL kppdecomp(a, ising)
1700   pivot(1) = 1
1701#endif
1702   istatus(ndec) = istatus(ndec) + 1
1703
1704  END SUBROUTINE ros_decomp
1705
1706
1707!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1708  SUBROUTINE ros_solve( a, pivot, b)
1709!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1710!  Template for the forward/backward substitution (using pre-computed LU decomposition)
1711!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1712!~~~> input variables
1713#ifdef full_algebra   
1714   REAL(kind=dp), INTENT(IN):: a(n, n)
1715   INTEGER :: ising
1716#else   
1717   REAL(kind=dp), INTENT(IN):: a(lu_nonzero)
1718#endif
1719   INTEGER, INTENT(IN):: pivot(n)
1720!~~~> inout variables
1721   REAL(kind=dp), INTENT(INOUT):: b(n)
1722
1723#ifdef full_algebra   
1724   CALL  DGETRS( 'N',N ,1,A,N,Pivot,b,N,ISING)
1725   IF (info < 0)THEN
1726      print* , "error in dgetrs. ising=", ising
1727   ENDIF 
1728#else   
1729   CALL kppsolve( a, b)
1730#endif
1731
1732   istatus(nsol) = istatus(nsol) + 1
1733
1734  END SUBROUTINE ros_solve
1735
1736
1737
1738!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1739  SUBROUTINE ros2
1740!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1741! --- AN L-STABLE METHOD,2 stages,order 2
1742!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1743
1744   double precision g
1745
1746    g = 1.0_dp + 1.0_dp/sqrt(2.0_dp)
1747    rosmethod = rs2
1748!~~~> name of the method
1749    ros_Name = 'ROS-2'
1750!~~~> number of stages
1751    ros_s = 2
1752
1753!~~~> the coefficient matrices a and c are strictly lower triangular.
1754!   The lower triangular (subdiagonal) elements are stored in row-wise order:
1755!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
1756!   The general mapping formula is:
1757!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
1758!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
1759
1760    ros_a(1) = (1.0_dp) /g
1761    ros_c(1) = (- 2.0_dp) /g
1762!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
1763!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
1764    ros_newf(1) = .TRUE.
1765    ros_newf(2) = .TRUE.
1766!~~~> m_i = coefficients for new step solution
1767    ros_m(1) = (3.0_dp) /(2.0_dp* g)
1768    ros_m(2) = (1.0_dp) /(2.0_dp* g)
1769! E_i = Coefficients for error estimator
1770    ros_e(1) = 1.0_dp/(2.0_dp* g)
1771    ros_e(2) = 1.0_dp/(2.0_dp* g)
1772!~~~> ros_elo = estimator of local order - the minimum between the
1773!    main and the embedded scheme orders plus one
1774    ros_elo = 2.0_dp
1775!~~~> y_stage_i ~ y( t + h* alpha_i)
1776    ros_alpha(1) = 0.0_dp
1777    ros_alpha(2) = 1.0_dp
1778!~~~> gamma_i = \sum_j  gamma_{i, j}
1779    ros_gamma(1) = g
1780    ros_gamma(2) = -g
1781
1782 END SUBROUTINE ros2
1783
1784
1785!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1786  SUBROUTINE ros3
1787!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1788! --- AN L-STABLE METHOD,3 stages,order 3,2 function evaluations
1789!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1790
1791   rosmethod = rs3
1792!~~~> name of the method
1793   ros_Name = 'ROS-3'
1794!~~~> number of stages
1795   ros_s = 3
1796
1797!~~~> the coefficient matrices a and c are strictly lower triangular.
1798!   The lower triangular (subdiagonal) elements are stored in row-wise order:
1799!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
1800!   The general mapping formula is:
1801!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
1802!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
1803
1804   ros_a(1) = 1.0_dp
1805   ros_a(2) = 1.0_dp
1806   ros_a(3) = 0.0_dp
1807
1808   ros_c(1) = - 0.10156171083877702091975600115545e+01_dp
1809   ros_c(2) =  0.40759956452537699824805835358067e+01_dp
1810   ros_c(3) =  0.92076794298330791242156818474003e+01_dp
1811!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
1812!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
1813   ros_newf(1) = .TRUE.
1814   ros_newf(2) = .TRUE.
1815   ros_newf(3) = .FALSE.
1816!~~~> m_i = coefficients for new step solution
1817   ros_m(1) =  0.1e+01_dp
1818   ros_m(2) =  0.61697947043828245592553615689730e+01_dp
1819   ros_m(3) = - 0.42772256543218573326238373806514_dp
1820! E_i = Coefficients for error estimator
1821   ros_e(1) =  0.5_dp
1822   ros_e(2) = - 0.29079558716805469821718236208017e+01_dp
1823   ros_e(3) =  0.22354069897811569627360909276199_dp
1824!~~~> ros_elo = estimator of local order - the minimum between the
1825!    main and the embedded scheme orders plus 1
1826   ros_elo = 3.0_dp
1827!~~~> y_stage_i ~ y( t + h* alpha_i)
1828   ros_alpha(1) = 0.0_dp
1829   ros_alpha(2) = 0.43586652150845899941601945119356_dp
1830   ros_alpha(3) = 0.43586652150845899941601945119356_dp
1831!~~~> gamma_i = \sum_j  gamma_{i, j}
1832   ros_gamma(1) = 0.43586652150845899941601945119356_dp
1833   ros_gamma(2) = 0.24291996454816804366592249683314_dp
1834   ros_gamma(3) = 0.21851380027664058511513169485832e+01_dp
1835
1836  END SUBROUTINE ros3
1837
1838!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1839
1840
1841!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1842  SUBROUTINE ros4
1843!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1844!     L-STABLE ROSENBROCK METHOD OF ORDER 4,WITH 4 STAGES
1845!     L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3
1846!
1847!      E. HAIRER AND G. WANNER,SOLVING ORDINARY DIFFERENTIAL
1848!      EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
1849!      SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
1850!      SPRINGER-VERLAG (1990)
1851!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1852
1853
1854   rosmethod = rs4
1855!~~~> name of the method
1856   ros_Name = 'ROS-4'
1857!~~~> number of stages
1858   ros_s = 4
1859
1860!~~~> the coefficient matrices a and c are strictly lower triangular.
1861!   The lower triangular (subdiagonal) elements are stored in row-wise order:
1862!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
1863!   The general mapping formula is:
1864!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
1865!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
1866
1867   ros_a(1) = 0.2000000000000000e+01_dp
1868   ros_a(2) = 0.1867943637803922e+01_dp
1869   ros_a(3) = 0.2344449711399156_dp
1870   ros_a(4) = ros_a(2)
1871   ros_a(5) = ros_a(3)
1872   ros_a(6) = 0.0_dp
1873
1874   ros_c(1) = -0.7137615036412310e+01_dp
1875   ros_c(2) = 0.2580708087951457e+01_dp
1876   ros_c(3) = 0.6515950076447975_dp
1877   ros_c(4) = -0.2137148994382534e+01_dp
1878   ros_c(5) = -0.3214669691237626_dp
1879   ros_c(6) = -0.6949742501781779_dp
1880!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
1881!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
1882   ros_newf(1) = .TRUE.
1883   ros_newf(2) = .TRUE.
1884   ros_newf(3) = .TRUE.
1885   ros_newf(4) = .FALSE.
1886!~~~> m_i = coefficients for new step solution
1887   ros_m(1) = 0.2255570073418735e+01_dp
1888   ros_m(2) = 0.2870493262186792_dp
1889   ros_m(3) = 0.4353179431840180_dp
1890   ros_m(4) = 0.1093502252409163e+01_dp
1891!~~~> e_i  = coefficients for error estimator
1892   ros_e(1) = -0.2815431932141155_dp
1893   ros_e(2) = -0.7276199124938920e-01_dp
1894   ros_e(3) = -0.1082196201495311_dp
1895   ros_e(4) = -0.1093502252409163e+01_dp
1896!~~~> ros_elo  = estimator of local order - the minimum between the
1897!    main and the embedded scheme orders plus 1
1898   ros_elo  = 4.0_dp
1899!~~~> y_stage_i ~ y( t + h* alpha_i)
1900   ros_alpha(1) = 0.0_dp
1901   ros_alpha(2) = 0.1145640000000000e+01_dp
1902   ros_alpha(3) = 0.6552168638155900_dp
1903   ros_alpha(4) = ros_alpha(3)
1904!~~~> gamma_i = \sum_j  gamma_{i, j}
1905   ros_gamma(1) = 0.5728200000000000_dp
1906   ros_gamma(2) = -0.1769193891319233e+01_dp
1907   ros_gamma(3) = 0.7592633437920482_dp
1908   ros_gamma(4) = -0.1049021087100450_dp
1909
1910  END SUBROUTINE ros4
1911
1912!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1913  SUBROUTINE rodas3
1914!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1915! --- A STIFFLY-STABLE METHOD,4 stages,order 3
1916!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1917
1918
1919   rosmethod = rd3
1920!~~~> name of the method
1921   ros_Name = 'RODAS-3'
1922!~~~> number of stages
1923   ros_s = 4
1924
1925!~~~> the coefficient matrices a and c are strictly lower triangular.
1926!   The lower triangular (subdiagonal) elements are stored in row-wise order:
1927!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
1928!   The general mapping formula is:
1929!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
1930!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
1931
1932   ros_a(1) = 0.0_dp
1933   ros_a(2) = 2.0_dp
1934   ros_a(3) = 0.0_dp
1935   ros_a(4) = 2.0_dp
1936   ros_a(5) = 0.0_dp
1937   ros_a(6) = 1.0_dp
1938
1939   ros_c(1) = 4.0_dp
1940   ros_c(2) = 1.0_dp
1941   ros_c(3) = -1.0_dp
1942   ros_c(4) = 1.0_dp
1943   ros_c(5) = -1.0_dp
1944   ros_c(6) = -(8.0_dp/3.0_dp)
1945
1946!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
1947!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
1948   ros_newf(1) = .TRUE.
1949   ros_newf(2) = .FALSE.
1950   ros_newf(3) = .TRUE.
1951   ros_newf(4) = .TRUE.
1952!~~~> m_i = coefficients for new step solution
1953   ros_m(1) = 2.0_dp
1954   ros_m(2) = 0.0_dp
1955   ros_m(3) = 1.0_dp
1956   ros_m(4) = 1.0_dp
1957!~~~> e_i  = coefficients for error estimator
1958   ros_e(1) = 0.0_dp
1959   ros_e(2) = 0.0_dp
1960   ros_e(3) = 0.0_dp
1961   ros_e(4) = 1.0_dp
1962!~~~> ros_elo  = estimator of local order - the minimum between the
1963!    main and the embedded scheme orders plus 1
1964   ros_elo  = 3.0_dp
1965!~~~> y_stage_i ~ y( t + h* alpha_i)
1966   ros_alpha(1) = 0.0_dp
1967   ros_alpha(2) = 0.0_dp
1968   ros_alpha(3) = 1.0_dp
1969   ros_alpha(4) = 1.0_dp
1970!~~~> gamma_i = \sum_j  gamma_{i, j}
1971   ros_gamma(1) = 0.5_dp
1972   ros_gamma(2) = 1.5_dp
1973   ros_gamma(3) = 0.0_dp
1974   ros_gamma(4) = 0.0_dp
1975
1976  END SUBROUTINE rodas3
1977
1978!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1979  SUBROUTINE rodas4
1980!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1981!     STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4,WITH 6 STAGES
1982!
1983!      E. HAIRER AND G. WANNER,SOLVING ORDINARY DIFFERENTIAL
1984!      EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
1985!      SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
1986!      SPRINGER-VERLAG (1996)
1987!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1988
1989
1990    rosmethod = rd4
1991!~~~> name of the method
1992    ros_Name = 'RODAS-4'
1993!~~~> number of stages
1994    ros_s = 6
1995
1996!~~~> y_stage_i ~ y( t + h* alpha_i)
1997    ros_alpha(1) = 0.000_dp
1998    ros_alpha(2) = 0.386_dp
1999    ros_alpha(3) = 0.210_dp
2000    ros_alpha(4) = 0.630_dp
2001    ros_alpha(5) = 1.000_dp
2002    ros_alpha(6) = 1.000_dp
2003
2004!~~~> gamma_i = \sum_j  gamma_{i, j}
2005    ros_gamma(1) = 0.2500000000000000_dp
2006    ros_gamma(2) = -0.1043000000000000_dp
2007    ros_gamma(3) = 0.1035000000000000_dp
2008    ros_gamma(4) = -0.3620000000000023e-01_dp
2009    ros_gamma(5) = 0.0_dp
2010    ros_gamma(6) = 0.0_dp
2011
2012!~~~> the coefficient matrices a and c are strictly lower triangular.
2013!   The lower triangular (subdiagonal) elements are stored in row-wise order:
2014!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
2015!   The general mapping formula is:  A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
2016!                  C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
2017
2018    ros_a(1) = 0.1544000000000000e+01_dp
2019    ros_a(2) = 0.9466785280815826_dp
2020    ros_a(3) = 0.2557011698983284_dp
2021    ros_a(4) = 0.3314825187068521e+01_dp
2022    ros_a(5) = 0.2896124015972201e+01_dp
2023    ros_a(6) = 0.9986419139977817_dp
2024    ros_a(7) = 0.1221224509226641e+01_dp
2025    ros_a(8) = 0.6019134481288629e+01_dp
2026    ros_a(9) = 0.1253708332932087e+02_dp
2027    ros_a(10) = -0.6878860361058950_dp
2028    ros_a(11) = ros_a(7)
2029    ros_a(12) = ros_a(8)
2030    ros_a(13) = ros_a(9)
2031    ros_a(14) = ros_a(10)
2032    ros_a(15) = 1.0_dp
2033
2034    ros_c(1) = -0.5668800000000000e+01_dp
2035    ros_c(2) = -0.2430093356833875e+01_dp
2036    ros_c(3) = -0.2063599157091915_dp
2037    ros_c(4) = -0.1073529058151375_dp
2038    ros_c(5) = -0.9594562251023355e+01_dp
2039    ros_c(6) = -0.2047028614809616e+02_dp
2040    ros_c(7) = 0.7496443313967647e+01_dp
2041    ros_c(8) = -0.1024680431464352e+02_dp
2042    ros_c(9) = -0.3399990352819905e+02_dp
2043    ros_c(10) = 0.1170890893206160e+02_dp
2044    ros_c(11) = 0.8083246795921522e+01_dp
2045    ros_c(12) = -0.7981132988064893e+01_dp
2046    ros_c(13) = -0.3152159432874371e+02_dp
2047    ros_c(14) = 0.1631930543123136e+02_dp
2048    ros_c(15) = -0.6058818238834054e+01_dp
2049
2050!~~~> m_i = coefficients for new step solution
2051    ros_m(1) = ros_a(7)
2052    ros_m(2) = ros_a(8)
2053    ros_m(3) = ros_a(9)
2054    ros_m(4) = ros_a(10)
2055    ros_m(5) = 1.0_dp
2056    ros_m(6) = 1.0_dp
2057
2058!~~~> e_i  = coefficients for error estimator
2059    ros_e(1) = 0.0_dp
2060    ros_e(2) = 0.0_dp
2061    ros_e(3) = 0.0_dp
2062    ros_e(4) = 0.0_dp
2063    ros_e(5) = 0.0_dp
2064    ros_e(6) = 1.0_dp
2065
2066!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
2067!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
2068    ros_newf(1) = .TRUE.
2069    ros_newf(2) = .TRUE.
2070    ros_newf(3) = .TRUE.
2071    ros_newf(4) = .TRUE.
2072    ros_newf(5) = .TRUE.
2073    ros_newf(6) = .TRUE.
2074
2075!~~~> ros_elo  = estimator of local order - the minimum between the
2076!        main and the embedded scheme orders plus 1
2077    ros_elo = 4.0_dp
2078
2079  END SUBROUTINE rodas4
2080!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2081  SUBROUTINE rang3
2082!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2083! STIFFLY-STABLE W METHOD OF ORDER 3,WITH 4 STAGES
2084!
2085! J. RANG and L. ANGERMANN
2086! NEW ROSENBROCK W-METHODS OF ORDER 3
2087! FOR PARTIAL DIFFERENTIAL ALGEBRAIC
2088!        EQUATIONS OF INDEX 1                                             
2089! BIT Numerical Mathematics (2005) 45: 761-787
2090!  DOI: 10.1007/s10543-005-0035-y
2091! Table 4.1-4.2
2092!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2093
2094
2095    rosmethod = rg3
2096!~~~> name of the method
2097    ros_Name = 'RANG-3'
2098!~~~> number of stages
2099    ros_s = 4
2100
2101    ros_a(1) = 5.09052051067020d+00;
2102    ros_a(2) = 5.09052051067020d+00;
2103    ros_a(3) = 0.0d0;
2104    ros_a(4) = 4.97628111010787d+00;
2105    ros_a(5) = 2.77268164715849d-02;
2106    ros_a(6) = 2.29428036027904d-01;
2107
2108    ros_c(1) = - 1.16790812312283d+01;
2109    ros_c(2) = - 1.64057326467367d+01;
2110    ros_c(3) = - 2.77268164715850d-01;
2111    ros_c(4) = - 8.38103960500476d+00;
2112    ros_c(5) = - 8.48328409199343d-01;
2113    ros_c(6) =  2.87009860433106d-01;
2114
2115    ros_m(1) =  5.22582761233094d+00;
2116    ros_m(2) = - 5.56971148154165d-01;
2117    ros_m(3) =  3.57979469353645d-01;
2118    ros_m(4) =  1.72337398521064d+00;
2119
2120    ros_e(1) = - 5.16845212784040d+00;
2121    ros_e(2) = - 1.26351942603842d+00;
2122    ros_e(3) = - 1.11022302462516d-16;
2123    ros_e(4) =  2.22044604925031d-16;
2124
2125    ros_alpha(1) = 0.0d00;
2126    ros_alpha(2) = 2.21878746765329d+00;
2127    ros_alpha(3) = 2.21878746765329d+00;
2128    ros_alpha(4) = 1.55392337535788d+00;
2129
2130    ros_gamma(1) =  4.35866521508459d-01;
2131    ros_gamma(2) = - 1.78292094614483d+00;
2132    ros_gamma(3) = - 2.46541900496934d+00;
2133    ros_gamma(4) = - 8.05529997906370d-01;
2134
2135
2136!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
2137!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
2138    ros_newf(1) = .TRUE.
2139    ros_newf(2) = .TRUE.
2140    ros_newf(3) = .TRUE.
2141    ros_newf(4) = .TRUE.
2142
2143!~~~> ros_elo  = estimator of local order - the minimum between the
2144!        main and the embedded scheme orders plus 1
2145    ros_elo = 3.0_dp
2146
2147  END SUBROUTINE rang3
2148!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2149
2150!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2151!   End of the set of internal Rosenbrock subroutines
2152!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2153END SUBROUTINE rosenbrock
2154 
2155SUBROUTINE funtemplate( t, y, ydot)
2156!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2157!  Template for the ODE function call.
2158!  Updates the rate coefficients (and possibly the fixed species) at each call
2159!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2160!~~~> input variables
2161   REAL(kind=dp):: t, y(nvar)
2162!~~~> output variables
2163   REAL(kind=dp):: ydot(nvar)
2164!~~~> local variables
2165   REAL(kind=dp):: told
2166
2167   told = time
2168   time = t
2169   CALL fun( y, fix, rconst, ydot)
2170   time = told
2171
2172END SUBROUTINE funtemplate
2173 
2174SUBROUTINE jactemplate( t, y, jcb)
2175!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2176!  Template for the ODE Jacobian call.
2177!  Updates the rate coefficients (and possibly the fixed species) at each call
2178!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2179!~~~> input variables
2180    REAL(kind=dp):: t, y(nvar)
2181!~~~> output variables
2182#ifdef full_algebra   
2183    REAL(kind=dp):: jv(lu_nonzero), jcb(nvar, nvar)
2184#else
2185    REAL(kind=dp):: jcb(lu_nonzero)
2186#endif   
2187!~~~> local variables
2188    REAL(kind=dp):: told
2189#ifdef full_algebra   
2190    INTEGER :: i, j
2191#endif   
2192
2193    told = time
2194    time = t
2195#ifdef full_algebra   
2196    CALL jac_sp(y, fix, rconst, jv)
2197    DO j=1, nvar
2198      DO i=1, nvar
2199         jcb(i, j) = 0.0_dp
2200      ENDDO
2201    ENDDO
2202    DO i=1, lu_nonzero
2203       jcb(lu_irow(i), lu_icol(i)) = jv(i)
2204    ENDDO
2205#else
2206    CALL jac_sp( y, fix, rconst, jcb)
2207#endif   
2208    time = told
2209
2210END SUBROUTINE jactemplate
2211 
2212  SUBROUTINE kppdecomp( jvs, ier)                                   
2213  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2214  !        sparse lu factorization                                   
2215  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2216  ! loop expansion generated by kp4                                   
2217                                                                     
2218    INTEGER  :: ier                                                   
2219    REAL(kind=dp):: jvs(lu_nonzero), w(nvar), a             
2220    INTEGER  :: k, kk, j, jj                                         
2221                                                                     
2222    a = 0.                                                           
2223    ier = 0                                                           
2224                                                                     
2225!   i = 1
2226    RETURN                                                           
2227                                                                     
2228  END SUBROUTINE kppdecomp                                           
2229 
2230SUBROUTINE chem_gasphase_integrate (time_step_len, conc, tempi, qvapi, fakti, photo, ierrf, xnacc, xnrej, istatus, l_debug, pe, &
2231                     icntrl_i, rcntrl_i) 
2232                                                                   
2233  IMPLICIT NONE                                                     
2234                                                                   
2235  REAL(dp), INTENT(IN)                  :: time_step_len           
2236  REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: conc                   
2237  REAL(dp), DIMENSION(:, :), INTENT(IN)   :: photo                   
2238  REAL(dp), DIMENSION(:), INTENT(IN)     :: tempi                   
2239  REAL(dp), DIMENSION(:), INTENT(IN)     :: qvapi                   
2240  REAL(dp), DIMENSION(:), INTENT(IN)     :: fakti                   
2241  INTEGER, INTENT(OUT), OPTIONAL        :: ierrf(:)               
2242  INTEGER, INTENT(OUT), OPTIONAL        :: xnacc(:)               
2243  INTEGER, INTENT(OUT), OPTIONAL        :: xnrej(:)               
2244  INTEGER, INTENT(INOUT), OPTIONAL      :: istatus(:)             
2245  INTEGER, INTENT(IN), OPTIONAL         :: pe                     
2246  LOGICAL, INTENT(IN), OPTIONAL         :: l_debug                 
2247  INTEGER, DIMENSION(nkppctrl), INTENT(IN), OPTIONAL  :: icntrl_i         
2248  REAL(dp), DIMENSION(nkppctrl), INTENT(IN), OPTIONAL  :: rcntrl_i         
2249                                                                   
2250  INTEGER                                 :: k   ! loop variable     
2251  REAL(dp)                               :: dt                     
2252  INTEGER, DIMENSION(20)                :: istatus_u               
2253  INTEGER                                :: ierr_u                 
2254  INTEGER                                :: istatf                 
2255  INTEGER                                :: vl_dim_lo               
2256                                                                   
2257                                                                   
2258  IF (PRESENT (istatus)) istatus = 0                             
2259  IF (PRESENT (icntrl_i)) icntrl  = icntrl_i                     
2260  IF (PRESENT (rcntrl_i)) rcntrl  = rcntrl_i                     
2261                                                                   
2262  vl_glo = size(tempi, 1)                                           
2263                                                                   
2264  vl_dim_lo = vl_dim                                               
2265  DO k=1, vl_glo, vl_dim_lo                                           
2266    is = k                                                         
2267    ie = min(k+ vl_dim_lo-1, vl_glo)                                 
2268    vl = ie-is+ 1                                                   
2269                                                                   
2270    c(:) = conc(is, :)                                           
2271                                                                   
2272    temp = tempi(is)                                             
2273                                                                   
2274    qvap = qvapi(is)                                             
2275                                                                   
2276    fakt = fakti(is)                                             
2277                                                                   
2278    CALL initialize                                                 
2279                                                                   
2280    phot(:) = photo(is, :)                                           
2281                                                                   
2282    CALL update_rconst                                             
2283                                                                   
2284    dt = time_step_len                                             
2285                                                                   
2286    ! integrate from t=0 to t=dt                                   
2287    CALL integrate(0._dp, dt, icntrl, rcntrl, istatus_u = istatus_u, ierr_u=ierr_u)
2288                                                                   
2289                                                                   
2290   IF (PRESENT(l_debug) .AND. PRESENT(pe)) THEN                       
2291      IF (l_debug) CALL error_output(conc(is, :), ierr_u, pe)         
2292   ENDIF                                                             
2293                                                                     
2294    conc(is, :) = c(:)                                               
2295                                                                   
2296    ! RETURN diagnostic information                                 
2297                                                                   
2298    IF (PRESENT(ierrf))   ierrf(is) = ierr_u                     
2299    IF (PRESENT(xnacc))   xnacc(is) = istatus_u(4)               
2300    IF (PRESENT(xnrej))   xnrej(is) = istatus_u(5)               
2301                                                                   
2302    IF (PRESENT (istatus)) THEN                                   
2303      istatus(1:8) = istatus(1:8) + istatus_u(1:8)               
2304    ENDIF                                                         
2305                                                                   
2306  END DO                                                           
2307 
2308                                                                   
2309! Deallocate input arrays                                           
2310                                                                   
2311                                                                   
2312  data_loaded = .FALSE.                                             
2313                                                                   
2314  RETURN                                                           
2315END SUBROUTINE chem_gasphase_integrate                                       
2316
2317END MODULE chem_gasphase_mod
2318 
Note: See TracBrowser for help on using the repository browser.