Home

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_smog/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: 91.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |