Home

Context Navigation

chem_gasphase_mod.f90 @ 4901

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

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

Context Navigation

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

Download in other formats: