Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |