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

Last change on this file since 3681 was 3681, checked in by hellstea, 3 years ago

Major update of pmc_interface_mod

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