Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |