Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |