source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_passive1/chem_gasphase_mod.f90 @ 3698

Last change on this file since 3698 was 3698, checked in by suehring, 5 years ago

Fix bad commit in gasphase_preproc

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