source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_passive/chem_gasphase_mod.f90 @ 4841

Last change on this file since 4841 was 4841, checked in by forkel, 4 years ago

updated copyright statements for chem_gasphase_mod.f90. This must be done in UTIL/chemistry/gasphase_preproc/kpp4palm/templates/module_header and NOT in SOURCE/chem_gasphase_mod.f90.

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