source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_phstat/chem_gasphase_mod.f90 @ 3789

Last change on this file since 3789 was 3789, checked in by forkel, 6 years ago

Removed unused variables from chem_gasphase_mod.f90

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