Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |