source: palm/trunk/SOURCE/chem_gasphase_mod.f90 @ 3703

Last change on this file since 3703 was 3655, checked in by knoop, 6 years ago

Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update

File size: 79.7 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 :: 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:19 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 = 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 Nov 30 13:52:19 2018
193! Working directory    : /home/forkel-r/palmstuff/work/trunk20181130/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! TSTART - Integration start time
221  REAL(kind=dp):: tstart
222! ATOL - Absolute tolerance
223  REAL(kind=dp):: atol(nvar)
224! RTOL - Relative tolerance
225  REAL(kind=dp):: rtol(nvar)
226! STEPMIN - Lower bound for integration step
227  REAL(kind=dp):: stepmin
228! CFACTOR - Conversion factor for concentration units
229  REAL(kind=dp):: cfactor
230
231! INLINED global variable declarations
232
233! QVAP - Water vapor
234  REAL(kind=dp):: qvap
235! FAKT - Conversion factor
236  REAL(kind=dp):: fakt
237
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 Nov 30 13:52:19 2018
259! Working directory    : /home/forkel-r/palmstuff/work/trunk20181130/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 Nov 30 13:52:19 2018
303! Working directory    : /home/forkel-r/palmstuff/work/trunk20181130/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 Nov 30 13:52:19 2018
365! Working directory    : /home/forkel-r/palmstuff/work/trunk20181130/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 Nov 30 13:52:19 2018
391! Working directory    : /home/forkel-r/palmstuff/work/trunk20181130/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 Nov 30 13:52:19 2018
449! Working directory    : /home/forkel-r/palmstuff/work/trunk20181130/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 Nov 30 13:52:19 2018
476! Working directory    : /home/forkel-r/palmstuff/work/trunk20181130/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 Nov 30 13:52:19 2018
503! Working directory    : /home/forkel-r/palmstuff/work/trunk20181130/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 Nov 30 13:52:19 2018
532! Working directory    : /home/forkel-r/palmstuff/work/trunk20181130/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 Nov 30 13:52:19 2018
558! Working directory    : /home/forkel-r/palmstuff/work/trunk20181130/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            chem_gasphase_integrate
665    MODULE PROCEDURE   chem_gasphase_integrate
666  END INTERFACE        chem_gasphase_integrate
667 
668 
669 CONTAINS
670 
671SUBROUTINE initialize()
672
673
674  INTEGER         :: j, k
675
676  INTEGER :: i
677  REAL(kind=dp):: x
678  k = is
679  cfactor = 1.000000e+00_dp
680
681  x = (0.) * cfactor
682  DO i = 1 , nvar
683  ENDDO
684
685  x = (0.) * cfactor
686  DO i = 1 , nfix
687    fix(i) = x
688  ENDDO
689
690! constant rate coefficients
691! END constant rate coefficients
692
693! INLINED initializations
694
695! End INLINED initializations
696
697     
698END SUBROUTINE initialize
699 
700SUBROUTINE integrate( tin, tout, &
701  icntrl_u, rcntrl_u, istatus_u, rstatus_u, ierr_u)
702
703
704   REAL(kind=dp), INTENT(IN):: tin  ! start time
705   REAL(kind=dp), INTENT(IN):: tout ! END time
706   ! OPTIONAL input PARAMETERs and statistics
707   INTEGER,      INTENT(IN), OPTIONAL :: icntrl_u(20)
708   REAL(kind=dp), INTENT(IN), OPTIONAL :: rcntrl_u(20)
709   INTEGER,      INTENT(OUT), OPTIONAL :: istatus_u(20)
710   REAL(kind=dp), INTENT(OUT), OPTIONAL :: rstatus_u(20)
711   INTEGER,      INTENT(OUT), OPTIONAL :: ierr_u
712
713   REAL(kind=dp):: rcntrl(20), rstatus(20)
714   INTEGER       :: icntrl(20), istatus(20), ierr
715
716   INTEGER, SAVE :: ntotal = 0
717
718   icntrl(:) = 0
719   rcntrl(:) = 0.0_dp
720   istatus(:) = 0
721   rstatus(:) = 0.0_dp
722
723    !~~~> fine-tune the integrator:
724   icntrl(1) = 0      ! 0 - non- autonomous, 1 - autonomous
725   icntrl(2) = 0      ! 0 - vector tolerances, 1 - scalars
726
727   ! IF OPTIONAL PARAMETERs are given, and IF they are >0,
728   ! THEN they overwrite default settings.
729   IF (PRESENT(icntrl_u))THEN
730     WHERE(icntrl_u(:)> 0)icntrl(:) = icntrl_u(:)
731   ENDIF
732   IF (PRESENT(rcntrl_u))THEN
733     WHERE(rcntrl_u(:)> 0)rcntrl(:) = rcntrl_u(:)
734   ENDIF
735
736
737   CALL rosenbrock(nvar, var, tin, tout,  &
738         atol, rtol,               &
739         rcntrl, icntrl, rstatus, istatus, ierr)
740
741   !~~~> debug option: show no of steps
742   ! ntotal = ntotal + istatus(nstp)
743   ! PRINT*,'NSTEPS=',ISTATUS(Nstp),' (',Ntotal,')','  O3=',VAR(ind_O3)
744
745   stepmin = rstatus(nhexit)
746   ! IF OPTIONAL PARAMETERs are given for output they
747   ! are updated with the RETURN information
748   IF (PRESENT(istatus_u))istatus_u(:) = istatus(:)
749   IF (PRESENT(rstatus_u))rstatus_u(:) = rstatus(:)
750   IF (PRESENT(ierr_u))  ierr_u       = ierr
751
752END SUBROUTINE integrate
753 
754SUBROUTINE fun(v, f, rct, vdot)
755
756! V - Concentrations of variable species (local)
757  REAL(kind=dp):: v(nvar)
758! F - Concentrations of fixed species (local)
759  REAL(kind=dp):: f(nfix)
760! RCT - Rate constants (local)
761  REAL(kind=dp):: rct(nreact)
762! Vdot - Time derivative of variable species concentrations
763  REAL(kind=dp):: vdot(nvar)
764
765
766! Computation of equation rates
767  a(1) = rct(1) * v(3)
768  a(2) = rct(2) * v(2) * v(4)
769
770! Aggregate function
771  vdot(1) = 0
772  vdot(2) = a(1) - a(2)
773  vdot(3) = - a(1) + a(2)
774  vdot(4) = a(1) - a(2)
775     
776END SUBROUTINE fun
777 
778SUBROUTINE kppsolve(jvs, x)
779
780! JVS - sparse Jacobian of variables
781  REAL(kind=dp):: jvs(lu_nonzero)
782! X - Vector for variables
783  REAL(kind=dp):: x(nvar)
784
785  x(3) = x(3) - jvs(5) * x(2)
786  x(4) = x(4) - jvs(8) * x(2) - jvs(9) * x(3)
787  x(4) = x(4) / jvs(10)
788  x(3) = (x(3) - jvs(7) * x(4)) /(jvs(6))
789  x(2) = (x(2) - jvs(3) * x(3) - jvs(4) * x(4)) /(jvs(2))
790  x(1) = x(1) / jvs(1)
791     
792END SUBROUTINE kppsolve
793 
794SUBROUTINE jac_sp(v, f, rct, jvs)
795
796! V - Concentrations of variable species (local)
797  REAL(kind=dp):: v(nvar)
798! F - Concentrations of fixed species (local)
799  REAL(kind=dp):: f(nfix)
800! RCT - Rate constants (local)
801  REAL(kind=dp):: rct(nreact)
802! JVS - sparse Jacobian of variables
803  REAL(kind=dp):: jvs(lu_nonzero)
804
805
806! Local variables
807! B - Temporary array
808  REAL(kind=dp):: b(4)
809
810! B(1) = dA(1)/dV(3)
811  b(1) = rct(1)
812! B(2) = dA(2)/dV(2)
813  b(2) = rct(2) * v(4)
814! B(3) = dA(2)/dV(4)
815  b(3) = rct(2) * v(2)
816! B(4) = dA(3)/dV(1)
817  b(4) = rct(3)
818
819! Construct the Jacobian terms from B's
820! JVS(1) = Jac_FULL(1,1)
821  jvs(1) = 0
822! JVS(2) = Jac_FULL(2,2)
823  jvs(2) = - b(2)
824! JVS(3) = Jac_FULL(2,3)
825  jvs(3) = b(1)
826! JVS(4) = Jac_FULL(2,4)
827  jvs(4) = - b(3)
828! JVS(5) = Jac_FULL(3,2)
829  jvs(5) = b(2)
830! JVS(6) = Jac_FULL(3,3)
831  jvs(6) = - b(1)
832! JVS(7) = Jac_FULL(3,4)
833  jvs(7) = b(3)
834! JVS(8) = Jac_FULL(4,2)
835  jvs(8) = - b(2)
836! JVS(9) = Jac_FULL(4,3)
837  jvs(9) = b(1)
838! JVS(10) = Jac_FULL(4,4)
839  jvs(10) = - b(3)
840     
841END SUBROUTINE jac_sp
842 
843  elemental REAL(kind=dp)FUNCTION k_arr (k_298, tdep, temp)
844    ! arrhenius FUNCTION
845
846    REAL,    INTENT(IN):: k_298 ! k at t = 298.15k
847    REAL,    INTENT(IN):: tdep  ! temperature dependence
848    REAL(kind=dp), INTENT(IN):: temp  ! temperature
849
850    intrinsic exp
851
852    k_arr = k_298 * exp(tdep* (1._dp/temp- 3.3540e-3_dp))! 1/298.15=3.3540e-3
853
854  END FUNCTION k_arr
855 
856SUBROUTINE update_rconst()
857 INTEGER         :: k
858
859  k = is
860
861! Begin INLINED RCONST
862
863
864! End INLINED RCONST
865
866  rconst(1) = (phot(j_no2))
867  rconst(2) = (arr2(1.8e-12_dp , 1370.0_dp , temp))
868  rconst(3) = (1.0_dp)
869     
870END SUBROUTINE update_rconst
871 
872!  END FUNCTION ARR2
873REAL(kind=dp)FUNCTION arr2( a0, b0, temp)
874    REAL(kind=dp):: temp
875    REAL(kind=dp):: a0, b0
876    arr2 = a0 * exp( - b0 / temp)
877END FUNCTION arr2
878 
879SUBROUTINE initialize_kpp_ctrl(status)
880
881
882  ! i/o
883  INTEGER,         INTENT(OUT):: status
884
885  ! local
886  REAL(dp):: tsum
887  INTEGER  :: i
888
889  ! check fixed time steps
890  tsum = 0.0_dp
891  DO i=1, nmaxfixsteps
892     IF (t_steps(i)< tiny(0.0_dp))exit
893     tsum = tsum + t_steps(i)
894  ENDDO
895
896  nfsteps = i- 1
897
898  l_fixed_step = (nfsteps > 0).and.((tsum - 1.0)< tiny(0.0_dp))
899
900  IF (l_vector)THEN
901     WRITE(*,*) ' MODE             : VECTOR (LENGTH=',VL_DIM,')'
902  ELSE
903     WRITE(*,*) ' MODE             : SCALAR'
904  ENDIF
905  !
906  WRITE(*,*) ' DE-INDEXING MODE :',I_LU_DI
907  !
908  WRITE(*,*) ' ICNTRL           : ',icntrl
909  WRITE(*,*) ' RCNTRL           : ',rcntrl
910  !
911  ! note: this is ONLY meaningful for vectorized (kp4)rosenbrock- methods
912  IF (l_vector)THEN
913     IF (l_fixed_step)THEN
914        WRITE(*,*) ' TIME STEPS       : FIXED (',t_steps(1:nfsteps),')'
915     ELSE
916        WRITE(*,*) ' TIME STEPS       : AUTOMATIC'
917     ENDIF
918  ELSE
919     WRITE(*,*) ' TIME STEPS       : AUTOMATIC '//&
920          &'(t_steps (CTRL_KPP) ignored in SCALAR MODE)'
921  ENDIF
922  ! mz_pj_20070531-
923
924  status = 0
925
926
927END SUBROUTINE initialize_kpp_ctrl
928 
929SUBROUTINE error_output(c, ierr, pe)
930
931
932  INTEGER, INTENT(IN):: ierr
933  INTEGER, INTENT(IN):: pe
934  REAL(dp), DIMENSION(:), INTENT(IN):: c
935
936  write(6,*) 'ERROR in chem_gasphase_mod ',ierr,C(1)
937
938
939END SUBROUTINE error_output
940 
941      SUBROUTINE wscal(n, alpha, x, incx)
942!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
943!     constant times a vector: x(1:N) <- Alpha*x(1:N)
944!     only for incX=incY=1
945!     after BLAS
946!     replace this by the function from the optimized BLAS implementation:
947!         CALL SSCAL(N,Alpha,X,1) or  CALL DSCAL(N,Alpha,X,1)
948!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
949
950      INTEGER  :: i, incx, m, mp1, n
951      REAL(kind=dp) :: x(n), alpha
952      REAL(kind=dp), PARAMETER  :: zero=0.0_dp, one=1.0_dp
953
954      IF (alpha .eq. one)RETURN
955      IF (n .le. 0)RETURN
956
957      m = mod(n, 5)
958      IF ( m .ne. 0)THEN
959        IF (alpha .eq. (- one))THEN
960          DO i = 1, m
961            x(i) = - x(i)
962          ENDDO
963        ELSEIF (alpha .eq. zero)THEN
964          DO i = 1, m
965            x(i) = zero
966          ENDDO
967        ELSE
968          DO i = 1, m
969            x(i) = alpha* x(i)
970          ENDDO
971        ENDIF
972        IF ( n .lt. 5)RETURN
973      ENDIF
974      mp1 = m + 1
975      IF (alpha .eq. (- one))THEN
976        DO i = mp1, n, 5
977          x(i)   = - x(i)
978          x(i + 1) = - x(i + 1)
979          x(i + 2) = - x(i + 2)
980          x(i + 3) = - x(i + 3)
981          x(i + 4) = - x(i + 4)
982        ENDDO
983      ELSEIF (alpha .eq. zero)THEN
984        DO i = mp1, n, 5
985          x(i)   = zero
986          x(i + 1) = zero
987          x(i + 2) = zero
988          x(i + 3) = zero
989          x(i + 4) = zero
990        ENDDO
991      ELSE
992        DO i = mp1, n, 5
993          x(i)   = alpha* x(i)
994          x(i + 1) = alpha* x(i + 1)
995          x(i + 2) = alpha* x(i + 2)
996          x(i + 3) = alpha* x(i + 3)
997          x(i + 4) = alpha* x(i + 4)
998        ENDDO
999      ENDIF
1000
1001      END SUBROUTINE wscal
1002 
1003      SUBROUTINE waxpy(n, alpha, x, incx, y, incy)
1004!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1005!     constant times a vector plus a vector: y <- y + Alpha*x
1006!     only for incX=incY=1
1007!     after BLAS
1008!     replace this by the function from the optimized BLAS implementation:
1009!         CALL SAXPY(N,Alpha,X,1,Y,1) or  CALL DAXPY(N,Alpha,X,1,Y,1)
1010!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1011
1012      INTEGER  :: i, incx, incy, m, mp1, n
1013      REAL(kind=dp):: x(n), y(n), alpha
1014      REAL(kind=dp), PARAMETER :: zero = 0.0_dp
1015
1016      IF (alpha .eq. zero)RETURN
1017      IF (n .le. 0)RETURN
1018
1019      m = mod(n, 4)
1020      IF ( m .ne. 0)THEN
1021        DO i = 1, m
1022          y(i) = y(i) + alpha* x(i)
1023        ENDDO
1024        IF ( n .lt. 4)RETURN
1025      ENDIF
1026      mp1 = m + 1
1027      DO i = mp1, n, 4
1028        y(i) = y(i) + alpha* x(i)
1029        y(i + 1) = y(i + 1) + alpha* x(i + 1)
1030        y(i + 2) = y(i + 2) + alpha* x(i + 2)
1031        y(i + 3) = y(i + 3) + alpha* x(i + 3)
1032      ENDDO
1033     
1034      END SUBROUTINE waxpy
1035 
1036SUBROUTINE rosenbrock(n, y, tstart, tend, &
1037           abstol, reltol,             &
1038           rcntrl, icntrl, rstatus, istatus, ierr)
1039!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1040!
1041!    Solves the system y'=F(t,y) using a Rosenbrock method defined by:
1042!
1043!     G = 1/(H*gamma(1)) - Jac(t0,Y0)
1044!     T_i = t0 + Alpha(i)*H
1045!     Y_i = Y0 + \sum_{j=1}^{i-1} A(i,j)*K_j
1046!     G *K_i = Fun( T_i,Y_i)+ \sum_{j=1}^S C(i,j)/H *K_j +
1047!         gamma(i)*dF/dT(t0,Y0)
1048!     Y1 = Y0 + \sum_{j=1}^S M(j)*K_j
1049!
1050!    For details on Rosenbrock methods and their implementation consult:
1051!      E. Hairer and G. Wanner
1052!      "Solving ODEs II. Stiff and differential-algebraic problems".
1053!      Springer series in computational mathematics,Springer-Verlag,1996.
1054!    The codes contained in the book inspired this implementation.
1055!
1056!    (C)  Adrian Sandu,August 2004
1057!    Virginia Polytechnic Institute and State University
1058!    Contact: sandu@cs.vt.edu
1059!    Revised by Philipp Miehe and Adrian Sandu,May 2006                 
1060!    This implementation is part of KPP - the Kinetic PreProcessor
1061!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1062!
1063!~~~>   input arguments:
1064!
1065!-     y(n)  = vector of initial conditions (at t=tstart)
1066!-    [tstart, tend]  = time range of integration
1067!     (if Tstart>Tend the integration is performed backwards in time)
1068!-    reltol, abstol = user precribed accuracy
1069!- SUBROUTINE  fun( t, y, ydot) = ode FUNCTION,
1070!                       returns Ydot = Y' = F(T,Y)
1071!- SUBROUTINE  jac( t, y, jcb) = jacobian of the ode FUNCTION,
1072!                       returns Jcb = dFun/dY
1073!-    icntrl(1:20)  = INTEGER inputs PARAMETERs
1074!-    rcntrl(1:20)  = REAL inputs PARAMETERs
1075!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1076!
1077!~~~>     output arguments:
1078!
1079!-    y(n)  - > vector of final states (at t- >tend)
1080!-    istatus(1:20) - > INTEGER output PARAMETERs
1081!-    rstatus(1:20) - > REAL output PARAMETERs
1082!-    ierr            - > job status upon RETURN
1083!                        success (positive value) or
1084!                        failure (negative value)
1085!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1086!
1087!~~~>     input PARAMETERs:
1088!
1089!    Note: For input parameters equal to zero the default values of the
1090!       corresponding variables are used.
1091!
1092!    ICNTRL(1) = 1: F = F(y)   Independent of T (AUTONOMOUS)
1093!              = 0: F = F(t,y) Depends on T (NON-AUTONOMOUS)
1094!
1095!    ICNTRL(2) = 0: AbsTol,RelTol are N-dimensional vectors
1096!              = 1: AbsTol,RelTol are scalars
1097!
1098!    ICNTRL(3)  -> selection of a particular Rosenbrock method
1099!        = 0 :    Rodas3 (default)
1100!        = 1 :    Ros2
1101!        = 2 :    Ros3
1102!        = 3 :    Ros4
1103!        = 4 :    Rodas3
1104!        = 5 :    Rodas4
1105!
1106!    ICNTRL(4)  -> maximum number of integration steps
1107!        For ICNTRL(4) =0) the default value of 100000 is used
1108!
1109!    RCNTRL(1)  -> Hmin,lower bound for the integration step size
1110!          It is strongly recommended to keep Hmin = ZERO
1111!    RCNTRL(2)  -> Hmax,upper bound for the integration step size
1112!    RCNTRL(3)  -> Hstart,starting value for the integration step size
1113!
1114!    RCNTRL(4)  -> FacMin,lower bound on step decrease factor (default=0.2)
1115!    RCNTRL(5)  -> FacMax,upper bound on step increase factor (default=6)
1116!    RCNTRL(6)  -> FacRej,step decrease factor after multiple rejections
1117!                          (default=0.1)
1118!    RCNTRL(7)  -> FacSafe,by which the new step is slightly smaller
1119!         than the predicted value  (default=0.9)
1120!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1121!
1122!
1123!    OUTPUT ARGUMENTS:
1124!    -----------------
1125!
1126!    T           -> T value for which the solution has been computed
1127!                     (after successful return T=Tend).
1128!
1129!    Y(N)        -> Numerical solution at T
1130!
1131!    IDID        -> Reports on successfulness upon return:
1132!                    = 1 for success
1133!                    < 0 for error (value equals error code)
1134!
1135!    ISTATUS(1)  -> No. of function calls
1136!    ISTATUS(2)  -> No. of jacobian calls
1137!    ISTATUS(3)  -> No. of steps
1138!    ISTATUS(4)  -> No. of accepted steps
1139!    ISTATUS(5)  -> No. of rejected steps (except at very beginning)
1140!    ISTATUS(6)  -> No. of LU decompositions
1141!    ISTATUS(7)  -> No. of forward/backward substitutions
1142!    ISTATUS(8)  -> No. of singular matrix decompositions
1143!
1144!    RSTATUS(1)  -> Texit,the time corresponding to the
1145!                     computed Y upon return
1146!    RSTATUS(2)  -> Hexit,last accepted step before exit
1147!    RSTATUS(3)  -> Hnew,last predicted step (not yet taken)
1148!                   For multiple restarts,use Hnew as Hstart
1149!                     in the subsequent run
1150!
1151!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1152
1153
1154!~~~>  arguments
1155   INTEGER,      INTENT(IN)  :: n
1156   REAL(kind=dp), INTENT(INOUT):: y(n)
1157   REAL(kind=dp), INTENT(IN)  :: tstart, tend
1158   REAL(kind=dp), INTENT(IN)  :: abstol(n), reltol(n)
1159   INTEGER,      INTENT(IN)  :: icntrl(20)
1160   REAL(kind=dp), INTENT(IN)  :: rcntrl(20)
1161   INTEGER,      INTENT(INOUT):: istatus(20)
1162   REAL(kind=dp), INTENT(INOUT):: rstatus(20)
1163   INTEGER, INTENT(OUT) :: ierr
1164!~~~>  PARAMETERs of the rosenbrock method, up to 6 stages
1165   INTEGER ::  ros_s, rosmethod
1166   INTEGER, PARAMETER :: rs2=1, rs3=2, rs4=3, rd3=4, rd4=5, rg3=6
1167   REAL(kind=dp):: ros_a(15), ros_c(15), ros_m(6), ros_e(6), &
1168                    ros_alpha(6), ros_gamma(6), ros_elo
1169   LOGICAL :: ros_newf(6)
1170   CHARACTER(len=12):: ros_name
1171!~~~>  local variables
1172   REAL(kind=dp):: roundoff, facmin, facmax, facrej, facsafe
1173   REAL(kind=dp):: hmin, hmax, hstart
1174   REAL(kind=dp):: texit
1175   INTEGER       :: i, uplimtol, max_no_steps
1176   LOGICAL       :: autonomous, vectortol
1177!~~~>   PARAMETERs
1178   REAL(kind=dp), PARAMETER :: zero = 0.0_dp, one  = 1.0_dp
1179   REAL(kind=dp), PARAMETER :: deltamin = 1.0e-5_dp
1180
1181!~~~>  initialize statistics
1182   istatus(1:8) = 0
1183   rstatus(1:3) = zero
1184
1185!~~~>  autonomous or time dependent ode. default is time dependent.
1186   autonomous = .not.(icntrl(1) == 0)
1187
1188!~~~>  for scalar tolerances (icntrl(2).ne.0) the code uses abstol(1)and reltol(1)
1189!   For Vector tolerances (ICNTRL(2) == 0) the code uses AbsTol(1:N) and RelTol(1:N)
1190   IF (icntrl(2) == 0)THEN
1191      vectortol = .TRUE.
1192      uplimtol  = n
1193   ELSE
1194      vectortol = .FALSE.
1195      uplimtol  = 1
1196   ENDIF
1197
1198!~~~>   initialize the particular rosenbrock method selected
1199   select CASE (icntrl(3))
1200     CASE (1)
1201       CALL ros2
1202     CASE (2)
1203       CALL ros3
1204     CASE (3)
1205       CALL ros4
1206     CASE (0, 4)
1207       CALL rodas3
1208     CASE (5)
1209       CALL rodas4
1210     CASE (6)
1211       CALL rang3
1212     CASE default
1213       PRINT *,'Unknown Rosenbrock method: ICNTRL(3) =',ICNTRL(3) 
1214       CALL ros_errormsg(- 2, tstart, zero, ierr)
1215       RETURN
1216   END select
1217
1218!~~~>   the maximum number of steps admitted
1219   IF (icntrl(4) == 0)THEN
1220      max_no_steps = 200000
1221   ELSEIF (icntrl(4)> 0)THEN
1222      max_no_steps=icntrl(4)
1223   ELSE
1224      PRINT *,'User-selected max no. of steps: ICNTRL(4) =',ICNTRL(4)
1225      CALL ros_errormsg(- 1, tstart, zero, ierr)
1226      RETURN
1227   ENDIF
1228
1229!~~~>  unit roundoff (1+ roundoff>1)
1230   roundoff = epsilon(one)
1231
1232!~~~>  lower bound on the step size: (positive value)
1233   IF (rcntrl(1) == zero)THEN
1234      hmin = zero
1235   ELSEIF (rcntrl(1)> zero)THEN
1236      hmin = rcntrl(1)
1237   ELSE
1238      PRINT *,'User-selected Hmin: RCNTRL(1) =',RCNTRL(1)
1239      CALL ros_errormsg(- 3, tstart, zero, ierr)
1240      RETURN
1241   ENDIF
1242!~~~>  upper bound on the step size: (positive value)
1243   IF (rcntrl(2) == zero)THEN
1244      hmax = abs(tend-tstart)
1245   ELSEIF (rcntrl(2)> zero)THEN
1246      hmax = min(abs(rcntrl(2)), abs(tend-tstart))
1247   ELSE
1248      PRINT *,'User-selected Hmax: RCNTRL(2) =',RCNTRL(2)
1249      CALL ros_errormsg(- 3, tstart, zero, ierr)
1250      RETURN
1251   ENDIF
1252!~~~>  starting step size: (positive value)
1253   IF (rcntrl(3) == zero)THEN
1254      hstart = max(hmin, deltamin)
1255   ELSEIF (rcntrl(3)> zero)THEN
1256      hstart = min(abs(rcntrl(3)), abs(tend-tstart))
1257   ELSE
1258      PRINT *,'User-selected Hstart: RCNTRL(3) =',RCNTRL(3)
1259      CALL ros_errormsg(- 3, tstart, zero, ierr)
1260      RETURN
1261   ENDIF
1262!~~~>  step size can be changed s.t.  facmin < hnew/hold < facmax
1263   IF (rcntrl(4) == zero)THEN
1264      facmin = 0.2_dp
1265   ELSEIF (rcntrl(4)> zero)THEN
1266      facmin = rcntrl(4)
1267   ELSE
1268      PRINT *,'User-selected FacMin: RCNTRL(4) =',RCNTRL(4)
1269      CALL ros_errormsg(- 4, tstart, zero, ierr)
1270      RETURN
1271   ENDIF
1272   IF (rcntrl(5) == zero)THEN
1273      facmax = 6.0_dp
1274   ELSEIF (rcntrl(5)> zero)THEN
1275      facmax = rcntrl(5)
1276   ELSE
1277      PRINT *,'User-selected FacMax: RCNTRL(5) =',RCNTRL(5)
1278      CALL ros_errormsg(- 4, tstart, zero, ierr)
1279      RETURN
1280   ENDIF
1281!~~~>   facrej: factor to decrease step after 2 succesive rejections
1282   IF (rcntrl(6) == zero)THEN
1283      facrej = 0.1_dp
1284   ELSEIF (rcntrl(6)> zero)THEN
1285      facrej = rcntrl(6)
1286   ELSE
1287      PRINT *,'User-selected FacRej: RCNTRL(6) =',RCNTRL(6)
1288      CALL ros_errormsg(- 4, tstart, zero, ierr)
1289      RETURN
1290   ENDIF
1291!~~~>   facsafe: safety factor in the computation of new step size
1292   IF (rcntrl(7) == zero)THEN
1293      facsafe = 0.9_dp
1294   ELSEIF (rcntrl(7)> zero)THEN
1295      facsafe = rcntrl(7)
1296   ELSE
1297      PRINT *,'User-selected FacSafe: RCNTRL(7) =',RCNTRL(7)
1298      CALL ros_errormsg(- 4, tstart, zero, ierr)
1299      RETURN
1300   ENDIF
1301!~~~>  check IF tolerances are reasonable
1302    DO i=1, uplimtol
1303      IF ((abstol(i)<= zero).or. (reltol(i)<= 10.0_dp* roundoff)&
1304         .or. (reltol(i)>= 1.0_dp))THEN
1305        PRINT *,' AbsTol(',i,') = ',AbsTol(i)
1306        PRINT *,' RelTol(',i,') = ',RelTol(i)
1307        CALL ros_errormsg(- 5, tstart, zero, ierr)
1308        RETURN
1309      ENDIF
1310    ENDDO
1311
1312
1313!~~~>  CALL rosenbrock method
1314   CALL ros_integrator(y, tstart, tend, texit,  &
1315        abstol, reltol,                         &
1316!  Integration parameters
1317        autonomous, vectortol, max_no_steps,    &
1318        roundoff, hmin, hmax, hstart,           &
1319        facmin, facmax, facrej, facsafe,        &
1320!  Error indicator
1321        ierr)
1322
1323!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1324CONTAINS !  SUBROUTINEs internal to rosenbrock
1325!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1326   
1327!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
1328 SUBROUTINE ros_errormsg(code, t, h, ierr)
1329!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
1330!    Handles all error messages
1331!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
1332 
1333   REAL(kind=dp), INTENT(IN):: t, h
1334   INTEGER, INTENT(IN) :: code
1335   INTEGER, INTENT(OUT):: ierr
1336   
1337   ierr = code
1338   print * , &
1339     'Forced exit from Rosenbrock due to the following error:' 
1340     
1341   select CASE (code)
1342    CASE (- 1) 
1343      PRINT *,'--> Improper value for maximal no of steps'
1344    CASE (- 2) 
1345      PRINT *,'--> Selected Rosenbrock method not implemented'
1346    CASE (- 3) 
1347      PRINT *,'--> Hmin/Hmax/Hstart must be positive'
1348    CASE (- 4) 
1349      PRINT *,'--> FacMin/FacMax/FacRej must be positive'
1350    CASE (- 5)
1351      PRINT *,'--> Improper tolerance values'
1352    CASE (- 6)
1353      PRINT *,'--> No of steps exceeds maximum bound'
1354    CASE (- 7)
1355      PRINT *,'--> Step size too small: T + 10*H = T',&
1356            ' or H < Roundoff'
1357    CASE (- 8) 
1358      PRINT *,'--> Matrix is repeatedly singular'
1359    CASE default
1360      PRINT *,'Unknown Error code: ',Code
1361   END select
1362   
1363   print * , "t=", t, "and h=", h
1364     
1365 END SUBROUTINE ros_errormsg
1366
1367!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1368 SUBROUTINE ros_integrator (y, tstart, tend, t, &
1369        abstol, reltol,                         &
1370!~~~> integration PARAMETERs
1371        autonomous, vectortol, max_no_steps,    &
1372        roundoff, hmin, hmax, hstart,           &
1373        facmin, facmax, facrej, facsafe,        &
1374!~~~> error indicator
1375        ierr)
1376!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1377!   Template for the implementation of a generic Rosenbrock method
1378!      defined by ros_S (no of stages)
1379!      and its coefficients ros_{A,C,M,E,Alpha,Gamma}
1380!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1381
1382
1383!~~~> input: the initial condition at tstart; output: the solution at t
1384   REAL(kind=dp), INTENT(INOUT):: y(n)
1385!~~~> input: integration interval
1386   REAL(kind=dp), INTENT(IN):: tstart, tend
1387!~~~> output: time at which the solution is RETURNed (t=tendIF success)
1388   REAL(kind=dp), INTENT(OUT)::  t
1389!~~~> input: tolerances
1390   REAL(kind=dp), INTENT(IN)::  abstol(n), reltol(n)
1391!~~~> input: integration PARAMETERs
1392   LOGICAL, INTENT(IN):: autonomous, vectortol
1393   REAL(kind=dp), INTENT(IN):: hstart, hmin, hmax
1394   INTEGER, INTENT(IN):: max_no_steps
1395   REAL(kind=dp), INTENT(IN):: roundoff, facmin, facmax, facrej, facsafe
1396!~~~> output: error indicator
1397   INTEGER, INTENT(OUT):: ierr
1398! ~~~~ Local variables
1399   REAL(kind=dp):: ynew(n), fcn0(n), fcn(n)
1400   REAL(kind=dp):: k(n* ros_s), dfdt(n)
1401#ifdef full_algebra   
1402   REAL(kind=dp):: jac0(n, n), ghimj(n, n)
1403#else
1404   REAL(kind=dp):: jac0(lu_nonzero), ghimj(lu_nonzero)
1405#endif
1406   REAL(kind=dp):: h, hnew, hc, hg, fac, tau
1407   REAL(kind=dp):: err, yerr(n)
1408   INTEGER :: pivot(n), direction, ioffset, j, istage
1409   LOGICAL :: rejectlasth, rejectmoreh, singular
1410!~~~>  local PARAMETERs
1411   REAL(kind=dp), PARAMETER :: zero = 0.0_dp, one  = 1.0_dp
1412   REAL(kind=dp), PARAMETER :: deltamin = 1.0e-5_dp
1413!~~~>  locally called FUNCTIONs
1414!    REAL(kind=dp) WLAMCH
1415!    EXTERNAL WLAMCH
1416!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1417
1418
1419!~~~>  initial preparations
1420   t = tstart
1421   rstatus(nhexit) = zero
1422   h = min( max(abs(hmin), abs(hstart)), abs(hmax))
1423   IF (abs(h)<= 10.0_dp* roundoff)h = deltamin
1424
1425   IF (tend  >=  tstart)THEN
1426     direction = + 1
1427   ELSE
1428     direction = - 1
1429   ENDIF
1430   h = direction* h
1431
1432   rejectlasth=.FALSE.
1433   rejectmoreh=.FALSE.
1434
1435!~~~> time loop begins below
1436
1437timeloop: DO WHILE((direction > 0).and.((t- tend) + roundoff <= zero)&
1438       .or. (direction < 0).and.((tend-t) + roundoff <= zero))
1439
1440   IF (istatus(nstp)> max_no_steps)THEN  ! too many steps
1441      CALL ros_errormsg(- 6, t, h, ierr)
1442      RETURN
1443   ENDIF
1444   IF (((t+ 0.1_dp* h) == t).or.(h <= roundoff))THEN  ! step size too small
1445      CALL ros_errormsg(- 7, t, h, ierr)
1446      RETURN
1447   ENDIF
1448
1449!~~~>  limit h IF necessary to avoid going beyond tend
1450   h = min(h, abs(tend-t))
1451
1452!~~~>   compute the FUNCTION at current time
1453   CALL funtemplate(t, y, fcn0)
1454   istatus(nfun) = istatus(nfun) + 1
1455
1456!~~~>  compute the FUNCTION derivative with respect to t
1457   IF (.not.autonomous)THEN
1458      CALL ros_funtimederivative(t, roundoff, y, &
1459                fcn0, dfdt)
1460   ENDIF
1461
1462!~~~>   compute the jacobian at current time
1463   CALL jactemplate(t, y, jac0)
1464   istatus(njac) = istatus(njac) + 1
1465
1466!~~~>  repeat step calculation until current step accepted
1467untilaccepted: do
1468
1469   CALL ros_preparematrix(h, direction, ros_gamma(1), &
1470          jac0, ghimj, pivot, singular)
1471   IF (singular)THEN ! more than 5 consecutive failed decompositions
1472       CALL ros_errormsg(- 8, t, h, ierr)
1473       RETURN
1474   ENDIF
1475
1476!~~~>   compute the stages
1477stage: DO istage = 1, ros_s
1478
1479      ! current istage offset. current istage vector is k(ioffset+ 1:ioffset+ n)
1480       ioffset = n* (istage-1)
1481
1482      ! for the 1st istage the FUNCTION has been computed previously
1483       IF (istage == 1)THEN
1484         !slim: CALL wcopy(n, fcn0, 1, fcn, 1)
1485       fcn(1:n) = fcn0(1:n)
1486      ! istage>1 and a new FUNCTION evaluation is needed at the current istage
1487       ELSEIF(ros_newf(istage))THEN
1488         !slim: CALL wcopy(n, y, 1, ynew, 1)
1489       ynew(1:n) = y(1:n)
1490         DO j = 1, istage-1
1491           CALL waxpy(n, ros_a((istage-1) * (istage-2) /2+ j), &
1492            k(n* (j- 1) + 1), 1, ynew, 1)
1493         ENDDO
1494         tau = t + ros_alpha(istage) * direction* h
1495         CALL funtemplate(tau, ynew, fcn)
1496         istatus(nfun) = istatus(nfun) + 1
1497       ENDIF ! IF istage == 1 ELSEIF ros_newf(istage)
1498       !slim: CALL wcopy(n, fcn, 1, k(ioffset+ 1), 1)
1499       k(ioffset+ 1:ioffset+ n) = fcn(1:n)
1500       DO j = 1, istage-1
1501         hc = ros_c((istage-1) * (istage-2) /2+ j) /(direction* h)
1502         CALL waxpy(n, hc, k(n* (j- 1) + 1), 1, k(ioffset+ 1), 1)
1503       ENDDO
1504       IF ((.not. autonomous).and.(ros_gamma(istage).ne.zero))THEN
1505         hg = direction* h* ros_gamma(istage)
1506         CALL waxpy(n, hg, dfdt, 1, k(ioffset+ 1), 1)
1507       ENDIF
1508       CALL ros_solve(ghimj, pivot, k(ioffset+ 1))
1509
1510   END DO stage
1511
1512
1513!~~~>  compute the new solution
1514   !slim: CALL wcopy(n, y, 1, ynew, 1)
1515   ynew(1:n) = y(1:n)
1516   DO j=1, ros_s
1517         CALL waxpy(n, ros_m(j), k(n* (j- 1) + 1), 1, ynew, 1)
1518   ENDDO
1519
1520!~~~>  compute the error estimation
1521   !slim: CALL wscal(n, zero, yerr, 1)
1522   yerr(1:n) = zero
1523   DO j=1, ros_s
1524        CALL waxpy(n, ros_e(j), k(n* (j- 1) + 1), 1, yerr, 1)
1525   ENDDO
1526   err = ros_errornorm(y, ynew, yerr, abstol, reltol, vectortol)
1527
1528!~~~> new step size is bounded by facmin <= hnew/h <= facmax
1529   fac  = min(facmax, max(facmin, facsafe/err** (one/ros_elo)))
1530   hnew = h* fac
1531
1532!~~~>  check the error magnitude and adjust step size
1533   istatus(nstp) = istatus(nstp) + 1
1534   IF ((err <= one).or.(h <= hmin))THEN  !~~~> accept step
1535      istatus(nacc) = istatus(nacc) + 1
1536      !slim: CALL wcopy(n, ynew, 1, y, 1)
1537      y(1:n) = ynew(1:n)
1538      t = t + direction* h
1539      hnew = max(hmin, min(hnew, hmax))
1540      IF (rejectlasth)THEN  ! no step size increase after a rejected step
1541         hnew = min(hnew, h)
1542      ENDIF
1543      rstatus(nhexit) = h
1544      rstatus(nhnew) = hnew
1545      rstatus(ntexit) = t
1546      rejectlasth = .FALSE.
1547      rejectmoreh = .FALSE.
1548      h = hnew
1549      exit untilaccepted ! exit the loop: WHILE step not accepted
1550   ELSE           !~~~> reject step
1551      IF (rejectmoreh)THEN
1552         hnew = h* facrej
1553      ENDIF
1554      rejectmoreh = rejectlasth
1555      rejectlasth = .TRUE.
1556      h = hnew
1557      IF (istatus(nacc)>= 1) istatus(nrej) = istatus(nrej) + 1
1558   ENDIF ! err <= 1
1559
1560   END DO untilaccepted
1561
1562   END DO timeloop
1563
1564!~~~> succesful exit
1565   ierr = 1  !~~~> the integration was successful
1566
1567  END SUBROUTINE ros_integrator
1568
1569
1570!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1571  REAL(kind=dp)FUNCTION ros_errornorm(y, ynew, yerr, &
1572                               abstol, reltol, vectortol)
1573!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1574!~~~> computes the "scaled norm" of the error vector yerr
1575!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1576
1577! Input arguments
1578   REAL(kind=dp), INTENT(IN):: y(n), ynew(n), &
1579          yerr(n), abstol(n), reltol(n)
1580   LOGICAL, INTENT(IN)::  vectortol
1581! Local variables
1582   REAL(kind=dp):: err, scale, ymax
1583   INTEGER  :: i
1584   REAL(kind=dp), PARAMETER :: zero = 0.0_dp
1585
1586   err = zero
1587   DO i=1, n
1588     ymax = max(abs(y(i)), abs(ynew(i)))
1589     IF (vectortol)THEN
1590       scale = abstol(i) + reltol(i) * ymax
1591     ELSE
1592       scale = abstol(1) + reltol(1) * ymax
1593     ENDIF
1594     err = err+ (yerr(i) /scale) ** 2
1595   ENDDO
1596   err  = sqrt(err/n)
1597
1598   ros_errornorm = max(err, 1.0d-10)
1599
1600  END FUNCTION ros_errornorm
1601
1602
1603!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1604  SUBROUTINE ros_funtimederivative(t, roundoff, y, &
1605                fcn0, dfdt)
1606!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1607!~~~> the time partial derivative of the FUNCTION by finite differences
1608!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1609
1610!~~~> input arguments
1611   REAL(kind=dp), INTENT(IN):: t, roundoff, y(n), fcn0(n)
1612!~~~> output arguments
1613   REAL(kind=dp), INTENT(OUT):: dfdt(n)
1614!~~~> local variables
1615   REAL(kind=dp):: delta
1616   REAL(kind=dp), PARAMETER :: one = 1.0_dp, deltamin = 1.0e-6_dp
1617
1618   delta = sqrt(roundoff) * max(deltamin, abs(t))
1619   CALL funtemplate(t+ delta, y, dfdt)
1620   istatus(nfun) = istatus(nfun) + 1
1621   CALL waxpy(n, (- one), fcn0, 1, dfdt, 1)
1622   CALL wscal(n, (one/delta), dfdt, 1)
1623
1624  END SUBROUTINE ros_funtimederivative
1625
1626
1627!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1628  SUBROUTINE ros_preparematrix(h, direction, gam, &
1629             jac0, ghimj, pivot, singular)
1630! --- --- --- --- --- --- --- --- --- --- --- --- ---
1631!  Prepares the LHS matrix for stage calculations
1632!  1.  Construct Ghimj = 1/(H*ham) - Jac0
1633!      "(Gamma H) Inverse Minus Jacobian"
1634!  2.  Repeat LU decomposition of Ghimj until successful.
1635!       -half the step size if LU decomposition fails and retry
1636!       -exit after 5 consecutive fails
1637! --- --- --- --- --- --- --- --- --- --- --- --- ---
1638
1639!~~~> input arguments
1640#ifdef full_algebra   
1641   REAL(kind=dp), INTENT(IN)::  jac0(n, n)
1642#else
1643   REAL(kind=dp), INTENT(IN)::  jac0(lu_nonzero)
1644#endif   
1645   REAL(kind=dp), INTENT(IN)::  gam
1646   INTEGER, INTENT(IN)::  direction
1647!~~~> output arguments
1648#ifdef full_algebra   
1649   REAL(kind=dp), INTENT(OUT):: ghimj(n, n)
1650#else
1651   REAL(kind=dp), INTENT(OUT):: ghimj(lu_nonzero)
1652#endif   
1653   LOGICAL, INTENT(OUT)::  singular
1654   INTEGER, INTENT(OUT)::  pivot(n)
1655!~~~> inout arguments
1656   REAL(kind=dp), INTENT(INOUT):: h   ! step size is decreased when lu fails
1657!~~~> local variables
1658   INTEGER  :: i, ising, nconsecutive
1659   REAL(kind=dp):: ghinv
1660   REAL(kind=dp), PARAMETER :: one  = 1.0_dp, half = 0.5_dp
1661
1662   nconsecutive = 0
1663   singular = .TRUE.
1664
1665   DO WHILE (singular)
1666
1667!~~~>    construct ghimj = 1/(h* gam) - jac0
1668#ifdef full_algebra   
1669     !slim: CALL wcopy(n* n, jac0, 1, ghimj, 1)
1670     !slim: CALL wscal(n* n, (- one), ghimj, 1)
1671     ghimj = - jac0
1672     ghinv = one/(direction* h* gam)
1673     DO i=1, n
1674       ghimj(i, i) = ghimj(i, i) + ghinv
1675     ENDDO
1676#else
1677     !slim: CALL wcopy(lu_nonzero, jac0, 1, ghimj, 1)
1678     !slim: CALL wscal(lu_nonzero, (- one), ghimj, 1)
1679     ghimj(1:lu_nonzero) = - jac0(1:lu_nonzero)
1680     ghinv = one/(direction* h* gam)
1681     DO i=1, n
1682       ghimj(lu_diag(i)) = ghimj(lu_diag(i)) + ghinv
1683     ENDDO
1684#endif   
1685!~~~>    compute lu decomposition
1686     CALL ros_decomp( ghimj, pivot, ising)
1687     IF (ising == 0)THEN
1688!~~~>    IF successful done
1689        singular = .FALSE.
1690     ELSE ! ising .ne. 0
1691!~~~>    IF unsuccessful half the step size; IF 5 consecutive fails THEN RETURN
1692        istatus(nsng) = istatus(nsng) + 1
1693        nconsecutive = nconsecutive+1
1694        singular = .TRUE.
1695        PRINT*,'Warning: LU Decomposition returned ISING = ',ISING
1696        IF (nconsecutive <= 5)THEN ! less than 5 consecutive failed decompositions
1697           h = h* half
1698        ELSE  ! more than 5 consecutive failed decompositions
1699           RETURN
1700        ENDIF  ! nconsecutive
1701      ENDIF    ! ising
1702
1703   END DO ! WHILE singular
1704
1705  END SUBROUTINE ros_preparematrix
1706
1707
1708!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1709  SUBROUTINE ros_decomp( a, pivot, ising)
1710!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1711!  Template for the LU decomposition
1712!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1713!~~~> inout variables
1714#ifdef full_algebra   
1715   REAL(kind=dp), INTENT(INOUT):: a(n, n)
1716#else   
1717   REAL(kind=dp), INTENT(INOUT):: a(lu_nonzero)
1718#endif
1719!~~~> output variables
1720   INTEGER, INTENT(OUT):: pivot(n), ising
1721
1722#ifdef full_algebra   
1723   CALL  dgetrf( n, n, a, n, pivot, ising)
1724#else   
1725   CALL kppdecomp(a, ising)
1726   pivot(1) = 1
1727#endif
1728   istatus(ndec) = istatus(ndec) + 1
1729
1730  END SUBROUTINE ros_decomp
1731
1732
1733!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1734  SUBROUTINE ros_solve( a, pivot, b)
1735!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1736!  Template for the forward/backward substitution (using pre-computed LU decomposition)
1737!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1738!~~~> input variables
1739#ifdef full_algebra   
1740   REAL(kind=dp), INTENT(IN):: a(n, n)
1741   INTEGER :: ising
1742#else   
1743   REAL(kind=dp), INTENT(IN):: a(lu_nonzero)
1744#endif
1745   INTEGER, INTENT(IN):: pivot(n)
1746!~~~> inout variables
1747   REAL(kind=dp), INTENT(INOUT):: b(n)
1748
1749#ifdef full_algebra   
1750   CALL  DGETRS( 'N',N ,1,A,N,Pivot,b,N,ISING)
1751   IF (info < 0)THEN
1752      print* , "error in dgetrs. ising=", ising
1753   ENDIF 
1754#else   
1755   CALL kppsolve( a, b)
1756#endif
1757
1758   istatus(nsol) = istatus(nsol) + 1
1759
1760  END SUBROUTINE ros_solve
1761
1762
1763
1764!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1765  SUBROUTINE ros2
1766!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1767! --- AN L-STABLE METHOD,2 stages,order 2
1768!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1769
1770   double precision g
1771
1772    g = 1.0_dp + 1.0_dp/sqrt(2.0_dp)
1773    rosmethod = rs2
1774!~~~> name of the method
1775    ros_Name = 'ROS-2'
1776!~~~> number of stages
1777    ros_s = 2
1778
1779!~~~> the coefficient matrices a and c are strictly lower triangular.
1780!   The lower triangular (subdiagonal) elements are stored in row-wise order:
1781!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
1782!   The general mapping formula is:
1783!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
1784!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
1785
1786    ros_a(1) = (1.0_dp) /g
1787    ros_c(1) = (- 2.0_dp) /g
1788!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
1789!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
1790    ros_newf(1) = .TRUE.
1791    ros_newf(2) = .TRUE.
1792!~~~> m_i = coefficients for new step solution
1793    ros_m(1) = (3.0_dp) /(2.0_dp* g)
1794    ros_m(2) = (1.0_dp) /(2.0_dp* g)
1795! E_i = Coefficients for error estimator
1796    ros_e(1) = 1.0_dp/(2.0_dp* g)
1797    ros_e(2) = 1.0_dp/(2.0_dp* g)
1798!~~~> ros_elo = estimator of local order - the minimum between the
1799!    main and the embedded scheme orders plus one
1800    ros_elo = 2.0_dp
1801!~~~> y_stage_i ~ y( t + h* alpha_i)
1802    ros_alpha(1) = 0.0_dp
1803    ros_alpha(2) = 1.0_dp
1804!~~~> gamma_i = \sum_j  gamma_{i, j}
1805    ros_gamma(1) = g
1806    ros_gamma(2) = -g
1807
1808 END SUBROUTINE ros2
1809
1810
1811!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1812  SUBROUTINE ros3
1813!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1814! --- AN L-STABLE METHOD,3 stages,order 3,2 function evaluations
1815!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1816
1817   rosmethod = rs3
1818!~~~> name of the method
1819   ros_Name = 'ROS-3'
1820!~~~> number of stages
1821   ros_s = 3
1822
1823!~~~> the coefficient matrices a and c are strictly lower triangular.
1824!   The lower triangular (subdiagonal) elements are stored in row-wise order:
1825!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
1826!   The general mapping formula is:
1827!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
1828!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
1829
1830   ros_a(1) = 1.0_dp
1831   ros_a(2) = 1.0_dp
1832   ros_a(3) = 0.0_dp
1833
1834   ros_c(1) = - 0.10156171083877702091975600115545e+01_dp
1835   ros_c(2) =  0.40759956452537699824805835358067e+01_dp
1836   ros_c(3) =  0.92076794298330791242156818474003e+01_dp
1837!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
1838!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
1839   ros_newf(1) = .TRUE.
1840   ros_newf(2) = .TRUE.
1841   ros_newf(3) = .FALSE.
1842!~~~> m_i = coefficients for new step solution
1843   ros_m(1) =  0.1e+01_dp
1844   ros_m(2) =  0.61697947043828245592553615689730e+01_dp
1845   ros_m(3) = - 0.42772256543218573326238373806514_dp
1846! E_i = Coefficients for error estimator
1847   ros_e(1) =  0.5_dp
1848   ros_e(2) = - 0.29079558716805469821718236208017e+01_dp
1849   ros_e(3) =  0.22354069897811569627360909276199_dp
1850!~~~> ros_elo = estimator of local order - the minimum between the
1851!    main and the embedded scheme orders plus 1
1852   ros_elo = 3.0_dp
1853!~~~> y_stage_i ~ y( t + h* alpha_i)
1854   ros_alpha(1) = 0.0_dp
1855   ros_alpha(2) = 0.43586652150845899941601945119356_dp
1856   ros_alpha(3) = 0.43586652150845899941601945119356_dp
1857!~~~> gamma_i = \sum_j  gamma_{i, j}
1858   ros_gamma(1) = 0.43586652150845899941601945119356_dp
1859   ros_gamma(2) = 0.24291996454816804366592249683314_dp
1860   ros_gamma(3) = 0.21851380027664058511513169485832e+01_dp
1861
1862  END SUBROUTINE ros3
1863
1864!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1865
1866
1867!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1868  SUBROUTINE ros4
1869!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1870!     L-STABLE ROSENBROCK METHOD OF ORDER 4,WITH 4 STAGES
1871!     L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3
1872!
1873!      E. HAIRER AND G. WANNER,SOLVING ORDINARY DIFFERENTIAL
1874!      EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
1875!      SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
1876!      SPRINGER-VERLAG (1990)
1877!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1878
1879
1880   rosmethod = rs4
1881!~~~> name of the method
1882   ros_Name = 'ROS-4'
1883!~~~> number of stages
1884   ros_s = 4
1885
1886!~~~> the coefficient matrices a and c are strictly lower triangular.
1887!   The lower triangular (subdiagonal) elements are stored in row-wise order:
1888!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
1889!   The general mapping formula is:
1890!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
1891!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
1892
1893   ros_a(1) = 0.2000000000000000e+01_dp
1894   ros_a(2) = 0.1867943637803922e+01_dp
1895   ros_a(3) = 0.2344449711399156_dp
1896   ros_a(4) = ros_a(2)
1897   ros_a(5) = ros_a(3)
1898   ros_a(6) = 0.0_dp
1899
1900   ros_c(1) = -0.7137615036412310e+01_dp
1901   ros_c(2) = 0.2580708087951457e+01_dp
1902   ros_c(3) = 0.6515950076447975_dp
1903   ros_c(4) = -0.2137148994382534e+01_dp
1904   ros_c(5) = -0.3214669691237626_dp
1905   ros_c(6) = -0.6949742501781779_dp
1906!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
1907!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
1908   ros_newf(1) = .TRUE.
1909   ros_newf(2) = .TRUE.
1910   ros_newf(3) = .TRUE.
1911   ros_newf(4) = .FALSE.
1912!~~~> m_i = coefficients for new step solution
1913   ros_m(1) = 0.2255570073418735e+01_dp
1914   ros_m(2) = 0.2870493262186792_dp
1915   ros_m(3) = 0.4353179431840180_dp
1916   ros_m(4) = 0.1093502252409163e+01_dp
1917!~~~> e_i  = coefficients for error estimator
1918   ros_e(1) = -0.2815431932141155_dp
1919   ros_e(2) = -0.7276199124938920e-01_dp
1920   ros_e(3) = -0.1082196201495311_dp
1921   ros_e(4) = -0.1093502252409163e+01_dp
1922!~~~> ros_elo  = estimator of local order - the minimum between the
1923!    main and the embedded scheme orders plus 1
1924   ros_elo  = 4.0_dp
1925!~~~> y_stage_i ~ y( t + h* alpha_i)
1926   ros_alpha(1) = 0.0_dp
1927   ros_alpha(2) = 0.1145640000000000e+01_dp
1928   ros_alpha(3) = 0.6552168638155900_dp
1929   ros_alpha(4) = ros_alpha(3)
1930!~~~> gamma_i = \sum_j  gamma_{i, j}
1931   ros_gamma(1) = 0.5728200000000000_dp
1932   ros_gamma(2) = -0.1769193891319233e+01_dp
1933   ros_gamma(3) = 0.7592633437920482_dp
1934   ros_gamma(4) = -0.1049021087100450_dp
1935
1936  END SUBROUTINE ros4
1937
1938!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1939  SUBROUTINE rodas3
1940!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1941! --- A STIFFLY-STABLE METHOD,4 stages,order 3
1942!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1943
1944
1945   rosmethod = rd3
1946!~~~> name of the method
1947   ros_Name = 'RODAS-3'
1948!~~~> number of stages
1949   ros_s = 4
1950
1951!~~~> the coefficient matrices a and c are strictly lower triangular.
1952!   The lower triangular (subdiagonal) elements are stored in row-wise order:
1953!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
1954!   The general mapping formula is:
1955!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
1956!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
1957
1958   ros_a(1) = 0.0_dp
1959   ros_a(2) = 2.0_dp
1960   ros_a(3) = 0.0_dp
1961   ros_a(4) = 2.0_dp
1962   ros_a(5) = 0.0_dp
1963   ros_a(6) = 1.0_dp
1964
1965   ros_c(1) = 4.0_dp
1966   ros_c(2) = 1.0_dp
1967   ros_c(3) = -1.0_dp
1968   ros_c(4) = 1.0_dp
1969   ros_c(5) = -1.0_dp
1970   ros_c(6) = -(8.0_dp/3.0_dp)
1971
1972!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
1973!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
1974   ros_newf(1) = .TRUE.
1975   ros_newf(2) = .FALSE.
1976   ros_newf(3) = .TRUE.
1977   ros_newf(4) = .TRUE.
1978!~~~> m_i = coefficients for new step solution
1979   ros_m(1) = 2.0_dp
1980   ros_m(2) = 0.0_dp
1981   ros_m(3) = 1.0_dp
1982   ros_m(4) = 1.0_dp
1983!~~~> e_i  = coefficients for error estimator
1984   ros_e(1) = 0.0_dp
1985   ros_e(2) = 0.0_dp
1986   ros_e(3) = 0.0_dp
1987   ros_e(4) = 1.0_dp
1988!~~~> ros_elo  = estimator of local order - the minimum between the
1989!    main and the embedded scheme orders plus 1
1990   ros_elo  = 3.0_dp
1991!~~~> y_stage_i ~ y( t + h* alpha_i)
1992   ros_alpha(1) = 0.0_dp
1993   ros_alpha(2) = 0.0_dp
1994   ros_alpha(3) = 1.0_dp
1995   ros_alpha(4) = 1.0_dp
1996!~~~> gamma_i = \sum_j  gamma_{i, j}
1997   ros_gamma(1) = 0.5_dp
1998   ros_gamma(2) = 1.5_dp
1999   ros_gamma(3) = 0.0_dp
2000   ros_gamma(4) = 0.0_dp
2001
2002  END SUBROUTINE rodas3
2003
2004!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2005  SUBROUTINE rodas4
2006!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2007!     STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4,WITH 6 STAGES
2008!
2009!      E. HAIRER AND G. WANNER,SOLVING ORDINARY DIFFERENTIAL
2010!      EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
2011!      SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
2012!      SPRINGER-VERLAG (1996)
2013!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2014
2015
2016    rosmethod = rd4
2017!~~~> name of the method
2018    ros_Name = 'RODAS-4'
2019!~~~> number of stages
2020    ros_s = 6
2021
2022!~~~> y_stage_i ~ y( t + h* alpha_i)
2023    ros_alpha(1) = 0.000_dp
2024    ros_alpha(2) = 0.386_dp
2025    ros_alpha(3) = 0.210_dp
2026    ros_alpha(4) = 0.630_dp
2027    ros_alpha(5) = 1.000_dp
2028    ros_alpha(6) = 1.000_dp
2029
2030!~~~> gamma_i = \sum_j  gamma_{i, j}
2031    ros_gamma(1) = 0.2500000000000000_dp
2032    ros_gamma(2) = -0.1043000000000000_dp
2033    ros_gamma(3) = 0.1035000000000000_dp
2034    ros_gamma(4) = -0.3620000000000023e-01_dp
2035    ros_gamma(5) = 0.0_dp
2036    ros_gamma(6) = 0.0_dp
2037
2038!~~~> the coefficient matrices a and c are strictly lower triangular.
2039!   The lower triangular (subdiagonal) elements are stored in row-wise order:
2040!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
2041!   The general mapping formula is:  A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
2042!                  C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
2043
2044    ros_a(1) = 0.1544000000000000e+01_dp
2045    ros_a(2) = 0.9466785280815826_dp
2046    ros_a(3) = 0.2557011698983284_dp
2047    ros_a(4) = 0.3314825187068521e+01_dp
2048    ros_a(5) = 0.2896124015972201e+01_dp
2049    ros_a(6) = 0.9986419139977817_dp
2050    ros_a(7) = 0.1221224509226641e+01_dp
2051    ros_a(8) = 0.6019134481288629e+01_dp
2052    ros_a(9) = 0.1253708332932087e+02_dp
2053    ros_a(10) = -0.6878860361058950_dp
2054    ros_a(11) = ros_a(7)
2055    ros_a(12) = ros_a(8)
2056    ros_a(13) = ros_a(9)
2057    ros_a(14) = ros_a(10)
2058    ros_a(15) = 1.0_dp
2059
2060    ros_c(1) = -0.5668800000000000e+01_dp
2061    ros_c(2) = -0.2430093356833875e+01_dp
2062    ros_c(3) = -0.2063599157091915_dp
2063    ros_c(4) = -0.1073529058151375_dp
2064    ros_c(5) = -0.9594562251023355e+01_dp
2065    ros_c(6) = -0.2047028614809616e+02_dp
2066    ros_c(7) = 0.7496443313967647e+01_dp
2067    ros_c(8) = -0.1024680431464352e+02_dp
2068    ros_c(9) = -0.3399990352819905e+02_dp
2069    ros_c(10) = 0.1170890893206160e+02_dp
2070    ros_c(11) = 0.8083246795921522e+01_dp
2071    ros_c(12) = -0.7981132988064893e+01_dp
2072    ros_c(13) = -0.3152159432874371e+02_dp
2073    ros_c(14) = 0.1631930543123136e+02_dp
2074    ros_c(15) = -0.6058818238834054e+01_dp
2075
2076!~~~> m_i = coefficients for new step solution
2077    ros_m(1) = ros_a(7)
2078    ros_m(2) = ros_a(8)
2079    ros_m(3) = ros_a(9)
2080    ros_m(4) = ros_a(10)
2081    ros_m(5) = 1.0_dp
2082    ros_m(6) = 1.0_dp
2083
2084!~~~> e_i  = coefficients for error estimator
2085    ros_e(1) = 0.0_dp
2086    ros_e(2) = 0.0_dp
2087    ros_e(3) = 0.0_dp
2088    ros_e(4) = 0.0_dp
2089    ros_e(5) = 0.0_dp
2090    ros_e(6) = 1.0_dp
2091
2092!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
2093!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
2094    ros_newf(1) = .TRUE.
2095    ros_newf(2) = .TRUE.
2096    ros_newf(3) = .TRUE.
2097    ros_newf(4) = .TRUE.
2098    ros_newf(5) = .TRUE.
2099    ros_newf(6) = .TRUE.
2100
2101!~~~> ros_elo  = estimator of local order - the minimum between the
2102!        main and the embedded scheme orders plus 1
2103    ros_elo = 4.0_dp
2104
2105  END SUBROUTINE rodas4
2106!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2107  SUBROUTINE rang3
2108!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2109! STIFFLY-STABLE W METHOD OF ORDER 3,WITH 4 STAGES
2110!
2111! J. RANG and L. ANGERMANN
2112! NEW ROSENBROCK W-METHODS OF ORDER 3
2113! FOR PARTIAL DIFFERENTIAL ALGEBRAIC
2114!        EQUATIONS OF INDEX 1                                             
2115! BIT Numerical Mathematics (2005) 45: 761-787
2116!  DOI: 10.1007/s10543-005-0035-y
2117! Table 4.1-4.2
2118!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2119
2120
2121    rosmethod = rg3
2122!~~~> name of the method
2123    ros_Name = 'RANG-3'
2124!~~~> number of stages
2125    ros_s = 4
2126
2127    ros_a(1) = 5.09052051067020d+00;
2128    ros_a(2) = 5.09052051067020d+00;
2129    ros_a(3) = 0.0d0;
2130    ros_a(4) = 4.97628111010787d+00;
2131    ros_a(5) = 2.77268164715849d-02;
2132    ros_a(6) = 2.29428036027904d-01;
2133
2134    ros_c(1) = - 1.16790812312283d+01;
2135    ros_c(2) = - 1.64057326467367d+01;
2136    ros_c(3) = - 2.77268164715850d-01;
2137    ros_c(4) = - 8.38103960500476d+00;
2138    ros_c(5) = - 8.48328409199343d-01;
2139    ros_c(6) =  2.87009860433106d-01;
2140
2141    ros_m(1) =  5.22582761233094d+00;
2142    ros_m(2) = - 5.56971148154165d-01;
2143    ros_m(3) =  3.57979469353645d-01;
2144    ros_m(4) =  1.72337398521064d+00;
2145
2146    ros_e(1) = - 5.16845212784040d+00;
2147    ros_e(2) = - 1.26351942603842d+00;
2148    ros_e(3) = - 1.11022302462516d-16;
2149    ros_e(4) =  2.22044604925031d-16;
2150
2151    ros_alpha(1) = 0.0d00;
2152    ros_alpha(2) = 2.21878746765329d+00;
2153    ros_alpha(3) = 2.21878746765329d+00;
2154    ros_alpha(4) = 1.55392337535788d+00;
2155
2156    ros_gamma(1) =  4.35866521508459d-01;
2157    ros_gamma(2) = - 1.78292094614483d+00;
2158    ros_gamma(3) = - 2.46541900496934d+00;
2159    ros_gamma(4) = - 8.05529997906370d-01;
2160
2161
2162!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
2163!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
2164    ros_newf(1) = .TRUE.
2165    ros_newf(2) = .TRUE.
2166    ros_newf(3) = .TRUE.
2167    ros_newf(4) = .TRUE.
2168
2169!~~~> ros_elo  = estimator of local order - the minimum between the
2170!        main and the embedded scheme orders plus 1
2171    ros_elo = 3.0_dp
2172
2173  END SUBROUTINE rang3
2174!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2175
2176!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2177!   End of the set of internal Rosenbrock subroutines
2178!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2179END SUBROUTINE rosenbrock
2180 
2181SUBROUTINE funtemplate( t, y, ydot)
2182!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2183!  Template for the ODE function call.
2184!  Updates the rate coefficients (and possibly the fixed species) at each call
2185!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2186!~~~> input variables
2187   REAL(kind=dp):: t, y(nvar)
2188!~~~> output variables
2189   REAL(kind=dp):: ydot(nvar)
2190!~~~> local variables
2191   REAL(kind=dp):: told
2192
2193   told = time
2194   time = t
2195   CALL fun( y, fix, rconst, ydot)
2196   time = told
2197
2198END SUBROUTINE funtemplate
2199 
2200SUBROUTINE jactemplate( t, y, jcb)
2201!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2202!  Template for the ODE Jacobian call.
2203!  Updates the rate coefficients (and possibly the fixed species) at each call
2204!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2205!~~~> input variables
2206    REAL(kind=dp):: t, y(nvar)
2207!~~~> output variables
2208#ifdef full_algebra   
2209    REAL(kind=dp):: jv(lu_nonzero), jcb(nvar, nvar)
2210#else
2211    REAL(kind=dp):: jcb(lu_nonzero)
2212#endif   
2213!~~~> local variables
2214    REAL(kind=dp):: told
2215#ifdef full_algebra   
2216    INTEGER :: i, j
2217#endif   
2218
2219    told = time
2220    time = t
2221#ifdef full_algebra   
2222    CALL jac_sp(y, fix, rconst, jv)
2223    DO j=1, nvar
2224      DO i=1, nvar
2225         jcb(i, j) = 0.0_dp
2226      ENDDO
2227    ENDDO
2228    DO i=1, lu_nonzero
2229       jcb(lu_irow(i), lu_icol(i)) = jv(i)
2230    ENDDO
2231#else
2232    CALL jac_sp( y, fix, rconst, jcb)
2233#endif   
2234    time = told
2235
2236END SUBROUTINE jactemplate
2237 
2238  SUBROUTINE kppdecomp( jvs, ier)                                   
2239  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2240  !        sparse lu factorization                                   
2241  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2242  ! loop expansion generated by kp4                                   
2243                                                                     
2244    INTEGER  :: ier                                                   
2245    REAL(kind=dp):: jvs(lu_nonzero), w(nvar), a             
2246    INTEGER  :: k, kk, j, jj                                         
2247                                                                     
2248    a = 0.                                                           
2249    ier = 0                                                           
2250                                                                     
2251!   i = 1
2252!   i = 2
2253!   i = 3
2254    jvs(5) = (jvs(5)) / jvs(2) 
2255    jvs(6) = jvs(6) - jvs(3) * jvs(5)
2256    jvs(7) = jvs(7) - jvs(4) * jvs(5)
2257!   i = 4
2258    jvs(8) = (jvs(8)) / jvs(2) 
2259    a = 0.0; a = a  - jvs(3) * jvs(8)
2260    jvs(9) = (jvs(9) + a) / jvs(6) 
2261    jvs(10) = jvs(10) - jvs(4) * jvs(8) - jvs(7) * jvs(9)
2262    RETURN                                                           
2263                                                                     
2264  END SUBROUTINE kppdecomp                                           
2265 
2266SUBROUTINE chem_gasphase_integrate (time_step_len, conc, tempi, qvapi, fakti, photo, ierrf, xnacc, xnrej, istatus, l_debug, pe, &
2267                     icntrl_i, rcntrl_i) 
2268                                                                   
2269  IMPLICIT NONE                                                     
2270                                                                   
2271  REAL(dp), INTENT(IN)                  :: time_step_len           
2272  REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: conc                   
2273  REAL(dp), DIMENSION(:, :), INTENT(IN)   :: photo                   
2274  REAL(dp), DIMENSION(:), INTENT(IN)     :: tempi                   
2275  REAL(dp), DIMENSION(:), INTENT(IN)     :: qvapi                   
2276  REAL(dp), DIMENSION(:), INTENT(IN)     :: fakti                   
2277  INTEGER, INTENT(OUT), OPTIONAL        :: ierrf(:)               
2278  INTEGER, INTENT(OUT), OPTIONAL        :: xnacc(:)               
2279  INTEGER, INTENT(OUT), OPTIONAL        :: xnrej(:)               
2280  INTEGER, INTENT(INOUT), OPTIONAL      :: istatus(:)             
2281  INTEGER, INTENT(IN), OPTIONAL         :: pe                     
2282  LOGICAL, INTENT(IN), OPTIONAL         :: l_debug                 
2283  INTEGER, DIMENSION(nkppctrl), INTENT(IN), OPTIONAL  :: icntrl_i         
2284  REAL(dp), DIMENSION(nkppctrl), INTENT(IN), OPTIONAL  :: rcntrl_i         
2285                                                                   
2286  INTEGER                                 :: k   ! loop variable     
2287  REAL(dp)                               :: dt                     
2288  INTEGER, DIMENSION(20)                :: istatus_u               
2289  INTEGER                                :: ierr_u                 
2290  INTEGER                                :: istatf                 
2291  INTEGER                                :: vl_dim_lo               
2292                                                                   
2293                                                                   
2294  IF (PRESENT (istatus)) istatus = 0                             
2295  IF (PRESENT (icntrl_i)) icntrl  = icntrl_i                     
2296  IF (PRESENT (rcntrl_i)) rcntrl  = rcntrl_i                     
2297                                                                   
2298  vl_glo = size(tempi, 1)                                           
2299                                                                   
2300  vl_dim_lo = vl_dim                                               
2301  DO k=1, vl_glo, vl_dim_lo                                           
2302    is = k                                                         
2303    ie = min(k+ vl_dim_lo-1, vl_glo)                                 
2304    vl = ie-is+ 1                                                   
2305                                                                   
2306    c(:) = conc(is, :)                                           
2307                                                                   
2308    temp = tempi(is)                                             
2309                                                                   
2310    qvap = qvapi(is)                                             
2311                                                                   
2312    fakt = fakti(is)                                             
2313                                                                   
2314    CALL initialize                                                 
2315                                                                   
2316    phot(:) = photo(is, :)                                           
2317                                                                   
2318    CALL update_rconst                                             
2319                                                                   
2320    dt = time_step_len                                             
2321                                                                   
2322    ! integrate from t=0 to t=dt                                   
2323    CALL integrate(0._dp, dt, icntrl, rcntrl, istatus_u = istatus_u, ierr_u=ierr_u)
2324                                                                   
2325                                                                   
2326   IF (PRESENT(l_debug) .AND. PRESENT(pe)) THEN                       
2327      IF (l_debug) CALL error_output(conc(is, :), ierr_u, pe)         
2328   ENDIF                                                             
2329                                                                     
2330    conc(is, :) = c(:)                                               
2331                                                                   
2332    ! RETURN diagnostic information                                 
2333                                                                   
2334    IF (PRESENT(ierrf))   ierrf(is) = ierr_u                     
2335    IF (PRESENT(xnacc))   xnacc(is) = istatus_u(4)               
2336    IF (PRESENT(xnrej))   xnrej(is) = istatus_u(5)               
2337                                                                   
2338    IF (PRESENT (istatus)) THEN                                   
2339      istatus(1:8) = istatus(1:8) + istatus_u(1:8)               
2340    ENDIF                                                         
2341                                                                   
2342  END DO                                                           
2343 
2344                                                                   
2345! Deallocate input arrays                                           
2346                                                                   
2347                                                                   
2348  data_loaded = .FALSE.                                             
2349                                                                   
2350  RETURN                                                           
2351END SUBROUTINE chem_gasphase_integrate                                       
2352
2353END MODULE chem_gasphase_mod
2354 
Note: See TracBrowser for help on using the repository browser.