Home

Context Navigation

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

Last change on this file since 3789 was 3789, checked in by forkel, 5 years ago
Removed unused variables from chem_gasphase_mod.f90
File size: 85.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |