Home

Context Navigation

source: palm/trunk/SOURCE/chem_gasphase_mod.f90 @ 4788

Last change on this file since 4788 was 4360, checked in by suehring, 5 years ago
Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update
File size: 81.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |