Home

Context Navigation

chem_gasphase_mod.f90 @ 4901

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

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

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_smog/chem_gasphase_mod.f90 @ 4901

Download in other formats: