Home

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_simplep/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.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |