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

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

Removed unused variables from chem_gasphase_mod.f90

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