Home

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_salsa+simple/chem_gasphase_mod.f90 @ 3606

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |