Home

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_simplep/chem_gasphase_mod.f90 @ 3585

Last change on this file since 3585 was 3585, checked in by forkel, 5 years ago
Comments in all chem_gasphase_mod.kpp, add def_salsagas/*.f90, removed executables
File size: 85.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |