Home

Context Navigation

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

Last change on this file since 3582 was 3566, checked in by monakurppa, 5 years ago
Branch salsa: included chemical mechanism including salsa gases
File size: 80.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |