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

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

Fix bad commit in gasphase_preproc

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