Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |