Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |