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

Last change on this file since 3458 was 3458, checked in by kanani, 3 years ago

Reintegrated fixes/changes from branch chemistry

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