Home

Context Navigation

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

Last change on this file since 3698 was 3698, checked in by suehring, 5 years ago
Fix bad commit in gasphase_preproc
File size: 85.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |