Home

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_salsa+simple/chem_gasphase_mod.f90 @ 3949

Last change on this file since 3949 was 3949, checked in by forkel, 5 years ago
fix in def_salsa+simple/chem_gasphase_mod.kpp and finaling def_salsa+phstat
File size: 87.8 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |