Home

Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

chem_gasphase_mod.f90

Last change on this file was 4841, checked in by forkel, 4 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: 81.5 KB

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

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

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_phstatp/chem_gasphase_mod.f90

Download in other formats: