source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_passive/chem_gasphase_mod.f90 @ 3789

Last change on this file since 3789 was 3789, checked in by forkel, 3 years ago

Removed unused variables from chem_gasphase_mod.f90

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