Home

Context Navigation

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

Last change on this file since 3780 was 3780, checked in by forkel, 5 years ago
removed read from unit 10 in chemistry_model_mod.f90, added get_mechanismname
File size: 85.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |