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: 81.2 KB

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

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

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_phstat/chem_gasphase_mod.f90 @ 4901

Download in other formats: