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

Last change on this file since 3800 was 3800, checked in by forkel, 3 years ago

added leading blanks to dummy statements

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