Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |