Home

Context Navigation

chem_gasphase_mod.f90 @ 4901

Last change on this file since 4901 was 4841, checked in by forkel, 3 years ago
updated copyright statements for chem_gasphase_mod.f90. This must be done in UTIL/chemistry/gasphase_preproc/kpp4palm/templates/module_header and NOT in SOURCE/chem_gasphase_mod.f90.
File size: 82.8 KB

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

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

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_salsa+phstat/chem_gasphase_mod.f90 @ 4901

Download in other formats: