source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_salsa+phstat/chem_gasphase_mod.f90 @ 3681

Last change on this file since 3681 was 3681, checked in by hellstea, 3 years ago

Major update of pmc_interface_mod

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